Multivariate And Multidimensional Quality Gain-Loss Function and Its Applications Based on Nonseparable Gaussian Processes

1. Introduction

In the 1960s, Taguchi [1] associated product quality with economic loss and proposed that the quality loss function attains its minimum value of zero when quality characteristics exactly match their target values. Any deviation of the quality characteristic value from the target value will result in quality loss. Spiring [23] proposed an asymmetric quality loss function to capture the asymmetric nature of losses in production processes. To address the unboundedness of the loss function, he further developed an inverted normal quality loss function. Building upon Taguchi's quality loss function, M Abolmohammadi et al. [4] extended the loss function by incorporating linear, quadratic, exponential, and Linex functions, thereby establishing the Linex predicted loss function. Artiles-leon [5] proposed a dimensionless “standardized” multi-parameter quality loss function. Pignatiello [6] introduced a multivariate quality loss function to address scenarios where the overall quality of products is jointly influenced by multiple quality characteristics. Error! Reference source not found. developed Type I and Type II exponential quality loss function models which is more suitable for the real scene, accounting for the practical tendency to replace rather than repair scrapped products and the asymmetric nature of quality losses. Wang et al. [7] optimized the quality cost curve by incorporating the Cobb-Douglas function into the Taguchi quality loss function. Cao et al. [8] developed a fuzzy quality loss function model addressing the fuzzy nature of product quality. Wang et al. [9] established a robust parameter optimization model for multiple characteristics based on Bayesian posterior probability theory. Ni et al. [11] conducted an in-depth research on process mean estimation under asymmetric loss functions.

To address mutual compensation and adaptive coupling phenomena between processes in engineering construction—where subsequent processes compensate for losses induced by preceding ones, and even processes within the same sequence complement each other to reduce overall quality loss, thereby producing a “quality compensation” effect.Wang et al. [13] defined the constant term in the Taylor series expansion as quality compensation and introduced the concept of a quality gain-loss function. Assuming constant compensation while retaining linear loss terms, they developed “larger the better” and “smaller the better” quality gain-loss function models. [14] Subsequently, fuzzy quality gain-loss function models [15,16] and grey quality gain-loss function models [17,18] were proposed based on fuzzy theory and grey system theory. To address the boundedness of quality gains and losses, an inverted normal quality gain-loss function was proposed, leading to the construction of a multivariate inverted normal model [19]. By retaining the third-order Taylor expansion, cubic and asymmetric cubic quality gain-loss functions were proposed [20,21]. To address asymmetry in quality gain-loss functions, a quadratic exponential quality gain-loss function model was developed [22]. To correct loss bias arising from omitting the cubic term, a cubic-term loss-considering “larger the better” and “smaller the better” quality gain-loss function model was proposed [23]. In engineering projects, quality losses are typically caused by the combined effects of multiple quality characteristics. Accordingly, a wide range of multivariate modeling and analysis methods are available, including Bayesian networks, multiple linear regression, random forests, nearest neighbor methods, and Gaussian processes.

Non-separable Gaussian Processes (NSGP) is an extension of Gaussian Processes (GP) designed to model complex spatial and temporal data where the correlation between spatial and temporal factors is relatively complex.Fricker Thomas E et al. [24] highlighted the flexibility and predictive strength of non-separable covariance structures; Ilias Bilionis et al. [25] proposed a posteriori probability measure for potential surrogate function spaces defined by Gaussian Process models; Datta Abhirup et al. [26] developed a scalable dynamic k-nearest neighbor Gaussian process model capable of providing sparse approximations for arbitrary spatiotemporal GPs. Dionissios T. et al. [27] employed a hybrid spectral approach to develop a new non-separable covariance kernel, validating it under various damping conditions; Gu et al. [28] proposed a fast non-separable Gaussian process, combined it with linear regression to predict methylation levels, and demonstrated its superiority; Zhang et al. [29] introduced a multi-output Gaussian process simulator with a non-separable auto-covariance function; Cui et al. [30] applied the fast non-separable Gaussian process model to robust parameter design. In summary, researchers worldwide have investigated non-separable Gaussian processes from the perspectives of kernel construction and algorithmic optimization, achieving several practical applications. Nevertheless, this line of research has not yet received widespread attention.

This paper extends the quality gain-loss function by incorporating the correlation between spatiotemporal factors and the weight distribution of multi-dimensional quality characteristics on the basis of overall quality gain-loss theory. In response to the limitations of existing related research—which often neglects factor interactions and struggles to balance robustness with computational efficiency when handling multidimensional quality characteristics—this study introduces the NSGP, thereby overcoming the assumption of separability between spatial and temporal factors. Within this framework, a multidimensional quality gain-loss function model reflecting spatiotemporal interaction effects is constructed. The proposed approach not only improves predictive accuracy and computational efficiency but also enables the optimized design of multiple quality characteristics in complex engineering projects. It provides new theoretical support and practical application pathways for quality control in hydraulic concrete construction and related engineering fields.

2. Nonseparable Gaussian Proess Regression Model

2.1. Quality Characteristic Data

Let $\left\{y_i(s_j)\mid i=1,2,\cdots ,K,j=1,2,\cdots ,N\right\}$ denote the quality characteristic dataest of $K=k+k^{\ast }$ samples across $N=n+n^{\ast }$ observation locations/time points. Here $y_i(s_j)$ represents the quality characteristic value of the ith sample at the jth location/time point.

Define two sets of observation locations or time points , denoted as $s^D=\{s_1^D,s_2^D,\cdots ,s_{\mathrm{n}}^D\}$ and $s^{\ast }=\{s_1^{\ast },s_2^{\ast },\cdots ,s_n^{\ast }\}$. For $S^D$ , the quality characteristic level $s^D=\{s_1^D,s_2^D,\cdots ,s_{\mathrm{n}}^D\}$ can be observed for all $K=k+k^{\ast }$ samples .In contrast , for $S^{\ast }$, the quality characteristic level $s^{\ast }=\{s_1^{\ast },s_2^{\ast },\cdots ,s_n^{\ast }\}$ is only observed for the first k samples , while it is unobservable for the remaining $k^{\ast }$ samples.Therefore, for $S^D$, the quality characteristic values $y{(s^D)}_{k\times n}$ and $y^{\ast }{(s^D)}_{k^{\ast }\times n}$ are available for all samples. For $S^{\ast }$, only the quality characteristic values $y{(s^{\ast })}_{k\times n^{\ast }}$ for the first k samples can be obtained, while the quality characteristic values $y^{\ast }{(s^{\ast })}_{k^{\ast }\times n^{\ast }}$ for the remaining $k^{\ast }$ samples are unobservable. The objective of this paper is to predict the distribution of $y^{\ast }(s^{\ast })$ , given the known quality characteristic values $y(s^D)$, $y(s^{\ast })$, and $y^{\ast }(s^D)$.

2.2. Nonseparable Gaussian Process Modeling

Starting from the common output region $Y(s)$ ,the $K\times 1$-dimensional output vector at each location or time point s ∈ φ can be expressed as:

$$Y(s)=A\tilde{v}(s)+{\epsilon }_0$$

(1)

Where ${\epsilon }_0~\mathrm{MN}(0,{\sigma }_0^2{I}_K)$ , and $A=(a_1,\cdots ,a_K)$ are the $K\times K$ matrix, and $a_i$ is the i-th factor loading vector. The vector $\tilde{v}(s)={\left({\tilde{v}}_1(s),\cdots ,{\tilde{v}}_K(s)\right)}^{{\rm T}}$ denotes the random factor vector at observation location or time point s. If each random factor process ${\tilde{v}}_i(·)$ is independently modeled as a zero-mean GP model , with its noise term belonging to the class $s\in \phi$ , then for any observation location or time point s in $\phi$ , we have:

