Three-dimensional Fuzzy Modeling for Nonlinear Distributed Parameter Systems Using Simultaneous Perturbation Stochastic Approximation

1. Introduction

The spatial distribution characteristics are widely present in many manufacturing industries such as petroleum, chemical, steelmaking, and rolling. The dynamic characteristics of the system change with time and space, and this type of system is called a distributed parameter system (DPS)[1]. The dynamic characteristics of DPS are usually described by partial differential equations and integral equations. Traditionally, the spatial distribution characteristics of such systems are often ignored, and the inherent DPS is approximated as lumped parameter systems (LPS), and then are modeled and controlled using the mature lumped parameter methods. With the development of the times, people’s requirements for product quality are gradually increasing. With the advancement of sensors, actuators, and computing technology, the modeling, control, and optimization techniques of DPS have received increasing attention in the fields of science and engineering. For practical use in areas like optimization, control design, analyzing, and numerical simulations of DPS[2], it is necessary to use an appropriate mathematical model to describe its dynamic characteristics. Therefore, the modeling of DPS has become a very important research field.

The modeling methods for DPS can usually be divided into: mechanism modeling with known partial differential equations and grey box/black box modeling with unknown partial differential equations. Mechanism modeling is the process of obtaining partial differential equations for DPS from the point of view of physical or chemical reaction knowledge of the system, such as thermal energy balance, material balance, or momentum balance relationships. Considering the complexity of calculations and the restricted quantity of sensors and actuators, modeling such systems further requires the use of spectral methods[1], finite element[3], finite difference techniques, etc., to simplify the essentially infinite dimensional system into a finite dimensional system[4]. Nonetheless, in real-world scenarios, it will be challenging to obtain accurate mathematical descriptions of the system through mechanism modeling method due to incomplete knowledge of physical and chemical processes. The second modeling method is based on data, including gray box modeling and black box modeling[5]. The model expressed by partial differential equations of the system is established through prior knowledge, and the unknown parameters of the model are obtained by identification from the input/output data of the system. This method is called grey box modeling. Directly identifying the structure and parameters of a system from input/output data is called black box modeling. In practical applications, black box modeling has a wider range of applications. How to extract the relationship between system inputs/outputs with spatial distribution characteristics from data is currently a challenge in modeling research.

In recent decades, researchers have conducted continuous research on black box modeling of DPS. Most methods decompose variables into a series of spatial functions and temporal coefficients[6], and then reconstruct the system. Traditional methods usually pre-set spatial functions, such as Jacobian polynomials, spline functions, etc., and prior knowledge also dictates the quantity of spatial functions. The time coefficients are estimated using traditional modeling methods for LPS. In the past 20 years, spatial functions have been estimated using the Karhunen-Loève (KL) decomposition method[7,8]. The time coefficient is still estimated using traditional modeling methods for LPS. In practical applications, essentially infinite dimensional systems require model simplification to obtain lower order ordinary differential equations. However, model simplification can lead to unmodeled dynamics and unknown nonlinearity[9].

Intelligent methods like neural networks and fuzzy systems have been incorporated into DPS modeling as well. References [8,10,11] combined neural networks with space-time separation methods and achieved good experimental results on parabolic DPS modeling problems. Wang et al. [10] used a nonlinear PCA method to obtain space base functions. After space-time separation and dimension reduction, a neural network was used to identify the time coefficients. Zhang et al.[11] used the PCA algorithm to obtain space base functions. After space-time separation and dimension reduction, the time model was decomposed into linear and nonlinear parts. The ARX model was used to model the linear part, and the RBF neural network was employed to model the nonlinear part. Fan et al.[8] used the KL decomposition method to obtain space base functions. After space-time separation and dimension reduction, a mixed time series model of long short-term memory network and multi-layer perceptron was used to model the time coefficients. However, the aforementioned modeling methods still require dimension reduction.

Three-dimensional fuzzy modeling[9,12] is a new type of DPS modeling method that has developed in the past five years. In three-dimensional fuzzy model, the functions of space-time separation and space-time reconstruction are naturally integrated into three-dimensional fuzzy rules. In each three-dimensional fuzzy rule, space-time separation is naturally achieved, where the preceding component is computed as the time factor and the resulting component is regarded as the space base function. In the combination operation of multiple excited three-dimensional fuzzy rules, space-time reconstruction is naturally achieved. Compared with traditional modeling methods of distributed parameter systems, three-dimensional fuzzy modeling method has two advantages: it does not rely on model dimension reduction and has language interpretability.

At present, the three-dimensional fuzzy model has been applied in the modeling of DPS such as rapid thermal processing and chemical reactors. Zhang et al.[9] introduced a three-dimensional fuzzy modeling approach utilizing clustering and support vector regression, where nearest neighbor clustering (NNC) method and similarity analysis technique were employed for mining and condensing the initial preceding components in the three-dimensional fuzzy rules, and support vector regression machines are employed to acquire knowledge about the space base functions of the resulting components in the three-dimensional fuzzy rules. However, the clustering results of the NNC algorithm are influenced by the initial representative point selection, and the structural simplification process is very cumbersome. Zhang et al.[12] proposed a method for three-dimensional fuzzy modeling that utilizes KL decomposition to derive space base functions for the resulting sets of the three-dimensional fuzzy rule. The resulting sets consist of a restricted amount of space base functions with higher weights. This method pertains to the scope of dimension reduction and results in dynamic characteristics that have not been modeled.

Based on the above discussion, this paper introduces an innovative three-dimensional fuzzy modeling method for DPS. Utilizing the Affinity Propagation (AP) clustering algorithm, the preceding components of the fuzzy rule are derived within the three-dimensional fuzzy modeling framework. By considering every data point as a potential class representative point, the AP algorithm circumvents the issue of clustering outcomes being constrained by the choice of initial class representative points. Unlike existing methods such as KL decomposition or support vector regression to learn space base functions, this paper sets the Fourier space base function as the resulting component of three-dimensional fuzzy rules. The construction of Fourier space base functions is simple and can better approximate the actual spatial distribution characteristics. Furthermore, the parameters of the Fourier space base function are fine-tuned using the simultaneous perturbation stochastic approximation (SPSA) algorithm.

