The nonlinear ARX model is used to represent a category of smooth nonlinear systems as below: [23]: where $y(a)$ represents output, $u(a)$ represents input, $\theta (a)$ represents white noise, $s_y$ and $s_u$ represent model orders. At an operating point ${\epsilon }_0$, $\eta (•)$ is expanded into the following Taylor polynomial

Then, substituting equation (2) into equation (1) yields the SD-ARX model as follows. where$\left\{{\rho }_{y,i}(\epsilon (a-1))|i=1, \mathrm{...} ,s_y\right\}$, $\left\{{\rho }_{u,j}(\epsilon (a-1))|j=1, \mathrm{...} ,s_u\right\}$ and ${\rho }_0(\epsilon (a-1))$ represent the regression coefficients that can be obtained by fitting using neural networks, including BRF neural networks [18,19,20,21,22,23] and wavelet neural networks [24]. $\epsilon (a-1)$ represents the operating state at time $a$, which may correspond to the system's input or/and output. When fixing $\epsilon (a-1)$, equation (3) is considered a locally linearized model, and when $\epsilon (a-1)$ varies with the system's operating point, equation (3) represents a model capable of globally describing the system's nonlinear properties. This model has local linear and global nonlinear features, making it advantageous for MPC design. In this article, we use LSTM or/and CNN to approximate the model’s regression coefficients and obtain a category of SD-ARX models that can effectively represent the system's nonlinear properties.

LSTM network serves as an upgraded variant of RNN, which can not only solve the problem that traditional RNN cannot handle the dependence of long sequence data, but also solve the issue that gradient exploding or vanishing is easy to occur with increasing training time and network layers of the neural networks. The LSTM network consists of multiple LSTM units, and its internal structure and chain structure are shown in Figure 1, with its core being the cell state $C$, represented by a horizontal line running through the entire cell at the top layer. It controls the addition or removal of information through gates. LSTM mainly consists of an input gate, a forget gate, and an output gate. The input gate controls which input information needs to be added to the current time step's storage unit, the forget gate controls which information needs to be forgotten in the storage unit of the previous time step, and the output gate controls which information in the current time step's storage unit needs to be output to the next time step. LSTM introduces a gate mechanism to control information transmission, remember content that needs to be memorized for a long time, and forget unimportant information, making it particularly effective in handling and forecasting critical events with long intervals and temporal delays in sequential data. The formula for an LSTM unit is as follows. where $l$ represents the time step of the input sequence, referred to as the number of rounds; $\odot$ represents matrix dot multiplication; $c_l$ and $h_l$ denotes cell state vector and hidden state vector, respectively; ${\sigma }_4(•)$ and ${\sigma }_5(•)$ denotes the nonlinear activation functions; $x_l$ represents input information; $f_l$ represents the forget gate’s output, which decides the amount of information from the preceding state $c_{l-1}$ needs to be forgotten. The value of the output ranges from 0 to 1, with a smaller value indicating more information to be forgotten, 0 means total forgetting, and 1 means total retention; $i_l$ indicates the input gate's output, which controls the current output state. Its output value is usually mapped to a range of 0 to 1 using the ${\sigma }_4(•)$ function, indicating the degree of activation of the output state, with 1 indicating full activation and 0 indicating no activation; $o_l$ represents the output gate's output, which decides the information to be output in this round; $W_{cx}$, $W_{ch}$, $W_{fx}$, $W_{fh}$, $W_{ix}$, $W_{ih}$, $W_{ox}$ and $W_{oh}$ denote weight coefficients, $b_f$, $b_i$, $b_o$ and $b_c$ denote offsets.

CNN has outstanding spatial feature extraction capabilities, mainly composed of the input layer, convolutional layer, pooling layer, fully connected layer, and output layer, as demonstrated in Figure 2. The input layer is employed for receiving raw data. The convolutional layer is utilized for extracting features from the input data. The pooling layer downsamples the convolution layer’s output to lower the data’s dimension. The fully connected layer turns the pooling layer’s output into the one-dimensional vector and outputs it through the output layer. Compared with the full connection mode of the traditional neural network, CNN uses a convolution kernel with appropriate size to extract the local features of this layer's input, that is, this CNN convolution layer's neurons are only connected to some upper layer's neurons. This local connectivity method reduces the neural network's parameter count and overfitting risk. At the same time, CNN employs weight sharing to enable convolutional kernels to capture identical features at varying places, resulting in reduced model training complexity and imparting the model with the characteristics of translation invariance. The convolution operation formula is below: where $\otimes$ represents convolution calculation; $x_{j_1}^{l_1}$ denotes the $j_1$-th feature map of the $l_1$-th convolution layer; $b_{j_1}^{l_1}$ represents offset; $W_{i_1{j}_1}^{l_1}$ represents convolution kernel matrix; ${\sigma }_1(•)$ represents the nonlinear activation functions; $g_1$ indicates the number of input feature maps.