$$\begin{array}{l}{\tilde{v}}_i(s)=v_i(s)+{\epsilon }_i\\ v_i(·)~GP(0,{\tau }_i^2{c}_i(·,·))\end{array}$$

(2)

Among these , ${ϵ}_i\sim N(0,{\sigma }_i^2)$ represents the independent noise , and ${\tau }_i^2{c}_i(s_a,s_b)$ denotes the covariance between any pair of location or time point $s_a,s_b\in \phi$ .The parameter ${\sigma }_i^2$, along with ${\tau }_i^2$ , correspond to their respective unknown varaances.

Let the signal-to-noise ratio parameter be denoted as ${\eta }_i={\displaystyle \frac{{\sigma }_i^2}{{\tau }_i^2}}$ . From equation (2), the marginal distribution of ${\tilde{v}}_i(s^D)={\left({\tilde{v}}_i(s_1^D),\cdots ,\widetilde{v_i}(s_n^D)\right)}^{{\rm T}}$ follows a multivariate normal distribution, expressed as:

$${\tilde{v}}_i(s^D)|{\tau }_i^2,{\eta }_i~MN(0,{\tau }_i^2(R_i+{\eta }_i{I}_n))$$

(3)

Where the correlation matrix Ri has the following entries for $(\mathrm{l},\mathrm{m})$ : $c_i(s_1^D,s_m^D)$, $i=1,2,\cdots ,K$ , $I_n$ is the identity matrix for $n\times n$ . Using the Matérn kernel to establish a correlation model yields:

$$c_i(d)={\displaystyle \frac{1}{2^{v_i-1}\Gamma (v_i)}}{\left({\displaystyle \frac{d}{{\gamma }_i}}\right)}^{v_i}{\kappa }_{v_i}\left({\displaystyle \frac{d}{{\gamma }_i}}\right)$$

(4)

Where $\Gamma (\nu )$ is the gamma function , ${\kappa }_{v_i}(·)$ is a modified Bessel function of the second kind with smoothing parameter , ${\gamma }_i$ is the range parameter . For any $s_a,s_b\in \phi$ , there is $d=\left|\left.s_a-s_b\right|\right.$.

Let $Y(s^D)={\left(y{(s^D)}^T；y^{\ast }{(s^D)}^T\right)}^T$. Following Higdon's method, apply singular value decomposition(SVD) to $Y\left(s^D\right)=UDV$. Furthermore, the computational order for estimating A is linearly related to n , and when $i\ne j$ , there is $a_i^T{a}_j=0$. Thus , $A^{T}A$becomes a diagonal matrix, and therefore this paper estimates A as $A={\displaystyle \frac{UD}{\sqrt{n}}}$.

For convenience , vectorizing the outputs of $Y_v(s^D)$ and ${\tilde{v}}_v(s^D)$ yields $Y_v(s^D):=\mathrm{vec}(Y(s^D))$ and ${\tilde{v}}_v(s^D):=\mathrm{vec}(\tilde{v}{(s^D)}^T)$ . Furthermore , defining a matrix $Kn\times Kn$ as $A_v:=\left[I_n\otimes a_1；\cdots ；I_n\otimes a_K\right]$ , $Y\left(s^D\right)$ can be expressed as:

$$Y_v(s^D)=A_v{\tilde{v}}_v(s^D)+{\epsilon }_{0v}$$

(5)

Among these ,${\epsilon }_{0v}~MN(0,{\sigma }_0^2{I}_{nK})$ , $({\tilde{v}}_v(s^D)|{\tau }_1^2,\cdots ,{\tau }_K^2,{\tilde{R}}_1,\cdots ,{\tilde{R}}_K)~MN(0,{\Sigma }_v)$ , and the $\left(l,m\right)$ item is defined as ${\tilde{c}}_i(s_i^D,s_m^D)=c_i(s_l^D,s_m^D)+{\eta }_i{1}_{l=1}$. In addition , ${\displaystyle {\sum }_v=}\mathrm{blkdiag}({\tau }_1^2{\tilde{R}}_1,\cdots ,{\tau }_k^2{\tilde{R}}_k)$ , and $\mathrm{blkdiag}(·)$ denotes the block-diagonal matrix associated with the observation locations or time points . Marginalizing directly over $\widetilde{v_v}(s^D)$ yields :

$$Y_v(s^D)|{\sigma }_0^2,{\tau }_1^2,\cdots ,{\tau }_K^2,{\tilde{R}}_1,\cdots ,{\tilde{R}}_K~MN(0,A_v{\Sigma }_v{A}_v^T+{\sigma }_0^2{I}_{Kn})$$

(6)

However , directly computing the likelihood function in (6) requires repeated manipulations of the $Kn\times Kn$ covariance matrix $A_v{\Sigma }_v{A}_v^T+{\sigma }_0^2{I}_{Kn}$ , resulting in computational complexity of $0\left(\left.K{n}^3\right)\right.$.

Lemma 1: Let $A=UD/\sqrt{n}$ , $A{(A^{{\rm T}}A)}^{-1}{A}^{{\rm T}}=I_K$. Then the marginal likelihood function of $Y_v(s^D)$ is :

(7)

Among these , $\widehat{v}(s^D)={(A_v^T{A}_v)}^{-1}{A}_v^T{Y}_v(s^D)$. Marginalizing over $v_v(s^D)$ yields :

(8)

Thus, the marginal likelihood in model (1) can be expressed as the product of K independent multivariate normal distributions:

$$p(Y(s^D)|{\sigma }_0^2,{\tau }_1^2,\cdots ,{\tau }_K^2,{\tilde{R}}_1,\cdots ,{\tilde{R}}_K)=|A_v^T{A}_v{|}^{-1/2}{\displaystyle \prod _{i=1}^K{p}_{MN}}({\widehat{v}}_i(s^D);0,{\tau }_i^2{\tilde{R}}_i+{\sigma }_0^2{(a_i^T{a}_i)}^{-1}-I_n)$$

(9)

Where $p_{MN}\text{}(·;\mu ,\Sigma )$ follows a multivariate normal distribution with mean μ and covariance matrix Σ , and ${\widehat{v}}_i(s^D)$ denotes the transpose of the i-th row of the data table:

$$\widehat{v}(s^D)={(A^{T}A)}^{-1}{A}^{T}Y(s^D)$$

(10)

Lemma 1 indicates that the marginal likelihood function of model (1) can be expressed as the product of K independent multivariate normal distributions, thereby greatly simplifying the computational burden. Specifically, instead of evaluating the density of a multivariate normal distribution with a covariance matrix of dimension Kn×Kn, the density can be estimated by computing K independent multivariate normal distributions , each with a covariance matrices of dimension n × n , rather than using an independent multivariate normal distribution with a covariance matrix of dimension Kn×Kn . However, when n is large, the computational cost of directly estimating the density via independent multivariate normal distributions with n × n covariance matrices remains prohibitively high. Suppose:

$$Y(s_j^{\ast })={(y{(s_j^{\ast })}^T;y^{\ast }{(s_j^{\ast })}^T)}^T$$

(11)

Among these , $y(s_j^{\ast })$ and $y^{\ast }(s_j^{\ast })$ denote the j-th row of $y(s^{\ast })$ and $y^{\ast }(s^{\ast })$ respectively.