The key innovations of this study are outlined below:

1) The AP clustering approach is utilized to dynamically determine the most suitable quantity of fuzzy rules and construct preceding sets for the three-dimensional fuzzy model.

2) The Fourier space base functions are initially introduced as the resulting sets of the three-dimensional fuzzy model.

3) The SPSA technique is employed to dynamically update the coefficients of Fourier space base functions for the three-dimensional fuzzy model.

This paper is outlined as follows. The second section offers problem description. The third section offers an in-depth explanation of the three-dimensional fuzzy modeling method. The fourth section applies the proposed modeling method to a rapid thermal chemical vapor deposition system. In the final section, the conclusion is given.

2. Problem Description

2.1. Mathematical Description of Distributed Parameter Systems and Its Characteristics

Numerous processes, like the industrial chemical reaction[1], semiconductor production[14], transmission in cure or reflow processes[15] and Lithium-ion Battery Thermal Process[16,17], are spatially distributed in practice. They can be depicted using the subsequent mathematical formula:

$$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle \frac{{\partial }^{m}w(f,t)}{\partial t^m}}={\rm Y}w(f,t)+\hslash \left(w\right)+\lambda b\left(f\right)Q\left(t\right)\hfill \\ l_a\le f\le l_b,t\ge 0\hfill \\ m=1\mathrm{or}2\hfill \end{array}\end{array}$$

(1)

where, $w(f,t)$ represents the output; $\hslash \left(w\right)$ is a nonlinear function; ${\rm Y}$ is linear space differential operator, which includes first or second order spatial derivatives ($\partial /\partial f,{\partial }^2/\partial f^2$), and is dense in Hilbert space; $\lambda$ is a constant vector; $Q=\left[q_1{q}_2\cdots q_{\underline{L}}\right]\in R^{1\times \underline{L}}$ is manipulating input vector (i.e. multiple control sources); $b\left(f\right)={\left[b_1\left(f\right)b_2\left(f\right)\cdots b_{\underline{L}}\left(f\right)\right]}^T\in R^{\underline{L}}$ is a known smooth function about f, which describes the distribution of Q in the spatial domain $\overline{F}\subset \left[l_a,l_b\right]$.

Usually, boundary and initial conditions also restrict the system depicted in the Equation (1). As an illustration, at t=1, the boundary and initial conditions can be depicted in the following equations:

$$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle \frac{\partial w\left(l_a,t\right)}{\partial f}}=g_a\left(w\left(l_a,t\right),Q\left(t\right)\right)\hfill \\ {\displaystyle \frac{\partial w\left(l_b,t\right)}{\partial f}}=g_b\left(w\left(l_b,t\right),Q\left(t\right)\right)\hfill \end{array}\end{array}$$

(2)

$$\begin{array}{c}\hfill \begin{array}{c}w(f,0)=w_0\left(f\right)\hfill \end{array}\end{array}$$

(3)

where, $g_a\left(w\left(l_a,t\right),Q\left(t\right)\right)$ is dependent on $w\left(l_a,t\right)$ and $Q\left(t\right)$, and $g_b\left(w\left(l_b,t\right),Q\left(t\right)\right)$ is dependent on $w\left(l_b,t\right)$ and $Q\left(t\right)$; $w_0\left(f\right)$ is a known function.

This type of system has spatiotemporal coupling characteristics and is an infinite dimensional system. $w(f,t)$ can be spatiotemporally separated into an infinite weighted sum of space base functions and time coefficients, as shown in Equation (4).

$$\begin{array}{c}\hfill \begin{array}{c}w(f,t)={\sum }_{i=1}^{\infty }{w}_i{\phi }_i\left(f\right)\hfill \end{array}\end{array}$$

(4)

Practically, it is necessary to utilize a finite quantity of sensors positioned at $f_1,f_2,\cdots ,f_p$; therefore, $\overline{F}=\left[f_1,f_2,\cdots ,f_P\right]$. $w(\overline{F},t)=\left[w\left(f_1,t\right)\right.,\left.w\left(f_2,t\right),\cdots ,w\left(f_p,t\right)\right]$ represents the output in space. The data-driven modeling challenge for the DPS, as outlined in Equations (1) - (3), entails the task of determining a spatiotemporal coupled model using the input ${\left\{Q\left(k\right)\right\}}_{k=1}^K$ and the output ${\left\{w\left(f_i,k\right)\right\}}_{i=1,k=1}^{P,K}$, where K is the time length. Essentially, systems with an infinite number of dimensions will be converted into the systems with a finite number of dimensions. When the mechanism model described in Equations (1) - (3) is unknown, utilizing the KL decomposition technique, the system’s space base function can be derived from the gathered data set $\overline{S}={\left\{Q\left(k\right)w\left(f_i,k\right)\right\}}_{i=1,k=1}^{P,K}$. Then, $w(f,t)$ can be approximately depicted as Equation (5).

$$\begin{array}{c}\hfill \begin{array}{c}{w}_n(f,t)=\sum _{i=1}^n{w}_i{\phi }_i\left(f\right)\hfill \end{array}\end{array}$$

(5)

The time coefficient in Equation (5) can be identified using traditional methods.

However, the modeling approach of DPS using KL decomposition unavoidably results in neglected dynamic features and uncertainties in the identified model due to the decrease of dimensions.

2.2. Three-dimensional Fuzzy Modeling

The three-dimensional fuzzy system[18] was first established in 2007 and has since been used in the domain of DPS control[19]. Recently, there has been effective development of the three-dimensional fuzzy system for modeling DPS. The relevant theory of three-dimensional fuzzy system and black box modeling theory were combined to form a three-dimensional fuzzy modeling framework as shown in Figure 1[9]. Unlike traditional modeling methods, the three-dimensional fuzzy modeling method can model a DPS without reduction dimension.