LSTM neural network adopts directional circulation, which is better than feedforward neural network to deal with the problem of before and after correlation between data, and is conducive to the establishment of a nonlinear system timing model. Simultaneously, a gating mechanism is introduced in LSTM networks to control the flow and discarding of historical information and input characteristics, which solves the long-term dependence problem in simple RNN networks. In the actual modeling process, single-layer LSTM is usually unable to express complex temporal nonlinear characteristics, so a series approach is often adopted to stack multiple LSTM layers to deepen the network and bolster modeling capabilities. The LSTM-ARX model can be obtained by approximating the function coefficient of the model (3) with LSTM. It combines the benefits of LSTM for handling long sequence data with the SD-ARX model's nonlinear expression capabilities to fully represent the nonlinear system’s dynamic features. Figure 3 describes the LSTM-ARX model's architecture and its expression below where $y(a)$ represents output; $u(a)$ represents input; $\theta (a)$ represents white noise; $s_y$$s_u$,$d$are model's orders; $s_d$ denotes time delay; ${\rho }_0(\epsilon (a-1))$, ${\rho }_{y,i}(\epsilon (a-1))$ and ${\rho }_{u,j}(\epsilon (a-1))$ represent the model's function coefficient; n denotes the number of hidden layers; $h_l^r(a)$ and $c_l^r(a)$ represent the hidden layer state and cell state of the $r$-th hidden layer of time step $l$ at time $a$; ${\sigma }_4(•)$,${\sigma }_5(•)$ and $f(•)$ represent activation functions (e.g., Hard Sigmoid and Tanh) are employed to improve the model's capability in capturing the system's nonlinearity; $\{b_f^k,b_i^k,b_o^k,b_c^k|k=1,\mathrm{...},n\}$,$b_y$ and $b$ represent offsets; $\{W_{fx}^k,W_{fh}^k,W_{ix}^k,W_{ih}^k,W_{ox}^k,W_{oh}^k,W_{cx}^k,W_{ch}^k|k=1,\mathrm{...},n\}$, $W_{hy}$ and $W$ represent weight coefficients.

CNN has strong nonlinear feature extraction capabilities. The CNN-ARX model is derived by fitting the SD-ARX model's coefficients by CNN. This model combines the ability of CNN's local feature extraction and the SD-ARX model's expression. Compared to traditional neural networks, CNN has the characteristics of local connection, weight sharing, pooling operation, etc., which can further decrease the complexity and overfitting risks of the model, and enable the model to possess a certain level of robustness and fault tolerance. Figure 4 describes the CNN-ARX model's architecture and its expression below where $x^0(a)$ represents input; ${\sigma }_1(•)$ and $f_1(•)$ represent activation functions; ${\sigma }_2(•)$ and ${\sigma }_3(•)$ indicates pooling operation and flattening operation; ${\tilde{x}}^{n_1}(a)$ denotes the one-dimensional vector obtained by flattening the output feature maps ${\widehat{x}}_{j_1}^{n_1}(a)$; $\left\{W_{i_1{j}_1}^{l_1}\right|l_1=1,\mathrm{...},n_1\}$, $W_{xy}$ and $W_1$ represent weight coefficients; $\left\{b_{j_1}^{l_1}\right|l_1=1,\mathrm{...},n_1\}$，$b_x$ and $b_1$ represent the offsets.

The LSTM-CNN-ARX model is obtained by combining LSTM and CNN to approximate SD-ARX model coefficients, as shown in Figure 5. First, employ LSTM for learning the temporal characteristics of the input data, followed by the utilization of CNN to capture the spatial characteristics, and finally calculate the SD-ARX model’s correlation coefficients through a full connection layer to obtain the LSTM-CNN-ARX model. This model integrates the spatiotemporal characteristic extraction ability of LSTM-CNN with the expressive power of the SD-ARX model, this approach effectively captures the nonlinear system’s dynamic features. The LSTM-CNN-ARX model's mathematical expression is obtained from Formulas (6) and (7), as follows

The aforementioned deep learning-based SD-ARX models adopt data-driven techniques for modeling. First, select appropriate and effective input and output data as identification data for modeling. Second, based on the complexity of the identified system, select the model’s initial orders, the network's layer count, and the node count per layer. Third, choose activation functions (such as Sigmoid, Tanh, and Relu) and optimizer (such as Adam, and SGD), and configure other hyperparameters of the model. Fourth, train the model and evaluate its performance using a test set, computing the mean square error (MSE). Fifthly, adjust the orders while keeping the other structures of the model unchanged to train the model with new orders again, and select the orders of the model with the minimum MSE as the optimal orders. Next, keep the other structures of the model unchanged, adjust the network's layer count, and select the network's layer count with the smallest MSE as the optimal number of layers. Using the same method, determine the other model structures, i.e., the node count, activation function, and hyperparameters, and ultimately choose the model structure with the smallest MSE as the optimal model.

WTS is a representative dual-input dual-output process control experimental equipment with strong coupling, large time delay, and strong nonlinearity. The water inflows of water tanks 1 and 2 are controlled by electric control valves ECV1 and ECV2, the water outflows are regulated by proportional valves PV1 and PV2, and the water level heights are detected by the liquid level sensors LLS1 and LLS2. Control experiments were performed on WTS and the models (6)-(8) are designed below: where $Y(a)={[y_1(a),y_2(a)]}^{\mathrm{T}}$represents the liquid level heights of tank 1 and tank 2; $U(a)={[u_1(a),u_2(a)]}^{\mathrm{T}}$ represents the openings (0%~100%) of electric control valves ECV1 and ECV2; $\left\{G_{i,a-1}\right|i=1,\mathrm{...},s_y\}$, $\left\{H_{j,a-1}\right|j=s_d,\mathrm{...},s_u+s_d-1\}$ and ${\rho }_{0,a-1}$ are the weight coefficients, which can be obtained by models (6)-(8); $\theta (a)$ indicates the white noise signal. Among them, we have set the working point $\epsilon (a-1)={[Y{(a-1)}^{\mathrm{T}},\cdots ,Y{(a-d)}^{\mathrm{T}}]}^{\mathrm{T}}$, because the change in water level is an important factor leading to the strongly nonlinear dynamic characteristics of WTS.