Lemma 2: Let $A=UD/\sqrt{n}$, $A{(A^{{\rm T}}A)}^{-1}{A}^{{\rm T}}=I_K$, $Y_v(s^D,s_j^{\ast }):=\mathrm{vec}(Y(s^D,s_j^{\ast }))$ and $v_v(s^D):=\mathrm{vec}(v{(s^D)}^T)$ , where all of these are $k\times (n+1)$-dimensional matrices. Vectorizing the outputs of $Y_v(s^D,s_j^{\ast })$ and $v_v(s^D)$ yields $Y_v(s^D,s_j^{\ast }):=\mathrm{vec}(Y(s^D,s_j^{\ast }))$ and $v_v(s^D):=\mathrm{vec}(v{(s^D)}^T)$, respectively . Accordingly , equation (1) can be written as:

$$Y_v(s^D,s_j^{\ast })=A_v{v}_v(s^D,s_j^{\ast })+\epsilon$$

(12)

Among these , $\epsilon ~MN(0,{\sigma }_0^2{I}_{(n+1)K})$.

1) For each $s_j^{\ast }$ , with $j=1,2,\cdots ,n^{\ast }$ . Setting $r_i(s_j^{\ast })={(c_i(s_j^{\ast },s_1^P),\cdots ,c_i(s_j^{\ast },s_n^D))}^T$ yields:

$$Y(s_j^{\ast })|Y(s^D),{\sigma }_0^2,{\tau }_{1:K}^2,{\eta }_{1:K},\kappa ~MN(\widehat{\mu }(s_j^{\ast }),\widehat{\Sigma }(s_j^{\ast }))$$

(13)

Among these , $\widehat{\mu }(s_j^{\ast })=A{\widehat{v}}^{\ast }(s_j^{\ast })$, and the i-th term of ${\widehat{v}}^{\ast }(s_j^{\ast })$is ${\widehat{v}}_i^{\ast }(s_j^{\ast })=r_i^T(s_j^{\ast })$. Additionally , for ${({\tilde{R}}_i+{\displaystyle \frac{{\sigma }_0^2{\left(a_i^T{a}_i\right)}^{-1}}{{\tau }_i}}{I}_n)}^{-1}{\widehat{v}}_i(s^D)$ , ${\widehat{v}}_i(s^D)$ is the transpose of the i-th term defined in formula(6). Furthermore , $\widehat{\Sigma }(s_j^{\ast })=A{D}^{\ast }(s_j^{\ast })A^T+{\sigma }_0^2{I}_K$ , where $D^{\ast }(s_j^{\ast })$ is a diagonal matrix whose i-th diagonal entry is:

$${\tau }_i^2\left.\left[{\widehat{c}}_i(s_j^{\ast },s_j^{\ast })-r_i^T(s_j^{\ast }){\left({\widehat{R}}_i+{\displaystyle \frac{{\sigma }_0^2(a_i^T{(a_i^T{a}_i)}^{-1}{\alpha }_i)}{{\tau }_i^2}}{I}_n\right)}^{-1}{r}_i(s_j^{\ast })\right.\right],i=1,2,\cdots ,K$$

(14)

2) Let the partition be defined as $\widehat{\mu }(s_j^{\ast })={({\widehat{\mu }}_0^T(s_j^{\ast }),{\widehat{\mu }}_{\ast }^T(s_j^{\ast }))}^T$ and $\widehat{\Sigma }(s_j^{\ast })=\left(\begin{array}{cc}{\widehat{\Sigma }}_{00}(s_j^{\ast })& {\widehat{\Sigma }}_{0\ast }(s_j^{\ast })\\ {\widehat{\Sigma }}_{\ast 0}(s_j^{\ast })& {\widehat{\Sigma }}_{\ast \ast }(s_j^{\ast })\end{array}\right)$ . Then, for any $s_j^{\ast }$ , the prediction distribution of $y^{\ast }(s_j^{\ast })$ is :

$$y^{\ast }(s_j^{\ast })|y(s^D),y(s^{\ast }),y^{\ast }(s^D),{\sigma }_0^2,{\tau }_{1:K}^2,{\gamma }_{1:K},{\eta }_{1:K}~MN({\widehat{\mu }}_{\ast |0}(s_j^{\ast }),{\widehat{\Sigma }}_{\ast |0}(s_j^{\ast }))$$

(15)

Among these ,${\widehat{\mu }}_{\ast }(s_j^{\ast })$ corresponds to${\widehat{\mu }}_{\ast |0}(s_j^{\ast })={\widehat{\mu }}_{\ast }(s_j^{\ast })+{\widehat{\Sigma }}_{\ast 0}(s_j^{\ast }){\widehat{\Sigma }}_{00}^{-1}(y(s_j^{\ast })-{\widehat{\mu }}_0(s_j^{\ast }))$ while ${\widehat{\Sigma }}_{\ast |0}(s_j^{\ast })$ corresponds to ${\widehat{\Sigma }}_{\ast |0}(s_j^{\ast })={\widehat{\Sigma }}_{\ast \ast }(s_j^{\ast })-{\widehat{\Sigma }}_{\ast 0}(s_j^{\ast }){\widehat{\Sigma }}_{00}^{-1}(s_j^{\ast }){\widehat{\Sigma }}_{0\ast }(s_j^{\ast })$.

2.3. Model Estimation and Computational Strategy

Computing the likelihood function in model (1) requires evaluating the inverse of Ri , where each inversion incurs a computational complexity cost of 0(n3) . Consequently , when n=K is large, the overall computational complexity escalates to 0(Kn3) , rendering direct implementation computationally infeasible. Therefore, this paper introduces a computationally efficient algorithm grounded in Gaussian random fields and Gaussian Markov random fields (GMRF), inspired by the work of Hartikainen and Särkkä . Compared with existing approaches, this proposed approach offers the notable advantage of avoiding any approximation of the likelihood function.

Consider a continuous regression model of order p defined by the following stochastic differential equation (SDE):

$$c_p{v}^{(p)}(s)+c_{p-1}{v}^{(p-1)}(s)+\cdots +c_{0}v(s)=b_{0}z(s)$$

(16)

Where $v^{(l)}(s)$ denotes the l-th derivative coefficient of the process $v(s)$, and$z(s)$ represents a standard Gaussian white noise process defined on the domain$s\in R$ . To avoid identifiability issues , we set $c_p=1$ without loss of generality. According to the work of Whittle , the spectral density of a Gaussian process equipped with a Matérn kernel is given by:

$$S_{\mathrm{Mat}}(t)\propto {\displaystyle \frac{1}{{({\lambda }^2+t^2)}^{(v+1/2)}}}$$

(17)

Among these, $\lambda ={\displaystyle \frac{\sqrt{2\nu }}{\gamma }}$ , $\gamma$is the coverage parameters , and $\nu$ is the smoothing parameter.

2.4. Computation of Continuous-Time Stochastic Processes

By Lemma1 , the likelihood of $Y(s^D)$ in model(1) is proportional to the product of K-multivariate normal distributions ${\widehat{v}}_i(s^D)$. Here , ${\widehat{v}}_i(s^D)$ denotes the transpose of the i-th row of $\widehat{v}(s^D)={(A^{T}A)}^{-1}{A}^{T}Y(s^D)$. Furthermore, due to the identifiable row property, we focus on $\widehat{v}(s^D)$ when $v_0^2=0$.

Choose The Matérn kernel function when $v=\raisebox{1ex}{$2$}\!\left/ \!\raisebox{-1ex}{$5$}\right.$ as the latent factor process:

$$c(d)=\left(1+{\displaystyle \frac{\sqrt{5}d}{\gamma }}+{\displaystyle \frac{5d^2}{3{\gamma }^2}}\right)\exp \left(-{\displaystyle \frac{\sqrt{5}d}{\gamma }}\right)$$

(18)

Define ${\theta }_i(s):={({\nu }_i(s),v_i^{(1)}(s),v_i^{(2)}(s))}^T$, and for $\forall s\in \phi$, each $v_i^{(l)}(s)$ corresponds to the l-th order derivative (with l=1,2) of the underlying latent process with respect to s .

Therefore , for any $s\in \phi$ , equation (2) can be expressed as:

$$\begin{array}{l}{\widehat{v}}_i(s)=v_i(s)+{\epsilon }_i\\ v_i(·)~\mathrm{GP}(0,{\tau }_i^2{c}_i(·,·))\end{array}$$

(19)

Combining the Gaussian process model in equation (21) with the correlation structure implied by the stochastic differential equation (SDE) in equation (18) yields the following unified representation:

$${\displaystyle \frac{d{\theta }_i(s)}{ds}}=J_i{\theta }_i(s)+L{z}_i(s)$$

(20)

Among these , $z_i(s)$ is a Gaussian white noise process with mean = 0 and variance =${\tau }_i^2$. And , $J_i=\left[\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\\ -{\lambda }_i^3& -{\lambda }_i^2& -3{\lambda }_i\end{array}\right]$, ${\lambda }_i=\sqrt{2v_i}/{\gamma }_i,i=1,2,\cdots ,K$. $L={(0,0,1)}^{{\rm T}}$.

For $j=1,\cdots ,N-1$, let $d_j=|s_{j+1}-s_j|$. Then the expression for the solution of the SDE in equation (20) is:

$$G_i(s_j)=e^{{\lambda }_i{d}_j}={\displaystyle \frac{e^{-{\lambda }_i{d}_j}}{2}}\left[\begin{array}{ccc}{\lambda }_i^2{d}_j^2+2{\lambda }_i+2& 2({\lambda }_i^2{d}_j^2+d_j)& d_j^2\\ -{\lambda }_i^3{d}_j^2& -2({\lambda }_i^3{d}_j^2-{\lambda }_i{d}_j-1)& 2d_j-{\lambda }_i{d}_j^2\\ {\lambda }_i^4-2{\lambda }_i^3{d}_j& 2({\lambda }_i^3{d}_j^2-3{\lambda }_i^2{d}_j)& {\lambda }_i^2{d}_j^2-4{\lambda }_i{d}_j+2\end{array}\right]$$

$$W_i(s_j)={\displaystyle \frac{4{\sigma }^2{\lambda }_i^5}{3}}\left[\begin{array}{ccc}{W}_{1i}(s_j)& W_{2i}(s_j)& W_{3i}(s_j)\\ W_{4i}(s_j)& W_{5i}(s_j)& W_{6i}(s_j)\\ W_{7i}(s_j)& W_{8i}(s_j)& W_{9i}(s_j)\end{array}\right]$$

Among these:

$$\begin{array}{l}{W}_{1i}(s_j)={\displaystyle \frac{e^{-2{\lambda }_i{d}_j}(3+6{\lambda }_i{d}_j+6{\lambda }_i^2{d}_j^2+4{\lambda }_i^3{d}_j^3+2{\lambda }_i^4{d}_j^4)-3}{-4{\lambda }_i^5}}\\ W_{2i}(s_j)=W_{4i}(s_j)={\displaystyle \frac{e^{-2{\lambda }_i{d}_j}{t}^4}{2}}\\ W_{3i}(s_j)=W_{7i}(s_j)={\displaystyle \frac{e^{-2{\lambda }_i{d}_j}(1+2{\lambda }_i{d}_j+2{\lambda }_i^2{d}_j^2+4{\lambda }_i^3{d}_j^3-2{\lambda }_i^4{d}_j^4)-1}{4{\lambda }_i^3}}\\ W_{5i}(s_j)={\displaystyle \frac{e^{-2{\lambda }_i{d}_j}(1+2{\lambda }_i{d}_j+2{\lambda }_i^2{d}_j^2-4{\lambda }_i^3{d}_j^3+2{\lambda }_i^4{d}_j^4)-1}{-4{\lambda }_i^3}}\\ W_{6i}(s_j)=W_{8i}(s_j)={\displaystyle \frac{e^{-2{\lambda }_i{d}_j}{d}_j^2(4-4{\lambda }_i{d}_j+{\lambda }_i^2{d}_j^2)}{2}}\\ W_i(s_0)=\left[\begin{array}{ccc}{\sigma }_i^2& 1& -{\sigma }_i^2{\lambda }_i^2/3\\ 0& {\sigma }_i^2{\lambda }_i^2/3& 1\\ -{\sigma }_i^2{\lambda }_i^2/3& 0& {\sigma }_i^2{\lambda }_i^4\end{array}\right]\end{array}$$

For $j=1,\cdots ,N$ , the joint distribution of ${\theta }_i(s_{0:n})$ is:

$$\left[\begin{array}{c}{\theta }_i(s_0)\\ {\theta }_i(s_1)\\ {\theta }_i(s_2)\\ ⋮\\ {\theta }_i(s_N)\end{array}\right]~MN\left(\left[\begin{array}{c}0\\ 0\\ 0\\ ⋮\\ 0\end{array}\right],{\left[\begin{array}{cccc}{W}_i{(s_0)}^{-1}& -G_i^T(s_1)W_i{(s_1)}^{-1}{G}_i(s_1)& \cdots & \cdots \\ -G_i^T(s_1)W_i{(s_1)}^{-1}{G}_i(s_1)& W_i^{-1}(s_1)& \cdots & \cdots \\ ⋮& -G_i^T(s_2)W_i{(s_2)}^{-1}{G}_i(s_2)& \cdots & \cdots \\ ⋮& \cdots & \cdots & W_i^{-1}(s_N)\end{array}\right]}^{-1}\right)$$

Let $q_i={\displaystyle \frac{16}{3}}{\sigma }_i^2{\lambda }_i^5$ , $F=(1,0,0)$. Then the SDE solution for model (2) can be expressed as:

$$\begin{array}{l}{{\tilde{v}}_i}_{(s_{j+1})}=F{\theta }_i(s+1)+{\epsilon }_i\\ {\theta }_i(s_j+1)=G_i(s_j){\theta }_i(s_j)+w_i(s_j)\end{array}$$

(21)

Among these , $G_i(s_j)=e^{J_i(s_{j+1}-s_j)}$. $w_i\text{}(s_j\text{})\sim MN(0,W_i\text{}(s_j\text{}))$( i =1, … , K , j = 1, … , N-1) and $W_i(s_j)={\displaystyle {\int }_0^{s_{j+1}-s_j}{e}^{J_{i}t}}L{q}_i{L}^{{\rm T}}{e}^{J_i^{T}t}dt$. Thus, the stationary distribution at ${\theta }_i$ is:

$${\theta }_i(s_0)~MN(0,W_i(s_0))$$

(22)

Among these , $W_i(s_0)={\displaystyle {\int }_0^{\infty }{e}^{J_{i}t}}L{q}_i{L}^T{e}^{J_i^{T}t}dt$.

Lemma 3: (Kalman Filter) Consider the continuous-time state-space model specified by equation (21), assume all distributions in this lemma are conditioned on $\left({\sigma }_i^2,{\tau }_i^2,{\gamma }_i\right)$