Assume that K and J represent the order of the input variable $Q\left(k\right)$ and output variable $w(\overline{F},k)$, respectively. Input variables of three-dimensional fuzzy system are $\left\{Q\right(k-1),Q(k-2),\cdots ,Q(k-K\left)\right\}$and $\{w(\overline{F},k),w(\overline{F},k-1),\cdots ,w(\overline{F},k-J)\}$. The three-dimensional fuzzy rule is given as shown in Equation (6).

$$\begin{array}{c}\hfill \begin{array}{c}{\overline{V}}_l:\hfill \\ \mathrm{if}w(\overline{F},k-1)\mathrm{is}{\overline{O}}_{1l}\mathrm{and}\cdots \mathrm{and}w(\overline{F},k-J)\mathrm{is}{\overline{O}}_{Jl}\hfill \\ \mathrm{and}{q}_1(k-1)\mathrm{is}{Q}_1^{1l}\mathrm{and}\cdots \mathrm{and}{q}_1(k-K)\mathrm{is}{Q}_1^{Kl}\hfill \\ ⋮\hfill \\ \mathrm{and}{q}_{\underline{L}}(k-1)\mathrm{is}{Q}_{\underline{L}}^{1l}\mathrm{and}\cdots \mathrm{and}{q}_{\underline{L}}(k-K)\mathrm{is}{Q}_{\underline{L}}^{Kl}\hfill \\ \mathrm{then}w(\overline{F},k)\mathrm{is}{\phi }_l\left(\overline{F}\right)\hfill \end{array}\end{array}$$

(6)

where ${\phi }_l\left(\overline{F}\right)$ is a space base function, ${\overline{O}}_{sl}$ is a three-dimensional fuzzy set, $Q_g^{hl}$ is a conventional fuzzy set, $s=1,2,\dots ,J,h=1,2,\dots ,K,g=1,2,\dots ,\underline{L},l=1,2,\dots ,\overline{Z}$ , $\overline{Z}$ indicates the count of fuzzy rules.

For any three-dimensional fuzzy rule ${\overline{V}}_l(l=1,2,\dots ,\overline{Z})$ , its preceding set ‘$\mathrm{if}w(\overline{F},k-1)is{\overline{O}}_{1l}\mathrm{and}$$\cdots \mathrm{and}{q}_{\underline{L}}(k-K)\mathrm{is}{Q}_{\underline{L}}^{Kl}$’ corresponds to the time coefficient $w_l$. The resulting set ‘$\mathrm{then}w(\overline{F},t)\mathrm{is}{\phi }_l\left(\overline{F}\right)$’ corresponds to the space base function ${\phi }_l\left(F\right)$, and space-time separation is ingeniously achieved through three-dimensional fuzzy rules. If $Z^{\prime }$ three-dimensional fuzzy rules are excited, the time coefficients and space base functions of $Z^{\prime }$ three-dimensional fuzzy rules will naturally realize space-time reconstruction. The final result is ${\sum }_{l=1}^{Z^{\prime }}{w}_l{\phi }_l\left(\overline{F}\right)$.

The three-dimensional fuzzy modeling framework can be divided into three parts.

(1) Acquire the preceding set of fuzzy rules in order to derive temporal coefficients.

(2) Acquire knowledge about the space base functions of the resulting set.

(3) Obtain predicted output through space-time reconstruction.

The study will explore an innovative approach to three-dimensional fuzzy modeling within this specific framework.

3. SPSA Learning Based Three-Dimensional Fuzzy Modeling

3.1. Modeling Methodology

As for an unknown DPS system, space base functions are frequently calculated by KL decomposition technique. However, implementing the KL decomposition technique is bound to result in the simplification of models. The three-dimensional fuzzy model has the innate ability to amalgamate space-time separation and space-time reconstruction into three-dimensional fuzzy rules and thereafter avoids the model reduction.

In this study, a novel SPSA learning based three-dimensional fuzzy model methodology is proposed, the scheme of which is shown in Figure 2.

The modeling methodology is divided into two parts, namely, learning for preceding components of three-dimensional fuzzy rules and learning for resulting components of three-dimensional fuzzy rules. In the first part, by using AP clustering, preceding components of three-dimensional fuzzy rules are learned. In the second part, the Fourier function is initially selected as the space base function as the resulting components of the three-dimensional fuzzy rule. Then the SPSA method is used to determine the coefficients of the Fourier space base function adaptively. Three distinct characteristics are evident in the SPSA learning based three-dimensional fuzzy modeling method as blow.

1) AP clustering is employed to manage a collection of spatiotemporal data and extract the preceding components of three-dimensional fuzzy rules in an adaptable manner, without the need to define the number of clusters.

2) The Fourier space base functions are investigated in three-dimensional fuzzy modeling.

3) SPSA algorithm is used to acquire optimal parameters of Fourier space base functions.

3.2. AP Clustering Learning for Preceding Components of Three-dimensional Fuzzy Rule

Frey and Dueck introduced Affinity Propagation (AP) Clustering[13], a novel clustering technique, in 2007. It can handle with large-scale data sets in a short period of time without initializations of cluster number and cluster core. The process of AP clustering involves initially selecting specific samples known as exemplars, followed by linking each sample on the left to its closest examples. Two kinds of message, namely responsibility and availability, are passed in AP clustering. Let $V={\left[r(i,\nu )\right]}_{n\times n}$ be responsibility matrix, $A={\left[a(i,\nu )\right]}_{n\times n}$ be the availability matrix, and $S={\left[s(i,\nu )\right]}_{n\times n}$ be the similarity matrix. The responsibility $r(i,\nu )$ reflects how strongly the sample i want to choose the candidate exemplar $\nu$ as its exemplar. The presence of $a(i,\nu )$ mirrors the collective data demonstrating the appropriateness of sample i in selecting candidate $\nu$ as the exemplar. The similarity $s(i,\nu )$ represents the degree of similarity between the sample i and its nearest exemplar $\nu$. Then, the message between V and A is propagated so that the sample similarity S in the cluster is the largest and the sample similarity between the clusters is the smallest. The detailed algorithm description is presented in Appendix A and the algorithm procedure is shown in Figure A.1 of Appendix A.