. Furthermore, for any j ≥ 1 and i = 1, 2, … , K , there is :

$${\theta }_i(s_{j-1}^D)|{\widehat{v}}_i(s_{1:j-1}^D)~MN\left(\left.m_i(s_{j-1}^D),C_i(s_{j-1}^D)\right)\right.$$

(23)

(1) Given ${\widehat{v}}_i(s_{1:j-1}^D)$, the one-step-ahead prediction distribution of ${\theta }_i(s_j^D)$ is:

$${\theta }_i(s_j^D)|{\widehat{v}}_i(s_{1:j-1}^D)~MN\left(\left.{\mu }_{{\theta }_i}(s_j^D),{\theta }_i(s_j^D)\right)\right.$$

(24)

Among these, ${\mu }_{{\theta }_i}(s_j^D)=E{\theta }_i\left(\left.(s_j^D)|{\widehat{v}}_i(s_{1:j-1}^D)\right)\right.=G_i(s_j^D)m_i(s_{j-1}^D)$ , ${\theta }_i(s_j^D)=\mathrm{var}\left(({\theta }_i(s_j^D)|\right.$ , $\left.{\widehat{v}}_i(s_{1:j-1}^D)\right)=G_i(s_j^D)C_i(s_{j-1}^D)G_i^T(s_j^D)+W_i(s_j^D)$.

(2) Given ${\widehat{v}}_i(s_j^D)$, the one-step-ahead prediction distribution of ${\widehat{v}}_i(s_j^D)$ is:

$${\widehat{v}}_i(s_j^D)|{\widehat{v}}_i(s_{1:j-1}^D)~N\left(\left.f_i(s_j^D),Q_i(s_j^D)\right)\right.$$

(25)

(3) Given ${\widehat{v}}_i(s_{1:j}^D)$, the filtering distribution of ${\theta }_i(s_j^D)$ is:

$${\theta }_i(s_j^D)|{\widehat{v}}_i(s_{1:j}^D)~MN(m_i(s_j^D),C_i(s_j^D))$$

(26)

Among these, $m_i(s_j^D)=E[{\theta }_i(s_j^D)|{\widehat{v}}_i(s_{1:j}^D)]={\mu }_{{\theta }_i}(s_j^D)+{\Sigma }_{{\theta }_i}(s_j^D)F^T{Q}_i^{-1}(s_j^D)e_i(s_j^D)$, $e_i(s_j^D)$$=var[{\theta }_i\text{}(s_j^D\text{})|{\widehat{v}}_i(s_{1:j}^D)={\Sigma }_{{\theta }_i}(s_j^D)-{\Sigma }_{{\theta }_i}(s_j^D)F^T{Q}_i^{-1}(s_j^D)F{\Sigma }_{{\theta }_i}(s_j^D)$. The prediction error is $e_i(s_j^D)$$={\widehat{v}}_i(s_j^D)-f_i(s_j^D)$.

To compute the likelihood function, the Kalman filter on $s^D$ is defined in Lemma3, the estimate of $\left({\sigma }_i^2,{\tau }_i^2,{\gamma }_i\right)$ when i=1,2,…K is given, then the Kalman filtering perform on all observation locations or time points $s_{1:N}$ (including $s^{\ast }$).

Lemma 4: (Kalman smoother) Consider the continuous-time state-space model specified by equation (21). For simplicity, all distributions in this lemma are assumed to be conditioned on $\left({\sigma }_i^2,{\tau }_i^2,{\gamma }_i\right)$

. Furthermore,for any i = 1, 2, … , K and j = 1, 2, … , N-1, let ${\theta }_i(s_{j+1})|{\widehat{v}}_i(s^D)~N({\mu }_{{\theta }_i}^{\ast }(s_{j+1}),{\Sigma }_{{\theta }_i}^{\ast }(s_{j+1}))$, then:

$${\theta }_i(s_j)|{\widehat{v}}_i(s_j^D)~N({\mu }_{{\theta }_i}^{\ast }(s_j),{\Sigma }_{{\theta }_i}^{\ast }(s_j))$$

(27)

Among these, ${\mu }_{{\theta }_i}^{\ast }(s_j)=m_i(s_j)+C_i(s_j)G_i^T(s_{j+1}){\Sigma }_{{\theta }_i}^{-1}(s_{j+1})({\mu }_{{\theta }_i}^{\ast }(s_{j+1})-{\mu }_{{\theta }_i}(s_j))$, ${\Sigma }_{{\theta }_i}^{\ast }(s_j)=C_i(s_j)-C_i(s_j)G_i^T(s_{j+1}){\Sigma }_{{\theta }_i}^{-1}(s_{j+1})({\Sigma }_{{\theta }_i}(s_{j+1})-{\Sigma }_{{\theta }_i}^{\ast }(s_{j+1})){\Sigma }_{{\theta }_i}^{-1}(s_{j+1})G_i(s_{j+1})C_i(s_j)$.

Based on the above settings, the posterior distribution of ${\theta }_i(s)(i=1,2,\cdots ,K)$can be efficiently computed via the forward filtering and backward smoothing (FFBS) algorithm, requiring only $0(n)$operations—substantially reducing the computational burden compared with $0(n^3)$ operations cost of directly evaluating of the Gaussian process model. Predictions at $s\ast$ are computed linearly with the number of observed locations or time points, resulting in an overall computational complexity of only $0(n)$.

It is worth noting that when $v_i=2/5$, both the Matérn covariance matrix and its inverse in equation(18) are rank-n matrices of size n × n. However, the covariance matrix ${\theta }_i(s)={(v_i(s),v_i^{(1)}(s),v_i^{(2)}(s))}^T$ of the latent states is sparse. Therefore, the result holds for all Matérn kernels when $v_i=(2m+1)/2(m\in N)$.

2.5. Parameter Estimation

The computationally intensive part of the FFBS algorithm involves sampling the latent states and marginalizing them:

$$p(\widehat{v}(s^D)|{\sigma }_{1：K}^2,{\tau }_{1：K}^2,{\gamma }_{1：K})=\left.{\displaystyle \sum _{i=1}^K\left\{p\left.\left({\widehat{v}}_i(s_1^D)|{\sigma }_i^2,{\tau }_i^2,{\gamma }_i\right.\right){\displaystyle \prod _{j=2}^{n}p}\left({\widehat{v}}_i(s_j^D\right.)\right.}|\left.{\widehat{v}}_i(s_{1:j-1}^D),{\sigma }_i^2,{\tau }_i^2,{\gamma }_i\right)\right\}$$

(28)

Among these, each term follows the normal distribution specified in the one-step-ahead prediction of the FFSE algorithm.

Lemma 5: (Likelihood by the Kalman filter) For each i = 1, 2, … , K, there is:

$$p({\widehat{v}}_i(s_{1:j}^D)|{\sigma }_i^2,{\tau }_i^2,{\gamma }_i)={(2\pi )}^{-n/2}{\displaystyle \prod _{j=1}^N{Q}_i^{-1/2}}(s_j^D)\exp \left\{-{\displaystyle \sum _{j=1}^N\frac{{({\widehat{v}}_i(s_j^D)-f_i(s_j^D))}^2}{Q_i(s_j^D)}}\right\}$$

(29)

Among these, $f_i(s_j^D)$ and $Q_i(s_j^D)$ are defined as in lemma 3. To complete the model, assume the following prior:

$$\pi ({\sigma }_{1:K}^2,{\tau }_{1:K}^2,{\gamma }_{1:K})\propto {\displaystyle \prod _{i=1}^K\left(\frac{\pi ({\tau }_i^2,{\gamma }_i)}{{\sigma }_i^2}\right)}$$

(30)

Let ${\zeta }_i:=1/{\gamma }_i$ and ${\eta }_i:={\tau }_i/{\sigma }_i^2$ . For i = 1, 2, … , K, assume the joint robust prior for the transformation parameters is : $\pi ({\zeta }_i,{\eta }_i)\propto {(C_i{\zeta }_i+{\eta }_i)}^{a_i}\exp (-b_i(C_i{\zeta }_i+{\eta }_i))$. Here, $a_i$ , $b_i$ and $C_i$ are prior parameters. If set $a_i=1/2$, $b_i=1$ and $C_i=\left|\phi \right|n$,where $\left|\phi \right|$ is the length of $\phi$. When n is large, sampling the posterior distribution via MCMC remains slow. Parameters can be estimated using the posterior mode:

$$({\widehat{\zeta }}_i,{\widehat{\eta }}_i)=\arg {\max }_{{\zeta }_i,{\eta }_i}\left\{p(\widehat{v}(s^D)|{\zeta }_i,{\eta }_i)\pi ({\zeta }_i,{\eta }_i)\right\}$$

(31)

By applying Lemma 5, the computation can be simplified to $0(n)$ operations.The NSGP model construction flowchart is shown in Figure 1.

3. Multivariate and Multidimensional Quality Profit and Loss Function Based on Non-separable Gaussian Process

3.1. Signal-to-Noise Ratio

The quality characteristic value of a product is a random variable influenced by multiple factors. Assuming the mathematical expectation is μ and the variance is σ2. For nominal the best, the closer μ is to the target value T, the better ; the smaller σ2 is, the better. In robust design, let the coefficient of variation be y and the signal-to-noise ratio be η , then :

$$\gamma =\frac{\sigma }{\mu }$$

(32)

$$\eta =\frac{{\mu }^2}{{\sigma }^2}$$

(33)

The signal-to-noise ratio ${{\displaystyle \eta }}_T$ for nominal the best is:

$${{\displaystyle \eta }}_T=10\log {\displaystyle \frac{{\mu }^2}{{\sigma }^2}}$$

(34)

Setting $\overline{y}={\displaystyle \frac{1}{n}}{\displaystyle \sum _{i=1}^n{y}_i}$,$S^2={\displaystyle \frac{1}{n-1}}{\displaystyle \sum _{i=1}^n(y_i}-\overline{y}{)}^2$, the estimated value $\widehat{{{\displaystyle \eta }}_T}$ for ${{\displaystyle \eta }}_T$ is:

$$\widehat{{{\displaystyle \eta }}_T}=10\log \left({\displaystyle \frac{{\overline{y}}^2}{S^2}}\right)$$

(35)

Similarly, larger the better signal-to-noise ratio ${{\displaystyle \eta }}_L$ and its estimated value $\widehat{{{\displaystyle \eta }}_L}$ are:

$${{\displaystyle \eta }}_L=10\log \left({\sigma }^2+{\mu }^2\right)$$

(36)

$$\widehat{{{\displaystyle \eta }}_L}=-10\log \left({\displaystyle \frac{1}{n}}{\displaystyle \sum _{i=1}^n\frac{1}{{y_i}^2}}\right)$$

(37)

Smaller the better signal-to-noise ratio ${{\displaystyle \eta }}_S$ and its estimate $\widehat{{{\displaystyle \eta }}_S}$ are:

$${{\displaystyle \eta }}_S=10\log {\displaystyle \frac{1}{{\mu }^2+{\sigma }^2}}$$

(38)

$$\widehat{{{\displaystyle \eta }}_S}=-10\log \left({\displaystyle \frac{1}{n}}{\displaystyle \sum _{i=1}^n{{{\displaystyle y}}_i}^2}\right)$$

(39)

The signal-to-noise ratio is a monotonic function of the coefficient of variation $\gamma$ ,serving as a crucial indicator of product quality stability. A higher $\eta$ value indicates greater product stability and reduced losses.

Both the coefficient of variation $\gamma$ and the signal-to-noise ratio $\eta$ are dimensionless data. They both simultaneously consider the effects of mean and variance on quality stability. However, $\eta$ is closer to a normal distribution, enabling linear additivity between quality characteristics and allowing interaction effects to be neglected.

3.2. Quality Gain-Loss Function Model

Leveraging the advantages of Taguchi's signal-to-noise ratio method and Pignatiello's multivariate quadratic loss function approach in robust parameter design, a multivariate quality gain-loss function suitable for spatiotemporal responses is constructed.

Based on the signal-to-noise ratio method, the loss weights ${\alpha }_i$ and ${\beta }_i$ for each spatial factor and temporal factor level under the quality characteristic $y_i$ are calculated as:

$${\alpha }_i={\displaystyle \frac{\frac{1}{{\eta }_i}}{{\displaystyle \sum _{l=1}^K\frac{1}{{\eta }_l}}}},i=1,2,\cdots K$$

(40)

$${\beta }_i={\displaystyle \frac{\frac{1}{{\eta }_j}}{{\displaystyle \sum _{m=1}^N\frac{1}{{\eta }_m}}}},j=1,2,\cdots ,N$$

(41)

Among these, ${\alpha }_i$ and ${\beta }_i$ represent the signal-to-noise ratios for the quality characteristics located in row i and column j respectively. Since a higher signal-to-noise ratio indicates greater stability in the quality characteristic and consequently smaller quality losses, ${\alpha }_i$ and ${\beta }_i$ can respectively determine the contribution levels of the spatial factor i and the temporal factor j to the fluctuations.

Construct the joint weight matrix using gain-loss weights:

$$\omega ={\left({\alpha }_1,{\alpha }_2,\cdots ,{\alpha }_K\right)}^T\times \left({\beta }_1,{\beta }_2,\cdots ,{\beta }_N\right)={\left(\begin{array}{cccc}{\alpha }_1{\beta }_1& {\alpha }_1{\beta }_2& \cdots & {\alpha }_1{\beta }_N\\ {\alpha }_2{\beta }_1& {\alpha }_2{\beta }_2& \cdots & {\alpha }_2{\beta }_N\\ ⋮& ⋮& \ddots & ⋮\\ {\alpha }_K{\beta }_1& {\alpha }_K{\beta }_2& \cdots & {\alpha }_K{\beta }_N\end{array}\right)}_{K\times N}$$

(42)

Among these , the elenment (i,j) of $\omega$ ,denoted as ${\omega }_{ij}={\alpha }_i{\beta }_j$ , represents the joint loss weight under the combined effect of the ith spatial factor level and the jth temporal factor level.

The multi-dimensional quality gain-loss functions for larger the better ($G_{GL}(y(x))$) , smaller the better ($G_{GS}(y(x))$) , and nominal the best ($G_{GT}(y(x))$) quality characteristics are respectively:

(43)

(44)

(45)

Let:

(46)

Then Equations (43) to (45) can be unified into the following form:

$$G(y(x))=={\displaystyle \sum _{i=1}^K{\displaystyle \sum _{j=1}^N{\omega }_{ij}}}G(y_i(x_j))$$

(47)

Among these , $y_i(x_j)$ represents the observed value of quality characteristic i at time point j for sample $y_i$ . ${\omega }_{ij}$ and $G(y_i(x_j))$ denote the joint gain-loss weight and the gain-loss function, respectively, under the combined influence of spatial factor i and time factor j .

3.3. Two-Stage Parameter Optimization

3.3.1. Time Factor Optimization

Using the NSGP model and the quality gain-loss function, the spatiotemporal response model is constructed as:

$$\arg \min G(y_{1:K}(s_{1:N}))={\displaystyle \sum _{i=1}^K{\displaystyle \sum _{j=1}^N{\omega }_{ij}}}G(y_i(s_j))$$

(48)

A nonlinear optimization algorithm is employed to find the globally optimal parameter design value $s^{\mathrm{opt}}=(s_1^{\mathrm{opt}},s_2^{\mathrm{opt}},\cdots ,s_N^{\mathrm{opt}})$ for the time factor s in equation(48). The FFBS algorithm calculates the optimal solution $s^{\mathrm{opt}}$ corresponding to the overall predicted mean.

3.3.2. Spatial Factor Optimization

Use the overall prediction mean from the first stage as the observation value for this section. Let $z_j$ denote the mass characteristic corresponding to the jth position/time point, and ${\beta }_j$ denote the gain/loss weight for mass characteristic $z_j$ . Then there is:

$$G(z_j(x))=\left\{\begin{array}{ll}{k}_i{(z_j(x)-T_j)}^2+g(z_j(x)-T_i)\hfill & NTB\hfill \\ k_i{(z_j(x))}^{-2}+g(z_j{(x)}^{-1})\hfill & LTB\hfill \\ k_i{(z_j(x))}^2+g(z_j(x))\hfill & STB\hfill \end{array}\right.$$

(49)

Employing multi-response robust parameter modeling and its design methodology, construct the multi-response optimization model:

$$\arg \min G(z(x))={\displaystyle \sum _{j=1}^N{\beta }_j}G(z_j(x))$$

(50)

Among these , ${\beta }_j$ and $G(z_j(x))$ represent the quality gain/loss weight and gain/loss function for the jth quality characteristic $z_j$ , respectively. Employ a nonlinear optimization algorithm to find the globally optimal parameter design values $x^{\mathrm{opt}}=(x_1^{\mathrm{opt}},x_2^{\mathrm{opt}},\cdots ,x_N^{\mathrm{opt}})$ for the space factor x in equation (50).

By integrating the optimal solutions from equations (48) and (50), the combined optimal parameter combination for the temporal and spatial factors is obtained: $(x^{opt},s^{opt})$.

The flowchart integrating NSGP with the multi-objective, multi-dimensional quality gain-loss function model is shown in Figure 2.

4. Simulation and Case Studies

4.1. Prediction Indicators

This paper employs simulation cases and concrete strength prediction examples to demonstrate the superiority of the proposed method. The following models were selected for comparison: univariate Gaussian process (GP), separable Gaussian process(SGP), two linear regression (LR) strategies (sample-based and time-point-based), two random decision forests (RF) strategies (sample-based and time-point-based), and nearest neighbor (NN) method. The predictive performance of these models was evaluated using the following metrics:

(1) Root mean square error:

$$\mathrm{RMSE}=\sqrt{{\displaystyle \frac{{\displaystyle \sum _{i=1}^{k^{\ast }}{\displaystyle \sum _{j=1}^{n^{\ast }}(y_i^{\ast }(}}{s}_j^{\ast })-{\widehat{y}}_i^{\ast }(s_j^{\ast }){)}^2}{k^{\ast }{n}^{\ast }}}}$$

(51)

(2) 95% posterior confidence interval coverage rate:

$$P_{\mathrm{CI}}(95\%)={\displaystyle \frac{1}{k^{\ast }{n}^{\ast }}}{\displaystyle \sum _{i=1}^{k^{\ast }}{\displaystyle \sum _{j=1}^{n^{\ast }}1}}\{y_i^{\ast }(s_j^{\ast })\in {\mathrm{CL}}_{ij}(95\%)\}$$

(52)

(3) 95% posterior confidence interval length:

$$L_{\mathrm{CI}}(95\%)={\displaystyle \frac{1}{k^{\ast }{n}^{\ast }}}{\displaystyle \sum _{i=1}^{k^{\ast }}{\displaystyle \sum _{j=1}^{n^{\ast }}\mathrm{length}}}\{{\mathrm{CI}}_{ij}(95\%)\}$$

(53)

(4) Multivariate quality gain-loss function:

$$G(y^{\ast }(s^{\ast }))={\displaystyle \sum _{i=1}^{k^{\ast }}{\displaystyle \sum _{j=1}^{n^{\ast }}{\omega }_{ij}}}G(y_i^{\ast }(s_j^{\ast }))$$

(54)

For $1\le i\le k^{\ast }$ and $1\le j\le n^{\ast }$ , $y_i^{\ast }(s_j^{\ast })$ is the observation value of the $s_j^{\ast }$th sample at the jth position/timestamp $s_j^{\ast }$, ${\widehat{y}}_i^{\ast }(s_j^{\ast })$ is the observation value of the ith sample at the jth position/timestamp $s_j^{\ast }$ , and ${\omega }_{ij}$ is the joint gain-loss weight of the ith sample at the jth position/timestamp $s_j^{\ast }$ . ${\mathrm{CI}}_{ij}(95\%)$ represents the 95% posterior confidence interval for the predicted mean.

4.2. Simulation Case

4.2.1. Simulation Function

In many complex engineering and scientific problems, the optimal response of a system —such as the maximum strength of concrete or the highest yield of crops—is often not a constant value but is concentrated within a specific region that dynamically varies with spatiotemporal factors. To capture the correlation between dynamic peak characteristics and process parameters, shape parameters and time parameters, this paper constructs a simulation function based on the Space-time Gaussian Pulse function, with the specific form as follows:

$$f(x_1,x_2,t)=[8e^{{-0.05\left.\left(x_1-x_{1_0}\right.\right)}^2}]e^{\frac{-{(t-t_0)}^2}{2w_t^2}+{\displaystyle \frac{{\left.-\left(x_2-x_{2_0}\right.\right)}^2}{2w_{x_2}^2}}}+0.02t-0.01x_2+\epsilon$$

(55)

Assume $x_1\in \left[5,10\right]$ is the process parameter, $x_2\in \left[40,100\right]$ is the shape parameter with $x_2\sim N\left(x_{2_0},w_{x_2}^2\right)$ , $t\in \left[120,180\right]$ is the time parameter with $t\sim N\left(t_0,w_t^2\right)$ , and the random noise is $\epsilon \sim N\left(0,{0.0001}^2\right)$ . Using a maximum-minimum Latin hypercube design, 60 spatial samples of x1 and x2 were randomly generated on [5, 10]× [40, 100] . Additionally, 600 temporal observation points of t were randomly generated on [120, 180] . Different values of (x1, x2 ) and t were randomly combined to generate spatio-temporal data of size 60× 600 . For this data, inference based on the Gaussian process model involves the inverse operation of a 600× 600-dimensional covariance matrix, resulting in an overall computational complexity of $O\left({\left(60\times 600\right)}^3\right)$ . Therefore, this simulation data demonstrates the feasibility and effectiveness of the proposed method.

To investigate the empirical correlation of output levels in both sample and time dimensions within the simulation data, R language was employed to generate and visualize the data: (a) correlation of output levels across samples; (b) empirical correlation of output levels aggregated by time intervals. Figure 3 reveals that, in the spatial dimension, output levels of the same spatial parameter are correlated with each other, and output levels between adjacent spatial parameters are highly correlated. In the temporal dimension, output levels at adjacent time points are correlated with each other, and this correlation weakens as the time interval increases.

To visually demonstrate the dynamic changes of the simulation function with respect to the spatial factor (x1,x2) and the temporal factor t , the surface plot and contour plots shown in Figure 4 were generated. As shown in Figure 4, the output level peaks at $x_1=x_{1_0}$ and decreases with increasing $\left|x-x_{1_0}\right|$ . Along the x2 dimension,the output exhibits a unimodal structure centered at $x_2=x_{2_0}$ , overlaid with a weak negative linear correlation trend overall. Along the t dimension, the output shows a unimodal structure centered at $x_2=x_{2_0}$ , overlaid with a weak positive linear correlation trend overall.

The results in Figure 3 and Figure 4 collectively indicate that the output level is influenced by the combined effects of spatial and temporal factors. Therefore, identifying and quantifying the correlation between dynamic peaks and spatio-temporal factors is crucial for regulating spatial and temporal factors and achieving optimal system response.

4.2.2. Simulation Results

This paper randomly divides 60 sets of simulation data into a training set (50 sets) and a test set (10 sets), with the former used for model fitting and the latter for out-of-sample prediction. Use R language to perform hold-out tests on each data set at 25%, 50%, 75%, and 90% ratios among 600 observation points respectively. Each proportion is randomly repeated three times, and the average is taken. Simulation results for each method are presented in Table 1.

The comparison results in Table 1 show that compared with other models, the NSGP model exhibits the lowest root mean square error and 95% posterior confidence interval length across all test scenarios. Compared to the NSGP model, the 95% posterior confidence interval coverage rate of the “point-in-time LR” model is closer to the nominal 95% level , primarily due to its wider 95% posterior confidence interval length. However, this model's RMSE is higher than NSGP's across all test cases. In summary, the NSGP model constructed in this study demonstrates superior performance in prediction accuracy, interval efficiency, and overall stability.

Figure 5 shows the predicted mean values obtained from direct covariance inversion and the FFBS algorithm. The points represent the output of ${\widehat{v}}_i(s_{1:250}^{\ast })|{\widehat{v}}_i(s_{1:750}^D)$ calculated using equation (21) for the given parameter set $\left(({\sigma }_i^2,{\gamma }_i,{\tau }_i\right)$, and the left panel represents the case when i = 1 and the left panel represents the case when i = 2 . The covariance inversion and FFBS algorithms compute the same quantity, so their predicted means are identical—as shown in Figure 5, the red triangle and blue dots overlap. Although the difference between the two computational methods is negligible, the FFBS algorithm enables precise calculation of numerous statistics, such as the posterior predicted mean, variance, and marginal likelihood function. Furthermore, these two algorithms exhibit significant differences in computational time. As shown in Figure 6, triangles and dots represent the FFBS algorithm and direct computation, respectively. When n=5000, estimating the likelihood function via the FFBS algorithm in R takes less than 1 second, whereas direct computation requires over 100 seconds. Clearly, the FFBS algorithm is more computationally efficient.

4.3. Case Study

Concrete production is a critical phase in dam construction. Concrete raw materials such as cement, coarse aggregates, and fine aggregates are mixed in a mixer or batching plant and then transported to the construction site for pouring. The slump of concrete mixture at the outlet reflects the workability of the concrete mixture in its initial state after mixing, and is an important indicator for evaluating fluidity, compactness, and construction performance. Resistance to segregation is a key indicator of concrete stability during transportation. Strong resistance prevents layering of aggregates and mortar, ensuring mixture uniformity. Mortar loss rate reflects the degree of mortar loss during transport; excessive loss leads to poor workability and insufficient strength in concrete. Therefore, properly controlling concrete slump, segregation resistance, and mortar loss rate is crucial for enhancing project quality.

This paper uses the relevant testing standards for concrete mixtures during the production of a dam as an example. Classified according to quality characteristics, concrete slump is the nominal the best, segregation resistance is the larger the better, and mortar loss rate is the smaller the better. The slump (response y1), segregation resistance (response y2), and mortar loss rate (response y3) of 70 concrete mixes were observed under varying transport and waiting times (time factor $s\in \left[20,40\right]$) . The two parameters in this experiment were water-cement ratio (spatial factor x1= 30%, 35%, 40%, 45%, 50%, 55%) and cement temperature (spatial factor x2= 12, 15, 18, 21, 24, 27, 30, 33).Assume the target value for response y1 is 5cm; the quality compensation resulting from interactions between spatiotemporal factors is 10 yuan/m³ . The objective of this experiment is to achieve optimal concrete strength by selecting appropriate transport and waiting times, water-cement ratio, and cement temperature.

Figure 7 displays the empirical strength correlation of 70 concrete samples across the transport and waiting time range of [20,40] . Generated using R-Language's construction- "Shiny," this figure visually illustrates the associative patterns of strength across samples and the temporal dimension. Results indicate that in the spatial dimension, output levels between adjacent spatial parameters are highly correlated. In the temporal dimension, the average mass strengths at proximate time points are correlated with each other, and this correlation weakens as the time interval increases.

As shown in Figure 7, quality intensity is jointly influenced by spatial and temporal parameters. This study randomly selected 18 samples to test the models under conditions where 25% of the measurement values were excluded from the analysis. To evaluate model performance, a comparative analysis of modeling and prediction capabilities was conducted within the R framework for the NSGP, GP, and SGP methods. The optimization results for each method are detailed in Table 2. The core functions invoked during model construction were consistent across all three approaches, utilizing the “ginv” function from the MASS package, the “dist” and ‘optim’ functions from the stats package, and the “expand.grid” functions from the base package in R.

A comparison of the optimization results in Table 2 clearly shows that the NSGP method achieves RMSE of 0.8847 and a 95% posterior confidence interval length of 1.8845, both lower than those of other methods. The 95% posterior confidence interval coverage rate of the NSGP method is 0.7233, which is closer to 0.95 compared to other methods. The proposed method yields the optimal combined parameter combination for spatio-temporal factors, with an optimized prediction cumulative sum of 150.4114, exceeding the other two methods. The multivariate multidimensional quality gain-loss value obtained by this method is 24.5306, significantly lower than the optimized values of other methods, enabling superior concrete strength design.

5. Conclusions

This paper establishes an engineering robust design scheme based on non-separable Gaussian processes and Quality Gain-Loss functions theory. First, considering the correlation between spatial and temporal factors, the Matérn kernel function is constructed to establish a functional model between temporal factors and quality characteristics. The fast and accurate algorithms of Kalman filters and Kalman smoothing are employed to estimate and predict model data. Second, the gain-loss weights for each quality characteristic were determined based on the signal-to-noise ratio, forming a joint weight matrix. Multivariate, multidimensional quality gain-loss function models were then constructed for the larger the better, smaller the better, and nominal the best, respectively. Finally, the time and spatial factors were optimized to obtain the combined optimal parameter combination. In simulation examples, simulation data were substituted into NSGP, GP, SGP, LR, FR, and NN models for comparative analysis, with the NSGP model yielding superior simulation results. Through calculating gains and losses based on slump at discharge, segregation resistance, mortar loss rate, and relevant spatio-temporal factors during dam concrete production, the optimization results combining NSGP with three quality gain-loss functions — Artiles-leon's multivariate quality loss function, Pignatiello's multivariate quadratic quality loss function, and the traditional multivariate quality gain-loss function — all yielded optimal outcomes. This methodology not only achieves joint optimization of spatiotemporal factors but also reduces reliance on the covariance separability assumption, offering novel insights and technical pathways for hydraulic concrete quality control and robust engineering design.