For a spatiotemporal data set $\overline{S}$ from DPS, AP clustering algorithm is employed to cluster the dataset $\overline{S}$ to $Z^{\prime }$ clusters. Each cluster corresponds to the preceding components of a three-dimensional fuzzy rule. Hence, the preceding components of $Z^{\prime }$ three-dimensional fuzzy rules are achieved. Assumed that $o^1,\cdots ,o^{Z^{\prime }}$ with $o^i=\left(o_{\chi }^i,o_c^i\right)$ are $Z^{\prime }$ cluster cores. These cores will guide to construct the preceding components of $Z^{\prime }$ three-dimensional fuzzy rules. Utilizing a Gaussian type membership function, the cluster’s core $o^i$ aligns with the Gaussian type MFs’ center in the preceding components, while $o_{\chi }^i$ and $o_c^i$ represent the centers of three-dimensional MFs and conventional MFs, respectively.

The breadth of Gaussian type membership function can be computed according to the interval range of the input variables from the input data set. Regarding the spatial output variable $w(\overline{F},k)$, the breadth of the Gaussian type three-dimensional MF at an identical spatial point is calculated as in (7).

$$\begin{array}{c}\hfill \begin{array}{c}{\sigma }_w\left(f_j\right)\approx max\left({\displaystyle \frac{w^{max}\left(f_j\right)-w^{min}\left(f_j\right)}{\eta }}\right)\hfill \\ \hfill \eta \in [0.01,2]\hfill \end{array}\end{array}$$

(7)

where $w^{max}\left(f_j\right)$ and $w^{min}\left(f_j\right)$ represent the highest and lowest $w(\overline{F},k)$ values at the sensor’s location $f=f_j$; $\eta$ is a variable set by the users.

Regarding the variable $Q\left(k\right)$, the breadth of the conventional Gaussian MF is outlined below.

$$\begin{array}{c}\hfill \begin{array}{c}{\sigma }_Q=\left[\begin{array}{ccc}{\sigma }_{{\mu }_1}\hfill & \cdots \hfill & {\sigma }_{{\mu }_{\underline{L}}}\hfill \end{array}\right]\hfill \end{array}\end{array}$$

(8)

where ${\sigma }_{{\mu }_i}={\displaystyle \frac{q_i^{max}-q_i^{min}}{{\eta }_i}}$, ${\eta }_i\in [0.01,2]$; $q_i^{max}$and $q_i^{min}$ represent the highest and lowest values of the ith input variable; ${\eta }_i$ is a variable set by the users; $i=1,2,\cdots ,\underline{L}$.

The coefficients $\eta$ in Equation (7) and ${\eta }_i$ in Equation (8) can be modified based on the present situation.

3.3. SPSA Learning for Resulting Components of Three-dimensional Fuzzy Rule

3.3.1. Fourier Space Base Function

Space base function is crucial to the performance of the modeling[6]. In three-dimensional fuzzy modeling, KL decomposition technique and SVR have been used to estimate space base functions. Traditional DPS modeling approaches allow for the selection of fixed form functions, including Spline functions, Legendre polynomials, and Jacobi polynomials[6], as space base functions, in addition to KL decomposition estimation. In this study, Fourier space base functions are taken as the resulting components of three-dimensional fuzzy rules. Spatial dynamics of a DPS can be represented by multiple Sine waves and Cosine waves with different amplitude and frequency. A Fourier space base function is expressed as in (9).

$$\begin{array}{c}\hfill \begin{array}{c}{\phi }_l\left(\overline{F}\right)=a_{i}sin\left(b_l\overline{F}\right)+c_{l}cos\left(b_l\overline{F}\right)+d_l\hfill \\ \hfill (l=1,2,\dots ,Z^{\prime })\hfill \end{array}\end{array}$$

(9)

Therefore, the variable $w(\overline{F},k)$ can be expressed as a linear sum of Fourier space base functions, with each space base function multiplied by its corresponding time coefficient, which is shown as in (10).

$$\begin{array}{c}\hfill \begin{array}{c}w(\overline{F},k)=\sum _{l=1}^{Z^{\prime }}{\psi }^l\left(k\right){\phi }_l\left(\overline{F}\right)\hfill \end{array}\end{array}$$

(10)

$$\begin{array}{c}\hfill \begin{array}{c}{\psi }^l\left(k\right)=\sum _{h=1}^P{a}_h^l\prod _{s=1}^J{E}_A\times \prod _{i=1}^{\underline{L}}\prod _{g=1}^K{E}_B\\ E_A=exp\left(-{\left(\left(w\left(f_h,k-s\right)-c_{sh}^l\right)/{\sigma }_{sh}^l\right)}^2\right)\\ E_B=exp\left(-{\left(\left(q_i(k-g)-d_i^{gl}\right)/{\delta }_i^{gl}\right)}^2\right)\end{array}\end{array}$$

(11)

where ${\psi }^l\left(t\right)$ is called as fuzzy basis function[19] in the context of three-dimensional fuzzy system, also called time coefficient in a three-dimensional fuzzy model; $c_{sh}^l$ and ${\sigma }_{sh}^l$ are denoted as central point and breadth of the Gaussian type three-dimensional fuzzy set ${\overline{O}}_{sl}$ at the hth sensor position; $d_i^{gl}$ and ${\delta }_i^{gl}$ represent the central point and breadth of the conventional Gaussian fuzzy set $Q_i^{gl}$; $a_h^l$ represents the weight assigned to the hth sensor location.

The SPSA method will be used to learn the coefficients of the Fourier space base function.

3.3.2. Parameter Learning Using SPSA Algorithm

The SPSA [20] is a highly effective technique for optimizing complex nonlinear problems. It possesses a robust theoretical basis and is frequently more effective and simpler to execute compared to widely recognized techniques like simulated annealing or genetic algorithms.In contrast with gradient descent method, SPSA does not necessitate the direct measurement of the gradient of the target function. Instead, it just requires two measurements of the target function at every iteration to approximate the gradient, regardless of the complexity of the optimization problem.The two measurements are obtained by randomly and simultaneously altering all the variables in the problem in a correct manner, which is known as “simultaneous perturbation”.

In this study, in each iteration the approximation gradients of all parameters of Fourier space base functions are computed by two randomly disturbed measurements of target functions. Let ${\widehat{\vartheta }}_k=[a_{1,k},b_{1,k},c_{1,k},d_{1,k},a_{2,k},b_{2,k},c_{2,k},d_{2,k},\dots ,a_{Z^{\prime },k,},b_{Z^{\prime },k},c_{Z^{\prime },k},d_{Z^{\prime },k}]\in R^{{4Z}^{\prime }}$ be all the parameters at the kth iteration, where $a_{i,k},b_{i,k},c_{i,k}$ and $d_{i,k}$ are parameters of ith space base function at the kth iteration. The target function is described as the L2 norm of the difference between the spatial output $w(\overline{F},t)$ and the predicted spatial output $\widehat{w}(\overline{F},t)$ as follows.

$$\begin{array}{c}\hfill \begin{array}{c}J\left(\widehat{\vartheta }\right)={\displaystyle \frac{1}{2}}\sum _{i=1}^n{\left(w_i(\overline{F},t)-{\overline{w}}_i(\overline{F},t)\right)}^2\hfill \end{array}\end{array}$$

(12)

The algorithm of SPSA is described as follows.

Step 1: Initialization and selection of coefficients

Initialize the timer number as $k=0$. Select an initial estimate ${\widehat{\vartheta }}_0$ and choose non-negative constants b, d, A, $\alpha$, and $\beta$ for SPSA algorithm. Then, calculate sequences $b_k$ and $d_k$ as in (13) and (14).

$$\begin{array}{c}\hfill \begin{array}{c}{b}_k=b/{(A+k+1)}^{\alpha }\hfill \end{array}\end{array}$$

(13)

$$\begin{array}{c}\hfill \begin{array}{c}{d}_k=d/{(k+1)}^{\beta }\hfill \end{array}\end{array}$$

(14)

Step 2: Simultaneous perturbation vector generation

Create a random perturbation vector ${\Delta }_k$ in a $4Z^{\prime }$-dimensional space, where each component of ${\Delta }_k$ is individually produced from a zero-mean probability distribution that satisfies the specified constraints in [20]. Here, we use Bernoulli $\pm 1$ distribution, which means that there is an equal probability of 0.5 for each event of $\pm 1$.

Step 3: Evaluation of the cost function

Calculate two values of the cost function $J\left(\widehat{\vartheta }\right)$ by perturbing the current ${\widehat{\vartheta }}_k$ simultaneously: $J\left({\widehat{\vartheta }}_k+d_k{\Delta }_k\right)$ and $J\left({\widehat{\vartheta }}_k-d_k{\Delta }_k\right)$, where $d_k$ and ${\Delta }_k$ are obtained from Step 1 and Step 2, respectively.

Step 4: Approximation of the gradient

Compute the simultaneous perturbation estimate for the gradient $\mathbf{g}\left({\widehat{\vartheta }}_k\right)$:

$$\begin{array}{c}\hfill \begin{array}{c}{\widehat{\mathbf{g}}}_k\left({\widehat{\vartheta }}_k\right)={\displaystyle \frac{J\left(\widehat{\vartheta }k+d_k{\Delta }_k\right)-J\left({\widehat{\vartheta }}_k-d_k{\Delta }_k\right)}{2d_k}}\left[\begin{array}{c}{{\Delta }_{k1}}^{-1}\\ {{\Delta }_{k2}}^{-1}\\ ⋮\\ {{\Delta }_{kp}}^{-1}\end{array}\right]\hfill \end{array}\end{array}$$

(15)

where ${\Delta }_{ki}$ represents the ith element of the ${\Delta }_k$ vector.

Step 5: Revision of the $\vartheta$ estimation

$$\begin{array}{c}\hfill \begin{array}{c}{\widehat{\vartheta }}_{k+1}={\widehat{\vartheta }}_k-b_k{\widehat{\mathbf{g}}}_k\left({\widehat{\vartheta }}_k\right)\hfill \end{array}\end{array}$$

(16)

Using (16) to modify ${\widehat{\vartheta }}_k$ to a fresh value ${\widehat{\vartheta }}_{k+1}$.

Step 6: Iterating or terminating

Go back to Step 2 and substitute k with $k+1$. Terminate the algorithm if there is minimal variation in several consecutive iterations or if the maximum number of iterations has been reached.

3.3.3. Three-dimensional Fuzzy Modeling Flowchart

The proposed three-dimensional fuzzy modeling is designed in seven stages, as illustrated in Figure 3.

Stage 1: Set up the elements of a three-dimensional fuzzy model, consisting of the structure of three-dimensional fuzzy rules, the method of three-dimensional fuzzification, the method of defuzzification, various types of fuzzy operators, and other related components.

Stage 2: Assign the Fourier space base functions, as depicted in Equation (9), to the resulting components of the three-dimensional fuzzy rules.

Stage 3: Create a collection of spatial and temporal input-output data from a DPS.

Stage 4: Cluster the spatiotemporal data using AP clustering algorithm.

Stage 5: Create the preceding components of three-dimensional fuzzy rules based on the clusters identified in Stage 4.

Stage 6: Utilizing the SPSA technique to determine the parameters of the Fourier space base functions in resulting components of three-dimensional fuzzy rules.

Stage 7: Combine the preceding components from Stage 5 with the resulting components of three-dimensional fuzzy rules to form the entire three-dimensional fuzzy rule foundation.

A three-dimensional fuzzy model is obtained after the completion of a seven-stage design process.

4. Application to RTCVD System

4.1. RTCVD System

Figure 4 displays the schematic of the rapid thermal chemical vapor deposition (RTCVD) system.

A 6-inch silicon wafer is positioned on a revolving support inside the system’s chamber and subjected to heat from a heating system. The heating system consists of three Lamp banks. Lamp bank 1 uniformly warms the whole surface area of the wafer, Lamp bank 2 specifically heats the periphery of the wafer, and Lamp bank 3 equally heats the wafer as a whole. The Figure 5 displays the incident radiation flux emitted by the heating lights. A 10% concentration of silane gas (SiH4) is fed into the reactor. SiH4 undergoes decomposition to produce silicon (Si) and hydrogen gas (H2). A layer of polysilicon, with a thickness of 0.5$\mu m$, is applied onto the wafer by depositing it at temperatures of about 800K or more, which takes about 1 minute. During wafer processing, the support is turned to ensure that the temperature is evenly distributed in the azimuthal direction. Due to the thinness of the silicon wafer, all changes in temperature along the azimuth direction are disregarded. Hence, it is necessary to maintain consistent temperature distribution over the wafer radius in order to achieve a homogeneous deposition of polysilicon on the wafer. This may be achieved by regulating the power supplied to the three zones of lighting groups. A one-dimensional space model of the thermal dynamics may be expressed as the following PDE [21,22]:

$$\begin{array}{c}\hfill \begin{array}{c}\begin{array}{cc}\hfill \partial T_f^{\prime }/\partial k^{\prime }=& {\kappa }_0\left[\left(1/f^{\prime }\right)\partial T_f^{\prime }/\partial f^{\prime }+{\partial }^2{T}_f^{\prime }/\partial f^{\prime 2}\right]\hfill \\ \hfill & +{\sigma }_0\left(1-T_f^{\prime 4}\right)+{\omega }_f{\varrho }_1\left(f^{\prime }\right)q_1\hfill \\ \hfill & +{\omega }_f{\varrho }_2\left(f^{\prime }\right)q_2+{\omega }_f{\varrho }_3\left(f^{\prime }\right)q_3\hfill \end{array}\hfill \end{array}\end{array}$$

(17)

Constrained by the specified border conditions.

$$\begin{array}{c}\hfill \begin{array}{c}\left\{\begin{array}{cc}\partial T_f^{\prime }/\partial f^{\prime }={\sigma }_{ed}\left(1-T_f^{\prime 4}\right)+{\varrho }_{ed}{q}_2\hfill & \mathrm{when}{f}^{\prime }=1\hfill \\ \partial T_f^{\prime }/\partial f^{\prime }=0\hfill & \mathrm{when}{f}^{\prime }=0\hfill \end{array}\right.\hfill \end{array}\end{array}$$

(18)

where the non-dimensional wafer temperature, denoted by $T_f^{\prime }=T_f/T_{amb}$, represents the ratio of the actual wafer temperature $\left(T_f\right)$ to the ambient temperature $\left(T_{amb}\right)$, which is 300K. The non-dimensional time, denoted by $k^{\prime }=k/\tau$, represents the ratio of the actual time $\left(k\right)$ to $\tau$, which is $2.9s$. The non-dimensional radius position, denoted by $f^{\prime }=f/V_f$, represents the ratio of the actual radius position $\left(f\right)$ to the wafer radius $\left(V_f\right)$, which is $7.6cm$. The variables $q_1$, $q_2$ and $q_3$ represent the percentage of the light source power. The variables ${\varrho }_1\left(f^{\prime }\right)$, ${\varrho }_2\left(f^{\prime }\right)$, ${\varrho }_3\left(f^{\prime }\right)$ represent the wafer incident radiation flux from the three-zone warming lights, respectively. The parameters in (17) through (18) are listed as follows:

$$\begin{array}{c}\hfill \begin{array}{c}{\kappa }_0=0.0021,{\sigma }_0=0.0012,{\sigma }_{ed}=0.0037\hfill \\ {\varrho }_{ed}=4.022,{\omega }_f=0.0256\hfill \end{array}\end{array}$$

The RTCVD system is a complex system that exhibits both spatial and temporal characteristics and has an unlimited number of dimensions. For practical purposes, a certain number of sensors must be used. Assume that sensors positioned at $f_1,f_2,\cdots ,f_P$ are employed to gauge the system’s output. Consider $\overline{F}$ as the spatial domain, where $\overline{F}$ is a vector containing the elements $f_1,f_2,\cdots ,f_P$. Let $w(\overline{F},k)$ be a vector denoted by $w\left(f_1,k\right),w\left(f_2,k\right),\cdots ,w\left(f_P,k\right)$, where $w(f,k)=T_f(f,k)$. The temporal input $Q\left(k\right)$ is defined as $\left[q_1\left(k\right)q_2\left(k\right)q_3\left(k\right)\right]$, where $k\in R$ is the time variable. The objective of the modeling problem is to determine a spatiotemporal model based on the input ${\left\{Q\left(k\right)\right\}}_{k=1}^K$ and the output ${\left\{w\left(f_i,k\right)\right\}}_{i=1,k=1}^{P,K}$, where K represents the duration of time.

In practical operation, 11 sensors are placed evenly along the radial direction. To simulate the impact of noise, 11 independent sets of Gaussian white noise signals were added to the data collected from 11 sensors. The average value is zero and the standard deviation is $\sigma \left(f\right)=A_d\left(f\right)n_d$, $A_d\left(f\right)=(max\left(w(f,t)\right)-min\left(w(f,t)\right))/3$ , $n_d=0.2\%$.

4.2. SPSA Learning Based Three-dimensional Fuzzy Modeling

In order to obtain sufficient dynamic information of the system, interference signals with amplitudes not exceeding 10% are added to the manipulated input variables $q_1$,$q_2$ , and $q_3$, respectively. Thus, the three manipulated input variables affected by interference signals, namely the excitation signals, can be given by the following Equation (19).

$$\begin{array}{c}\hfill \begin{array}{c}{q}_1\left(t\right)=q_{1a}+{\Delta }_{1a}*q_{1a}*rand\left(\right)*S_r\\ q_2\left(t\right)=q_{2a}+{\Delta }_{2a}*q_{2a}*rand\left(\right)*S_r\\ q_3\left(t\right)=q_{3a}+{\Delta }_{3a}*q_{3a}*rand\left(\right)*S_r\\ S_r=sign(rand\left(\right)-0.5)\end{array}\end{array}$$

(19)

where $q_{1a}$, $q_{2a}$ and $q_{3a}$ are steady-state inputs at a furnace temperature of $1000K$. ${\Delta }_{1a}$,${\Delta }_{2a}$ and ${\Delta }_{3a}$ are the disturbance amplitude of $q_{1a}$, $q_{2a}$ and $q_{3a}$. In this study, $q_{1a}$, $q_{2a}$ and $q_{3a}$ are 0.2028, 0.1008 and 0.2245, respectively, ${\Delta }_{1a}$,${\Delta }_{2a}$ and ${\Delta }_{3a}$ are set to $10\%$.

A sample interval of 0.5 seconds was established in this investigation. The entire simulation time period was 7000 seconds. Therefore, 14000 samples were generated, in which 600 samples were randomly selected for training experiments and then 300 samples were randomly selected for testing experiments.

In this study, for simplicity, $J=1$ and $K=1$ were chosen. Thus, the dataset can be written as Equation (20).

$$\begin{array}{c}\hfill \begin{array}{c}S=\left\{{\chi }^k,w_f^k\right\}=\left\{\left({\chi }_f^k,{\chi }_q^k\right),w_f^k\mid {\chi }_f^k\in R^{11\times 1},\right.\hfill \\ \left.{\chi }_q^k\in R^{1\times 3},w_f^k\in {\mathrm{R}}^{11\times 1},k=1,\cdots ,600\right\}\hfill \end{array}\end{array}$$

(20)

where

$$\begin{array}{c}\hfill \begin{array}{c}{\chi }^k=\left({\chi }_f^k,{\chi }_q^k\right)\hfill \\ {\chi }_f^k=\left[w(\overline{F},k-1)\right]\hfill \\ {\chi }_q^k=\left[q_1(k-1),q_2(k-1),q_3(k-1)\right]\hfill \\ w_f^k=w(\overline{F},k)={\left[w\left(f_1,k\right),\cdots ,w\left(f_{11},k\right)\right]}^T\hfill \end{array}\end{array}$$

For ${\chi }^k$ in the dataset S, AP clustering is used. The AP method sets the coefficient of dampness $\lambda$ to 0.9 and the reference degree P to 0.5. Following the process of clustering learning, a total of Seven fuzzy rules and their accompanying clustering cores were successfully acquired. Equation (7) shows the breadth setting of the Gaussian spatial membership function, whereas Equation (8) shows the breadth setting of the conventional Gaussian membership function.

This paper uses Fourier space base functions to represent the fuzzy rule resulting set, and applies the SPSA method to improve the coefficients of Fourier space base functions.

In the SPSA algorithm, the relevant parameter settings are as follows: $\alpha =0.602$, $\beta =0.101$, $ite{r}_{max}=50000$, $A=5000$, $d=0.003$, $b=0.015$.

Based on the learning results of the AP algorithm, it is evident that there are a total of 7 fuzzy rules, so the parameter dimension that needs to be adjusted is 28, and the optimization range control of the parameters is $[-10,10]$. The end point criterion is defined as either reaching the maximum number of iterations or when the change in iteration becomes less than 0.00001 for five consecutive occasions.

The initial value of ${\widehat{\vartheta }}_0$ is obtained through gradient descent algorithm, as shown in Equation (21).

$$\begin{array}{c}\hfill \begin{array}{cccc}\hfill [0.0376& \hfill 29.4590& \hfill -0.0307& \hfill -1.1126\\ \hfill -0.0343& \hfill 26.6246& \hfill -0.0037& \hfill 0.7921\\ \hfill -0.0272& \hfill 26.6246& \hfill -0.0464& \hfill 0.1736\\ \hfill 0.0048& \hfill 25.0165& \hfill 0.0339& \hfill 0.4068\\ \hfill 0.0007& \hfill 18.4397& \hfill -0.0079& \hfill 0.9764\\ \hfill 0.0016& \hfill 18.4623& \hfill -0.0041& \hfill -0.0672\\ \hfill -0.0323& \hfill -8.0553& \hfill -0.0381& \hfill 0.2669]\end{array}\end{array}$$

(21)

By using the SPSA algorithm, the optimal parameters for the Fourier space base function of the fuzzy rule resulting set can be obtained.

For the convenience of observation, Figure 6 shows seven space base functions.

The seven three-dimensional fuzzy rules derived from the data are listed as follows:

${\overline{\mathbf{V}}}_1:$ if $w(\overline{F},k-1)$ is very Positive Large and $q_1(k-1)$ is Positive Small and $q_2(k-1)$ is Positive Large and $q_3(k-1)$ is very Positive Large, then $w(\overline{F},k)$ is [0.8147 0.9058 0.1270 0.9134 0.6324 0.0975 0.2785 0.5469 0.9575 0.9649 0.1576].

${\overline{\mathbf{V}}}_2:$ if $w(\overline{F},k-1)$ is more than Zero and $q_1(k-1)$ is less than Positive Medium and $q_2(k-1)$ is Positive Small and $q_3(k-1)$ is Positive Medium, then $w(\overline{F},k)$ is [6.7941 6.7002 3.3976 5.6020 0.9932 2.9523 6.4101 5.5455 6.7164 4.5902 0.2500].

${\overline{\mathbf{V}}}_3:$ if $w(\overline{F},k-1)$ is less than Positive Small and $q_1(k-1)$ is more than Postive Medium and $q_2(k-1)$ is more than Zero and $q_3(k-1)$ is very Positive Large, then $w(\overline{F},k)$ is [12.7369 -14.0099 -10.1810 -11.3661 -11.1470 -5.8834 -9.8322 -2.5678 -10.5907 -0.4775 -4.1538].

${\overline{\mathbf{V}}}_4:$ if $w(\overline{F},k-1)$ is very Positive Large and $q_1(k-1)$ is Positive Large and $q_2(k-1)$ is slightly Positive Small and $q_3(k-1)$ is Positive Large, then $w(\overline{F},k)$ is [0.4617 0.9713 8.2346 6.9483 3.1710 9.5022 0.3445 4.3874 3.8156 7.6552 7.9520].

${\overline{\mathbf{V}}}_5:$ if $w(\overline{F},k-1)$ is more than Positive Medium and $q_1(k-1)$ is Positive Large and $q_2(k-1)$ is more than Zero and $q_3(k-1)$ is Positive Medium, then $w(\overline{F},k)$ is [1.8687 4.8976 4.4559 6.4631 7.0936 7.5469 2.7603 6.7970 6.5510 1.6261 1.1900].

${\overline{\mathbf{V}}}_6:$ if $w(\overline{F},k-1)$ is more than Positive Medium and $q_1(k-1)$ is Positive Medium and $q_2(k-1)$ is slightly Positive Small and $q_3(k-1)$ is very Positive Large, then $w(\overline{F},k)$ is [-2.4918 -4.7987 -1.7019 -2.9263 -1.1191 -3.7563 -1.2755 -2.5298 -2.5673 -2.3571 -3.125].

${\overline{\mathbf{V}}}_7:$ if $w(\overline{F},k-1)$ is Positive Large and $q_1(k-1)$ is Positive Small and $q_2(k-1)$ is Positive Small and $q_3(k-1)$ is less than Positive Medium, then $w(\overline{F},k)$ is [-1.2263 -0.4110 -1.1889 -0.7443 -1.0026 -0.4407 -1.3535 -1.6484 -2.6840 -1.4376 -0.6583].

Finally, the complete three-dimensional fuzzy modeling of RTCVD has been achieved. The abbreviation for the suggested three-dimensional modeling technique, which combines AP clustering with Fourier space base functions, is AP Fourier SPSA-3D fuzzy modeling.

The simulation results of AP-Fourier-SPSA-3D are shown in the following figures. The predicted result of the model in the training data is displayed in Figure 7, while the discrepancies between the expected and actual values are illustrated in Figure 8. Figure 9 displays the anticipated outcome of the procedure in the test data, whereas Figure 10 illustrates the flaws in the predictions. The root mean square error (RMSE) performance indices are shown in Table 1.

4.3. Simulation Comparison

In this study, three modeling methods were selected for comparison: KL-LS[23] and NNC-SVR-3D[9]. The KL-LS modeling approach utilizes the KL decomposition technique to compute space functions and decrease model dimensions. Subsequently, the least squares identification method is employed to determine low dimensional temporal coefficients. NNC-SVR-3D fuzzy modeling method is to use nearest neighbor clustering algorithm to learn the preceding components of three-dimensional fuzzy rules, which is reduced with the use of a similarity assessing approach, and use multiple single output support vector machines to learn the space base functions of the resulting components. Additionally, gradient descent approach is utilized to determine the coefficients of Fourier space base function for the suggested three-dimensional fuzzy modeling based on AP clustering and Fourier space base function, which is abbreviated as AP-Fourier GD-3D.

The prediction errors of KL-LS model, NNC-SVR model on training data and test data can be seen in Figure 11 and Figure 12, Figure 13 and Figure 14, respectively. The time normalized absolute errors $TNAE\left(f\right)=\sum \left|e\right(f,k\left)\right|/\sum \Delta k$ of these models on training and test data are shown in Figure 15 and Figure 16, respectively. The RMSE performance indexes of these models are shown in Table 1.

From Table 1 and Figure 7 to Figure 16, it can be evident that the three-dimensional fuzzy models outperform the KL-LS model. The KL-LS modeling method is traditional modeling methods based on KL decomposition, which require dimension reduction. Consequently, this will unavoidably result in the emergence of unmodeled dynamic characteristics and uncertainties. The suggested three-dimensional modeling approach seamlessly merges space-time separation and space-time reconstruction, eliminating the necessity for reducing dimensions. In three three-dimensional fuzzy models, the proposed AP-Fourier SPSA-3D model has better performance. It discloses that the Fourier space base function has better spatial characteristics to disclose the spatial relationship of the DPS system, the AP clustering algorithm can adaptively select cluster cores and its numbers based on the distribution characteristics of data, and SPSA algorithm can successfully avoid getting stuck in local optima.

Table 1. Comparison of RMSE among AP-Fourier-SPSA-3D, AP-Fourier-GD-3D,

NNC-SVR-3D and KL-LS

model	AP-Fourier-SPSA-3D	AP-Fourier-GD-3D	NNC-SVR-3D	KL-LS
Training data	0.9229	0.9899	1.2862	2.4481
Test data	0.9027	0.9558	1.2214	2.2903

Table 1. Comparison of RMSE among AP-Fourier-SPSA-3D, AP-Fourier-GD-3D,

NNC-SVR-3D and KL-LS

model	AP-Fourier-SPSA-3D	AP-Fourier-GD-3D	NNC-SVR-3D	KL-LS
Training data	0.9229	0.9899	1.2862	2.4481
Test data	0.9027	0.9558	1.2214	2.2903

5. Conclusions

The study introduced a three-dimensional fuzzy modeling technique using SPSA learning for DPS, which widely exists in the manufacturing industry. On the three-dimensional fuzzy modeling framework, AP clustering was used to learn the preceding components, and adaptively determine the quantity of fuzzy rules. The three-dimensional fuzzy model used the Fourier space base function as the resulting components, while the SPSA method was employed to dynamically acquire the coefficients of the Fourier space base function. The suggested approach was implemented on a conventional distributed parameter system in the semiconductor manufacturing industry. The simulation experiment demonstrates the superiority of the three-dimensional fuzzy modeling method over the conventional space-time separation modeling method. Additionally, it shows that the Fourier space base function is more effective in revealing the spatial distribution characteristics of distributed parameter systems compared to the spatial function obtained through KL decomposition.