1 Introduction

Microwave heating with the unique superiority of high efficiency has been widely applied in many fields over the past few decades, especially in recent industrial applications, high-power continuous microwave heating system (HPCMHS) has attracted great appeal and is successfully applied in heating or drying of wood [1], cement [2], ceramic [3] and food [4] etc. To improve the product quality and energy efficiency, increasing studies are focusing on the modeling of microwave heating process, since an accurate mathematical model can provide a platform for control and decision making. Mathematically, microwave heating process is modeled by a partial differential equations (PDE) with nonhomogeneous conditions. Such PDE is not only involved the characteristics of thermodynamics, but also coupled the conversion of microwave energy, which is computationally expensive [5]. Especially in HPCMHS, the system integrates multiple microwave sources by means of power synthesis to achieve large scale synchronous heating purpose. In such systems, each independent microwave source is capable of automatically adjusting its power, which contributes to more complex coupling of multi-physical fields and stronger interaction among system states. Consequently, such models normally consist of sets of PDEs with infinite-dimensional nature, although some numerical techniques can solve these PDEs, it can be computationally more demanding.

An attractive solution to deal with these drawbacks is to use system identification techniques that allow for solving such mathematically intractable and computationally complex problems. One such tool is artificial neural network (ANN). These tools when equipped with learning algorithms can modeling and control of any given systems [6,7,8,9,10,11,12]. Attracted by this, increasing researches have applied ANNs for modeling and prediction purposes in microwave heating or drying systems [13,14,15,16]. Most of these models have employed feed-forward neural network (FNN). In FNN the inputs feed forward through network layers to output and hence, only the forward connections are presented between neurons. However, in the case of nonlinear dynamic systems, especially for HPCMHS, the system output is a function of past history of inputs and outputs and the problem of identification gets prohibitive when FNN are used for this purpose. This is because they do not contain any memory feature in their structure and hence, are limited by static system mapping [17]. Responding to this deficiency, recurrent neural network (RNN) with self-feedback loops is able to capture the rich dynamics of nonlinear dynamic systems, as recurrent nature of neurons in RNN gives them ability to store information which can be used for later stage [18,19,20,21]. In [22], RNN can be classified into two categories: fully connected recurrent neural network (FCRNN) and partially connected recurrent neural network (PCRNN). The PCRNN has fewer parameters and requires shorter training time than the previous one. Since [23] has mentioned that, for system identification, the identifier should be chosen to have few parameters, as the more model parameters we use, the higher model randomness involved. With this goal in mind, a diagonal recurrent neural network (DRNN) is employed in this paper to identify HPCMHS, as DRNN has considerably fewer weights than the FCRNN and the network is simplified significantly [23].

The typical training process with RNNs determines all network parameters by minimizing sum-squared error between true output and simulated output from the network with supervised algorithms. The early popular algorithm is the gradient descend (GD) based backpropagation which usually converges to local minima due to its strong sensitivity on initial values of weights. To mitigate this problem, many evolutionary algorithms (EAs) can be used to search the optimal network parameters, such as genetic algorithm (GA) [24], particle swarm optimization (PSO) [25], differentiation evolution [26] and ant colony optimization [27]. The main deficiencies include: (1) low efficiency due to the high complexity of EAs; (2) the probability of finding the global optimum parameters will decrease exponentially as the dimensionality of the search space increases, which implies that such EAs are typically more difficult to find the optimal parameters when facing a large scale network structure. Further, filter-based methods are also applied to train RNNs for nonlinear system identification [28, 29]. But Kalman filter based learning algorithm is very sensitive to noises entered in system, which may yield inaccurate results. Besides, to guarantee faster convergence and stability of the system, Kumar et al. [30] have presented a DRNN based identification scheme with lyapunov stability based adaptive learning rate. Although accuracy can be improved compared with fixed learning rate, the conventional GD algorithm still meets the difficulty of trapping in local minima. To improve this, Kumar et al. [31] further proposed a lyapunov stability based gradient descend (LSGD) algorithm for DRNN-based adaptive control system design. The stability and robustness of the whole system can be improved using the lyapunov stability criterion, however, the use of new learning method for a DRNN identifier has not been tested. Additionally, these learning methods are only applicable to RNNs of fixed structure. Kumar et al. [30] has pointed out that the selection of optimal number of neurons in the hidden layer is still a very difficult task for DRNN.

Regarding to network construction, one simple way is to repeat these learning procedures with different network size until the optimal one is acquired based on some termination criterion [32]. This method is, however, computationally too expensive. More recent works on RNN training have focused on simultaneous structure and parameter learning using EAs. The objective function to be minimized could be the combination of training error and a measure of complexity of network structure [33,34,35]. Although such methods can simultaneously determine the structure and parameters of RNNs, the computational cost could be prohibitive for large problems because the increase in the dimension of the optimization problem is at least proportional to the square of the number of hidden nodes, which could directly deteriorate the efficiency and robustness of EAs. Another way of determining the model size is to use the subset selection algorithms, which has been widely used to construction of polynomial nonlinear autoregressive models with exogenous inputs (NARX) networks [36], radial-basis function (RBF) networks [37,38,39] and wavelet networks [40]. Such neural model structure is a linear combination of nonlinear functions, such as polynomials, Gaussian functions, and wavelets. Particularly noted that DRNN also have a standard network structure consisting of one recurrent layer with nonlinear functions (an activation function with tapped delay of its output) and a linear output, which can be formulated as the linear-in-the-parameters model. Therefore, candidate neuron pool can be first chosen a priory, a parsimonious network is then determined from these candidates using some subset selection methods, among which orthogonal least squares (OLS) is perhaps the most popular [36]. OLS approaches are derived from an orthogonal decomposition (OD) of the regression matrix. In OLS each regressor has to orthogonal to previously orthogonaliazed regressors and this operation is repeated until all regressors in the matrix are orthogonal to each other, which requires expensive computation. Although, a number of fast OLS algorithms have been proposed to improve the efficiency [41,42,43]. Their numerical stability can not be guaranteed.

In this paper, a novel DRNN construction is developed with two-stage method. In the first stage, fast recursive algorithm (FRA) is employed to construct DRNN. Unlike OLS that uses OD on the regression matrix, the more recent FRA uses a regression context where fast selection of model terms can be achieved on the original regression matrix. It has been shown that the FRA requires much less computational effort and is also numerically more stable than OLS [44]. Besides, in a traditional way such forward model selection algorithm is presented as an independent training process to replace the typical ANN training process, as the linear output parameters are updated by solving least squares problem [45]. The drawback is that such method does not fully use the information contained in a recurrent neuron, failing to exploit the full potential of the recurrent network structure. Thus, we propose using FRA method as just first stage within the context of an overarching methodology. And in our second stage, based on the reduced model size, the LSGD algorithm is used for tuning the parameters of DRNN to ensure the faster convergence and improved stability, which is not guaranteed by the conventional GD principle [31]. The combination of using FRA for an initial solution to set the number of recurrent neurons and their weights, while then using LSGD as complement to set optimal network parameters, is a novelty of this work. The proposed model has sufficient simplicity and accuracy that qualify it to be easily adopted in control system design for HPCMHS. The main contributions of this paper are summarized as follows.

  1. 1.

    Propose to use FRA at the first stage for DRNN structural optimization to replace the conventional OLS and EAs, which significantly reduce the computational complexity.

  2. 2.

    Unlike updating output weights by solving least squares, we propose to employ LSGD at the second stage to ensure system stability and mitigate the problem of stucking at local minima.

  3. 3.

    Integrate seamlessly the LSGD weight adaptation and FRA structural optimization into one DRNN-based identification scheme to handle the nonlinearity and complexity of HPCMHS.

The paper is organized as follows. Section 2 gives a brief introduction of HPCMHS with a data collecting experiment. Section 3 proposes a novel method for DRNN construction in terms of two key aspects: (1) forward model selection (2) backward error propagation. In Sect. 4 this method is then applied to modeling HPCMHS and the identification framework is presented. Section 5 gives simulation experiment, the performances of five identifiers are tested and compared. Finally, in Sect. 6, the contributions are summarized, broader implications are discussed, and potential work is noted.

2 High-Power Continuous Microwave Heating System

The experimental HPCMHS is shown in Fig. 1a. The system is mainly constructed of four parts, a power generation part, a control part constructed by a host computer and PLC, microwave transmission part, and multimode cavity integrated with a conveyor belt. Five microwave sources constitute the power generation part to achieve the large-scale heating purpose. The position of each microwave source is shown in Fig. 1b, in which two pairs of microwave source are connected in series, while a single one is employed in the last cavity. The magnetron integrated in each source enables generate microwave power ranging from 0 to 3000 W, which is capable of continuous adjustment. The microwave power is then directed into the cavity by transmission part. In the microwave heating process, materials distributed evenly on the conveyor belt is continuously transported from the inlet to outlet. An optical fiber temperature sensor is employed at the outlet of the cavity to measure the temperature of exported material. PLC is used to establish the required connection between device and computer, the power setting can be send to PLC to lead the power generation part by Modbus protocol, while the microwave powers, temperature measurement and conveyor speed are real time recorded until the heating process is finished.

In this context, a data collecting experiment is conducted to obtain the required input/output data set by a real-time microwave rice drying experiment. The objective of the experiment is to acquire the dried rice in the outlet through appropriate regulation of microwave powers and conveyor speed. While the heating process began, continuous rice distributed evenly on the conveyor belt is continuously transported through the cavity. During the process, microwaves penetrate the rice sample, heating the kernel water until it defused to surface. The humid air is then drawn out of the cavity by a suction system. Temperature variation of exported rice is detected by optical sensor. At the end of heating process, a vast volume of data including microwave powers, rice temperatures and conveyor speed are saved in host computer, thus a set of input/output data sets ready to be used in system identification to develop the desired model are obtained.

Fig. 1
figure 1

Schematic of experimental high-power continuous microwave heating system

Full size image

3 Two-Stage Algorithm for DRNN Construction

3.1 DRNN Structure

Figure 2 shows the structure of DRNN. Input weight vector connected external inputs to hidden layer is denoted as \({W^I}\). Neuron presented in the hidden layer is known as recurrent neuron that has an internal feedback loop to itself. These feedback loops are called diagonal weights which provide a weighted unit delay output of recurrent neuron. This is the reason that DRNN can store information for later use and is better at dealing with time-varying input and output. The diagonal weights are represented as \({W^D} =\{ w_1^D,w_2^D\ldots w_m^D\} \), where m denotes the number of recurrent neurons in the hidden layer. The output weight is denoted by \({W^O} = \{ w_1^O,w_2^O\ldots w_m^O\} \), while the input vector is \(X = \{ {x_1}(t),{x_2}(t)\ldots {x_n}(t)\} \), where n denotes the input dimension. Then the mathematical model for DRNN is shown below

$$\begin{aligned} y(t)= & {} \sum \limits _m {W_m^O{H_m}(t)} \end{aligned}$$
(1)
$$\begin{aligned} {H_m}(t)= & {} \varphi ({S_m}(t)) \end{aligned}$$
(2)
$$\begin{aligned} {S_m}(t)= & {} W_m^D{H_m}(t - 1) + \sum \limits _n {W_{nm}^I{x_n}(t)} \end{aligned}$$
(3)

where for each time step t, \({W^I}\) denotes the input weights, \({S_m}(t)\) is the sum of inputs of mth recurrent neuron, \({H_m}(t)\) is the output of mth recurrent neuron, and y(t) is the output of DRNN. It should be noted that activation function \(\varphi ( \cdot )\) chosen in the recurrent neuron is usually nonlinear tangent hyperbolic function which enables the output have both positive and negative value, while the activation in the output neuron is considered to be linear so that there is no restriction on its value.

Fig. 2
figure 2

Structure of DRNN

Full size image

3.2 Stage 1-Forward Model Selection by FRA

A linear-in-the-parameters structure for DRNN can be formulated as a linear combination of nonlinear functions, which is given by

$$\begin{aligned} y(t) = \sum \limits _{i = 1}^m {{\varphi _i}\left( w_i^D{H_i}(t - 1) + \sum \limits _{j = 1}^n {w_{ji}^I{x_j}(t)} \right) } W_m^O + \varepsilon (t) \end{aligned}$$
(4)

where \({\varphi _i}( \cdot ),\; i = 1,\ldots ,m\) are all candidate recurrent neurons, and \(\varepsilon (t)\) is the model residual sequence.

Remark 1

the total number (m) of candidate neurons can initially be significantly large, which may cause over-fitting. Therefore, it is important to find a parsimonious model with a much smaller number, say \(k\;(k < m)\) of terms for DRNN construction.

Suppose N data samples \(\{ x(t),y(t)\} _{t = 1}^N\) are used for model identification, (4) can then be formulated as

$$\begin{aligned} y = \varPhi \varTheta + \varXi \end{aligned}$$
(5)

where

\(\begin{array}{l} \varTheta = W_m^O = {[w_1^O,w_2^O,\ldots ,w_m^O]^T} \in {R^m}\\ \varPhi = [{\varphi _1},{\varphi _2},\ldots ,{\varphi _m}] \in {R^{N \times m}}\\ {\varphi _i} = [{\varphi _i}(w_i^D{H_i}(0) + \sum \limits _{j = 1}^n {w_{ji}^I{x_j}(1)} ),{\varphi _i}(w_i^D{H_i}(1) + \sum \limits _{j = 1}^n {w_{ji}^I{x_j}(2)} )),\ldots ,{\varphi _i}(w_i^D{H_i}(N)\\ \qquad \quad + \sum \limits _{j = 1}^n {w_{ji}^I{x_j}(N)} ){]^T} \in {R^N}\quad i = 1,2,\ldots ,m\\ \varXi = {[\varepsilon (1),\varepsilon (2),\ldots ,\varepsilon (N)]^T} \in {R^N} \end{array}\)

In order to select the recurrent neurons and to build a compact DRNN model, the net contribution (the contribution of recurrent neuron) \(\delta {E_k}\) of each term chosen for DRNN needs to be first computed explicitly. In OLS, this is done by an OD to the regression matrix \(\varPhi \), which requires heavy computation. In this paper, FRA is employed to compute this contribution by solving the least-squares problem recursively.

First, the cost function \({E_k}\) of first k columns is given as [44]

$$\begin{aligned} \left. \begin{array}{l} {E_k} = {y^T}y - {\hat{\varTheta }} _k^T\varPhi _k^Ty\\ {{{\hat{\varTheta }} }_k} = {(\varPhi _k^T{\varPhi _k})^{ - 1}}\varPhi _k^Ty \end{array} \right\} ,\ \ k = 1,\ldots ,m. \end{aligned}$$
(6)

where \({\varPhi _k} \in {R^{N \times k}}\) contains the first k columns of the full regression \(\varPhi \) in (5). \({{\hat{\varTheta }} _k} \in {R^k}\) is the estimate of the first k parameters in \(\varTheta \).

To introduce the forward model selection procedure, a matrix series is defined

$$\begin{aligned} {R_k} = \left\{ \begin{array}{l@{\quad }l} I - {\varPhi _k}{(\varPhi _k^T{\varPhi _k})^{ - 1}}\varPhi _k^T, &{} 0 <k \le m,\\ I, &{} k = 0. \end{array} \right. \end{aligned}$$
(7)

The matrix \({R_k}\) is a residue matrix, which holds the following properties [44]

$$\begin{aligned} {R_{k + 1}}= & {} {R_k} - \frac{{{R_k}{\varphi _{k + 1}}\varphi _{k + 1}^TR_k^T}}{{\varphi _{k + 1}^T{R_k}{\varphi _{k + 1}}}},\ \ k = 0,1,\ldots ,(m - 1). \end{aligned}$$
(8)
$$\begin{aligned} R_k^T= & {} {R_k}; ({R_k})^{2} ={R_k},\quad k = 0,1,\ldots ,m. \end{aligned}$$
(9)
$$\begin{aligned} {R_i}{R_j}= & {} {R_j}{R_i},\ \ i \ge j;\ i,j = 0,1,\ldots ,m. \end{aligned}$$
(10)
$$\begin{aligned} {R_k}\varphi= & {} \left\{ \begin{array}{l} 0, 0 < k \le m,\quad rank([{\varPhi _k},\varphi ]) = k\\ {\varphi ^{(k)}} \ne 0,\quad rank([{\varPhi _k},\varphi ]) = k + 1 \end{array} \right. \end{aligned}$$
(11)

where \({\varphi ^{(k)}} = {R_k}\varphi \).

Now, from (8)–(11) and the definition of \({R_k}\) in (7), cost function \({E_k}\) can be expressed as

$$\begin{aligned} {E_k} = {y^T}y - {\hat{\varTheta }} _k^T\varPhi _k^Ty = {y^T}{R_k}y \end{aligned}$$
(12)

Then, the net contribution of \({\varphi _{k + 1}}\) to the cost function is given by

$$\begin{aligned} \delta {E_{k + 1}} = {E_{k + 1}} - {E_k} = - \frac{{{y^T}{R_k}{\varphi _{k + 1}}\varphi _{k + 1}^TR_k^Ty}}{{\varphi _{k + 1}^T{R_k}{\varphi _{k + 1}}}} \end{aligned}$$
(13)

Further, defining

$$\begin{aligned} \varphi _i^{(k)} = {R_k}{\varphi _i},\ \ i = 1,\ldots ,m,\ k = 1,\ldots ,m. \end{aligned}$$
(14)

(13) then becomes

$$\begin{aligned} \delta {E_{k + 1}} = - \frac{{{{({y^T}\varphi _{k + 1}^{(k)})}^2}}}{{{{(\varphi _{k + 1}^{(k)})}^T}\varphi _{k + 1}^{(k)}}},\ \ k = 1,\ldots ,m - 1. \end{aligned}$$
(15)

To further simplify (15), two new notations are needed

$$\begin{aligned} \left. \begin{array}{l} {a_{k,i}} = {(\varphi _k^{(k - 1)})^T}\varphi _i^{(k - 1)}\\ {a_{k,y}} = {(\varphi _k^{(k - 1)})^T}y,\\ i = 1,\ldots ,n,k = 1,\ldots ,n \end{array} \right\} . \end{aligned}$$
(16)

Finally, according to [44], the net contribution of \({\varphi _{k + 1}},\ k = 1,\ldots ,m - 1\) to the cost function can be explicitly expressed as

$$\begin{aligned} \delta {E_{k + 1}} = - \frac{{{{\left( {y^T}{\varphi _{k + 1}} - \sum \limits _{j = 1}^k {({{{a_{j,y}}{a_{j,k + 1}}} \Big / {{a_{j,j}}}})} \right) ^2}}}}{{{{({\varphi _{k + 1}})}^T}{\varphi _{k + 1}} - \sum \limits _{j = 1}^k {({{a_{j,k + 1}^2} \Big / {{a_{j,j}}}})} }} \end{aligned}$$
(17)

(17) expresses the net contribution of a selected regressor term to the cost function. In general, the error reduction ratio is employed to represent the proportion of the dependent neuron variance, which can be defined by [36]

$$\begin{aligned} {[{e_{rr}}]_{k + 1}} = {{ - \delta {E_{k + 1}}} \Big / {{y^T}y}} \end{aligned}$$
(18)

This can be used to calculate the contribution of each recurrent neuron in candidates pool. The detail neuron selection procedure is given as follows:

Step 1: To add the first recurrent neuron in DRNN. For \(i = 1,\ldots ,m\), where m denotes the number of candidate recurrent neurons that are generated randomly by DRNN. The error reduction ratio in the first step can be defined as

$$\begin{aligned} {[{e_{rr}}]_{1i}} = \frac{{{{({y^T}{\varphi _{1i}})}^2}}}{{\varphi _{1i}^T{\varphi _{1i}}{y^T}y}} \end{aligned}$$
(19)

where \({\varphi _{1i}}\) represents candidate recurrent neuron in the first step. Assume that \([{e_{rr}}]_1^j = max\{ {[{e_{rr}}]_{1i}},i = 1,\ldots ,m\} \). j denotes the sequence number of candidate neuron. Then the corresponding neuron \({\varphi _1} = {\varphi _{1j}}\) is selected as the first neuron chosen in the DRNN model, and \(w_{1j}^D\) and \(w_{1j}^I\) are selected as the initial diagonal weight and input weight of the first selected recurrent neuron.

Step 2: To add the second neuron in DRNN. The best neuron is selected from the rest of candidate neurons (the first selected neuron has been removed from the candidate pool). For \(i = 1,\ldots ,m - 1\), compute

$$\begin{aligned} {[{e_{rr}}]_{2i}}= & {} \frac{{{{({y^T}{\varphi _{2i}} - {{{a_{1,y}}{a_{1,2i}}} \Big / {{a_{1,1}}}})}^2}}}{{(\varphi _{2i}^T{\varphi _{2i}} - {{a_{1,2i}^2} \Big / {{a_{1,1}})({y^T}y)}}}}\nonumber \\= & {} \frac{{{{({y^T}{\varphi _{2i}} - {{\varphi _1^Ty\varphi _1^T{\varphi _{2i}}} \Big / {\varphi _1^T{\varphi _1}}})}^2}}}{{(\varphi _{2i}^T{\varphi _{2i}} - {{{{(\varphi _1^T{\varphi _{2i}})}^2}} \Big / {\varphi _1^T{\varphi _1})({y^T}y)}}}} \end{aligned}$$
(20)

Again, assume that \([{e_{rr}}]_2^j = max\{ {[{e_{rr}}]_{2i}}, i = 1,\ldots ,m - 1\} \). Then \({\varphi _2} = {\varphi _{2j}}\) is selected as the second neuron, and \(w_{2j}^D\) and \(w_{2j}^D\) are selected as the corresponding weights of the second selected neuron. The calculation of error reduction ratio at kth stage is generalized and simplified as follow

$$\begin{aligned}{}[{e_{rr}}]_{ki} = = \frac{{{{\left( {y^T}{\varphi _{ki}} - \sum \nolimits _{j = 1}^{k - 1} {{{\varphi _j^T{R_{j - 1}}y\varphi _j^T{R_{j - 1}}{\varphi _{ki}}} \Big / {\varphi _j^T{R_{j - 1}}{\varphi _j}}}} \right) }^2}}}{\left( \varphi _{ki}^T{\varphi _{ki}} - \sum \nolimits _{j = 1}^{k - 1} (\varphi _j^T{R_{j - 1}}{\varphi _{ki}})^{2} \Big / \varphi _j^T{R_{j - 1}}{\varphi _j}\right) ({y^T}y)} \end{aligned}$$
(21)

The selection procedure is continue until the kth stage when

$$\begin{aligned} 1 - \sum \limits _{i = 1}^k {{{[{e_{rr}}]}_i} < \delta } \end{aligned}$$
(22)

where \(\delta \) is a desired tolerance for error reduction ratio (set to small value, e.g., 0.01).

Remark 2

The FRA procedure works iteratively to calculate \({[{e_{rr}}]_i}\), \({\varphi _i}\), and then add the new recurrent neuron after each iteration until the termination criteria is satisfied. This offers several advantages. Firstly, a parsimonious compact DRNN structure can be obtained as the size is determined automatically. Secondly, since FRA always selects the neuron that maximizes the decrease to the variance of the desired output among all possible choices, a smaller number of recurrent neurons will be constructed, or fewer parameters of the model will be adjusted, which considerably reduce the randomness involved in an identifier [23]. Thirdly, since the candidate neuron pool is generated by propagating the input to the recurrent layer with some random weights, the initial input and diagonal weights can be simultaneously determined as the suitable neurons are selected, which significantly improve the convergence speed at the second training stage and mitigate the problem of trapping in local minima.

To this end, the optimal DRNN structure and the initial input and diagonal weights have been determined through FRA. In the following, the overall training based on LSGD will be introduced in detail.

3.3 Stage 2-Backward Error Propagation by LSGD

Stage 2 illustrates the parameter learning procedure. To lay the foundation of weight learning algorithm, the generalized form of weight update equation is given as

$$\begin{aligned} {W^g}(t + 1) = {W^g}(t) - \eta \varDelta {W^g}(t) \end{aligned}$$
(23)

where \({W^g}\) denotes the generalized weight of DRNN (including input, diagonal and output weights), \(\varDelta {W^g}\) denotes the required adjustment to be made and \(\eta \) represents learning rate. For conventional GD method, the purpose is to reduce the mean square error which is defined as

$$\begin{aligned} E(t) = \frac{1}{2}e{(t)^2} \end{aligned}$$
(24)

where e(t) is the error between true output and simulated output from the network. Then \(\varDelta {W^g}\) according to the conventional GD principle is given as

$$\begin{aligned} \varDelta {W^g}(t) = \frac{{\partial E(t)}}{{\partial e(t)}}\frac{{\partial e(t)}}{{\partial {W^g}(t)}} \end{aligned}$$
(25)

The details of update rules for input, diagonal and output weights can be seen in [30]. As (25) may get into the problem of local minima, a Lyapunov stability criterion is used to derive the weight adjustment rule in this paper and a Lyapunov function is chosen as

$$\begin{aligned} {L_c}(t) = \frac{1}{2}[e{(t)^2} + {W^g}{(t)^2}] \end{aligned}$$
(26)

Since this defined function is positive, the second condition for asymptotic stability must be satisfied in turns to give us the condition on \(\varDelta {W^g}\). The discrete-time error of Lyapunov function is written as

$$\begin{aligned} \varDelta {L_c}(t) = \frac{1}{2}\left\{ {[e{{(t + 1)}^2} + {W^g}{{(t + 1)}^2}] - [e{{(t)}^2} + {W^g}{{(t)}^2}]} \right\} \end{aligned}$$
(27)

To further simplify it

$$\begin{aligned} \varDelta {L_c}(t) = \frac{1}{2}\left\{ {[\varDelta e(t) + 2e(t)]\varDelta e(t) + [\varDelta {W^g}(t) + 2{W^g}(t)]\varDelta {W^g}(t)} \right\} \end{aligned}$$
(28)

where \(\varDelta e(t) = e(t + 1) - e(t)\) and \(\varDelta {W^g}(t) = {W^g}(t + 1) - {W^g}(t)\). To rearrangement of terms in (28), we can get

$$\begin{aligned} \varDelta {L_c}(t)= & {} \frac{1}{2}\varDelta {W^g}{(t)^2}\left\{ {1 + {{\left[ \frac{{\varDelta e(t)}}{{\varDelta {W^g}(t)}}\right] }^2}} \right\} \nonumber \\&+ \varDelta {W^g}(t)\left\{ {{W^g}(t) + e(t)\left[ \frac{{\varDelta e(t)}}{{\varDelta {W^g}(t)}}\right] } \right\} \end{aligned}$$
(29)

In (29), we usually set \(\frac{{\varDelta e(t)}}{{\varDelta {W^g}(t)}} = \frac{{\partial e(t)}}{{\partial {W^g}(t)}}\), thus

$$\begin{aligned} \varDelta {L_c}(t)= & {} \frac{1}{2}\varDelta {W^g}{(t)^2}\left\{ {1 + {{\left[ \frac{{\partial e(t)}}{{\partial {W^g}(t)}}\right] }^2}} \right\} \nonumber \\&+ \varDelta {W^g}(t)\left\{ {{W^g}(t) + e(t)\left[ \frac{{\partial e(t)}}{{\partial {W^g}(t)}}\right] } \right\} \end{aligned}$$
(30)

Theorem 1

For Lyapunov stability criterion, a given system is said to be asymptotic stable if \({L_c}(t) > 0,\ \varDelta {L_c}(t) \le 0\). This condition is satisfied only if

$$\begin{aligned} \varDelta {W^g}(t) = - \frac{{\left\{ {{W^g}(t) + e(t)\left[ \frac{{\partial e(t)}}{{\partial {W^g}(t)}}\right] } \right\} }}{{1 + {{\left[ \frac{{\partial e(t)}}{{\partial {W^g}(t)}}\right] }^2}}} \end{aligned}$$
(31)

Proof

From (30), let

$$\begin{aligned} \varDelta {L_c}(t)= & {} \frac{1}{2}\varDelta {W^g}{(t)^2}\left\{ {1 + {{\left[ \frac{{\partial e(t)}}{{\partial {W^g}(t)}}\right] }^2}} \right\} \nonumber \\&+ \varDelta {W^g}(t)\left\{ {{W^g}(t) + e(t)\left[ \frac{{\partial e(t)}}{{\partial {W^g}(t)}}\right] } \right\} = - \frac{1}{2}l \end{aligned}$$
(32)

where \(l \ge 0\) in order to satisfy \(\varDelta {L_c}(t) \le 0\). So (32) becomes

$$\begin{aligned} \varDelta {L_c}(t)= & {} \frac{1}{2}\varDelta {W^g}{(t)^2}\left\{ {1 + {{\left[ \frac{{\partial e(t)}}{{\partial {W^g}(t)}}\right] }^2}} \right\} \nonumber \\&+ \varDelta {W^g}(t)\left\{ {{W^g}(t) + e(t)\left[ \frac{{\partial e(t)}}{{\partial {W^g}(t)}}\right] } \right\} + \frac{1}{2}l = 0 \end{aligned}$$
(33)

Considering \(\varDelta {L_c}(t)\) as a quadratic equation form \(\varDelta {L_c}(t) = a\varDelta {W^g}{(t)^2} + b\varDelta {W^g}(t) + c = 0\), where \(a = \frac{1}{2}\left\{ {1 + {{[\frac{{\partial e(t)}}{{\partial {W^g}(t)}}]}^2}} \right\} \), \(b = \left\{ {{W^g}(t) + e(t)[\frac{{\partial e(t)}}{{\partial {W^g}(t)}}]} \right\} \) and \(c = \frac{1}{2}l\). To have a single unique solution of \(\varDelta {L_c}(t)\), \(\sqrt{{b^2} - 4ac} \) must be equal to zero. Putting values into \(\sqrt{{b^2} - 4ac} = 0\), we get

$$\begin{aligned} \sqrt{{{\left\{ {{W^g}(t) + e(t)\left[ \frac{{\partial e(t)}}{{\partial {W^g}(t)}}\right] } \right\} }^2} - \left\{ {1 + {{\left[ \frac{{\partial e(t)}}{{\partial {W^g}(t)}}\right] }^2}} \right\} l} = 0 \end{aligned}$$
(34)

Simplify (34) and l comes out to be

$$\begin{aligned} l = \frac{{\left\{ {{W^g}(t) + e(t)\left[ \frac{{\partial e(t)}}{{\partial {W^g}(t)}}\right] } \right\} }}{{1 + {{\left[ \frac{{\partial e(t)}}{{\partial {W^g}(t)}}\right] }^2}}} \end{aligned}$$
(35)

As mentioned before \(l \ge 0\), which means

$$\begin{aligned} \frac{{\left\{ {{W^g}(t) + e(t)\left[ \frac{{\partial e(t)}}{{\partial {W^g}(t)}}\right] } \right\} }}{{1 + {{\left[ \frac{{\partial e(t)}}{{\partial {W^g}(t)}}\right] }^2}}} \ge 0 \end{aligned}$$
(36)

Thus, the unique root of \(\varDelta {L_c}(t)\) will be given as

$$\begin{aligned} \varDelta {W^g}(t) = - \frac{b}{{2a}} = - \frac{{\left\{ {{W^g}(t) + e(t)\left[ \frac{{\partial e(t)}}{{\partial {W^g}(t)}}\right] } \right\} }}{{1 + {{\left[ \frac{{\partial e(t)}}{{\partial {W^g}(t)}}\right] }^2}}} \end{aligned}$$
(37)

\(\square \)

Using the Lyapunov stability based weight adjustment, the weight update equation in (23) can be rewritten as

$$\begin{aligned} {W^g}(t + 1) = {W^g}(t) + \eta \frac{{\left\{ {{W^g}(t) + e(t)\left[ \frac{{\partial e(t)}}{{\partial {W^g}(t)}}\right] } \right\} }}{{1 + {{\left[ \frac{{\partial e(t)}}{{\partial {W^g}(t)}}\right] }^2}}} \end{aligned}$$
(38)

(38) gives LSGD weight update rule for DRNN. Because the design of weight adjustment is based on the Lyapunov stability criterion, the stability of the whole system can be guaranteed and the problem of local minima is mitigated.

4 DRNN Based HPMHS Identification

System identification is based on applying a suitable input signal to the system to be identified in order to excite its dynamics and consequently to observe its response to the applied input signal. The result of this process is a set of input/output data that can be used to develop the desired system model. In this context, a data collecting experiment has been introduced in Sect. 2. In this data set, the input signals to HPCMHS are five microwave powers and conveyor speed. Turning up the output power can enhance the electromagnetic field, while turn down the conveyor speed can prolong the microwave radiation time and thus raising the temperature, vice versa. Table 1 shows the process data acquired in the experiment.

Table 1 Experimental data during microwave rice drying process
Full size table

In [30], author has proposed a DRNN-based identification scheme for nonlinear dynamic systems. It only requires the plants present value \({y_p}(t)\) and one past value \({y_p}(t - 1)\) along with present value of externally applied input u(t) in order to compute its next value \({y_{DRNN}}(t + 1)\). This identification scheme can significantly reduce the number of parameters to be tuned in a DRNN identifier compared with FNN identifiers, thus shortening training time. Its mathematical equation is as follow:

$$\begin{aligned} {y_p}(t + 1) = {F_{DRNN}}[{y_p}(t),{y_p}(t - 1), u(t)] \end{aligned}$$
(39)

In HPCMHS, u(t) consists of five microwave powers and conveyor speed, while \({y_p}(t)\) denotes the exported rice temperature.

Before training, the input normalization is necessary. Although normalization step is not novel in and of itself, including it within the identification scheme is crucial. As a preprocessing step, normalization is performed to decrease the relative variance of each element in the input vector x. In addition, normalization is helpful to facilitate the convergence speed of training process. Equation (40) shows the normalization equation applied to the ith component of x, N denotes the number of training set. The consequent normalized value is between 0 and 1.

$$\begin{aligned} \left. \begin{array}{l} {{{\bar{x}}}_i} = \frac{{{x_i} - {x_{{i_{\min }}}}}}{{{x_{{i_{\max }}}} - {x_{{i_{\min }}}}}} \in (0,1)\\ {x_{{i_{\min }}}} = \min [x_i^1,x_i^2,\ldots x_i^N]\\ {x_{{i_{\max }}}} = \max [x_i^1,x_i^2,\ldots x_i^N] \end{array} \right\} . \end{aligned}$$
(40)

The DRNN based identification involves a two-stage training process for determining network structure and necessary parameters, each of stage is designed to improve accuracy and thus enable effective for representing HPCMHS. The identification scheme is presented in Fig. 3. First, the network training inputs x are normalized by (40). Then, a candidate neuron pool is generated randomly and the best recurrent neuron is selected by FRA at each iteration. Finally, LSGD method is employed to train the fixed size network model until the desired error is reached.

Fig. 3
figure 3

Flow chart for the new training process for DRNN

Full size image

5 Simulation Study

In order to validate the effectiveness of the newly developed identifier, this section compares the new approach with other common ANNs. The new network, henceforth referred to as FRA-LSGD-DRNN, is evaluated in terms of its computational complexity and accuracy. The computation involved is dominated by the forward model selection process and the final total count of network parameters to be tuned. The accuracy is determined by testing error which is calculated based on two common methods (1) mean square error (MSE), which measures the average magnitude of the error over the test sample, and (2) mean absolute error (MAE), which measures the average magnitude of the errors in a set of forecast predictions. The proposed identifier is first tested with an experimental regression problem with different inputs. Then, its performance is evaluated in HPCMHS identification. The work of training networks was completed on a Windows 7 computer with AMD A8-4500M processor and 4GB RAM.

In the two examples, FRA-LSGD-DRNN is compared against four ANNs, namely OLS-LSGD-DRNN, LSGD-DRNN, DRNN and FNN. In the OLS based DRNN, the forward model selection process is done by a matrix decomposition and orthogonal transformation. Its computational complexity will be compared with FRA based DRNN. For LSGD-DRNN, DRNN and FNN based identifiers, the optimal number of hidden neurons cannot be determined automatically. In order to find the number of hidden neurons, several methods have been developed. Among them, grid-search is probably the most commonly used and liable searching methods [32]. In grid-search method, all possible values are checked and the one that gives the best performance is picked up. The grid based search method can offer a suitable solution for hidden neuron selection but with high computational cost. In terms of parameter learning, the LSGD method is compared with conventional GD principle (DRNN and FNN identifier) in terms of identification accuracy and system stability.

Learning rate is also a crucial factor in ANN training. It decides the speed at which the parameters of identifier will adjust during the training. If \(\eta \) is chosen closer to 1 then it may leads instability and trap into local minimal. On the other hand a small value of \(\eta \) will have a very slow learning process as little improvement will be made in parameter values from one iteration to the next. Generally, a relatively large learning rate will be adopted at the start of training to shorten the convergent time, while a small value of \(\eta \) is chosen at the end of training to prevent trapping into local minimal. In our study, \(\eta \) is chosen to decrease along the time steps.

$$\begin{aligned} \eta = {\eta _0}{e^{ - \alpha t}} \end{aligned}$$
(41)

where \({\eta _0}\) denotes the initial learning rate, \(\alpha \) denotes the decrease rate and t refers to simulation time.

5.1 Experimental Regression Problem

5.1.1 Validation with Fixed System Input

Consider a nonlinear dynamical system [30], shown as follows

$$\begin{aligned} \left. \begin{array}{l} {y_p}(t + 1) = \frac{{{y_p}(t){y_p}(t - 1)}}{{1+{y_p}{{(t)}^2}+ {y_p}{{(t- 1)}^2}}}+u(t)\\ u(t) = \sin \left( \frac{{2\pi t}}{{25}}\right) \end{array} \right\} . \end{aligned}$$
(42)

The initial parameters set for five identifiers are shown in Table 2.

Table 2 Parameter setting for various models
Full size table
Table 3 Computational complexity of FRA and OLS based DRNN
Full size table

Given a desired tolerance, both FRA-LSGD-DRNN and OLS-LSGD-DRNN can determine the optimal recurrent neurons automatically. Suppose m candidate neurons are initially generated randomly, from which only k terms are to be selected. The forward model selection is to select the k most significant recurrent neurons when N data samples are used for identification. Then, the total computational complexity for OLS is as follow

$$\begin{aligned} {C_{OLS}} = 2mN({k^2} + 2k - 1) - 4kN({k^2} - 1)/3 + k(k - 1)/2. \end{aligned}$$
(43)

Similarly, the computational complexity for FRA is

$$\begin{aligned} {C_{FRA}} = [(2m - k)(k + 1) + 2m]N + km(k + 4) - 6m - k(2{k^2} + 3k - 17)/3 - 2.\nonumber \\ \end{aligned}$$
(44)

Table 3 gives the exact computational complexity when both algorithms are applied to select DRNN recurrent neurons under different circumstances. It is clearly to seen that only 5 most significant neurons are chosen when tolerance is set to \(\delta = 0.01\), thus a parsimonious compact DRNN structure can be obtained. Along with the increase of data samples N and initial candidate neurons m, the computational complexity of OLS based DRNN is approximately 5 times higher than FRA based DRNN, the computational effort of later one is significantly reduced. For LSGD-DRNN, the suitable network structure can be found by the grid-search method. Figure 4 shows the grid-searching process, it is obvious that more recurrent neurons lead to a smaller testing error at the beginning, but the generalization performance becomes poor for large number of hidden neurons. The best generalization ability was found when 5 neurons were used, after that, the network became over fitted. Although grid-search method yields to the same result with forward model selection algorithm, requires heavy computational cost, which again demonstrates the effectiveness and efficiency of the FRA based structure learning method.

Fig. 4
figure 4

Test MSE versus the growth of hidden neurons number for LSGD-DRNN

Full size image

Forward model selection process is followed by the parameter learning stage. The testing error of five identifiers are shown in Fig. 5. It is clear that the LSGD based DRNNs outperform other identifiers with relatively high accuracy. The conventional GD-based ANN experience difficulties in predicting this example, as evidenced by the large error. Although OLS-LSGD-DRNN and LSGD-DRNN have the closet performance with the proposed model, requires demanding computation when determining the structure. The various details for five identifiers averaged over 3 times are shown in Table 4. It can be seen that the proposed network structure offers fewer number of parameters to be tuned as compared to conventional DRNN and FNN model used in [30] with the same example, which directly contributes to the significantly reduced training time in the second training stage. Also, the average MSE and MAE obtained by proposed identifier is least among the other identifiers, which validates the outperformance achieved by FRA-LSGD-DRNN. This improved performance can also be demonstrated by the speeder error reduction that the proposed identifier compared with other identifiers, which is shown in Fig. 6. This trend is especially prominent when LSGD-based optimization is used, while the conventional GD-based DRNN and FNN even fail to become convergent within 30 training epochs.

Fig. 5
figure 5

Testing result of various models with fixed system input

Full size image
Table 4 Performance comparison of various identifiers
Full size table
Fig. 6
figure 6

Test MSE along the epoch of various identifiers

Full size image

5.1.2 Validation with External System Input

To further demonstrate the stability and robustness of the LSGD-based training method, a new external input is given to the system to test the resulting identifiers. This new external input is given by

$$\begin{aligned} u(t) = \left\{ \begin{array}{ll} \sin \left( \frac{{2\pi t}}{{250}}\right) &{} t \le 100,\\ 0.8 \times \sin \left( \frac{{2\pi t}}{{250}}\right) + 0.2 \times \sin \left( \frac{{2\pi t}}{{25}}\right) &{} 100 < t \le 500. \end{array} \right. \end{aligned}$$
(45)

The corresponding test results in Fig. 7 show superior performance for the LSGD-based identifiers. The averaged MSE and MAE of five identifiers are shown in Table 5, which again indicates that the proposed method have better generalization ability of approximating the unknown dynamics as compared to conventional GD-based training algorithm.

Fig. 7
figure 7

Testing result of various models with external system input

Full size image

5.2 HPCMHS Identification

The proposed FRA-LSGD-DRNN is then employed to identify HPCMHS, and its computational complexity and accuracy is evaluated and compared with other identifiers. In this context, the new network helps provide a computationally efficient tool for temperature prediction in microwave heating process. Especially in the field of control system design, it is critical to establish a simple and control-oriented model with high fidelity and low complexity to capture the rich dynamics of HPCMHS. ANNs provide a powerful tool for creating fast running models.

In this case, 3000 data sets are obtained during the rich drying experiment, where 2100 data sets are used for training DRNN with proposed two-stage method, while the rest of data is used to evaluate the generalization capacity of the resulting model. The same numbers of training and testing cases are used with both candidate identifiers.

Given a desired tolerance \(\delta = 0.01\), both FRA-LSGD-DRNN and OLS-LSGD-DRNN can find the suitable number of recurrent neurons from 50 candidate neurons generated initially. After selection process, 10 neurons were added into the network. By (43) and (44), the computational complexity of OLS based constructive method is approximately 10 times higher than FRA, which again validates superior efficiency of the employed algorithm. For the rest of identifiers, the grid-search method can help to find relatively appropriate number of hidden neurons. Figure 8 shows the grid searching process, which is clear to see that 10, 13 and 18 neurons were finally included in the LSGD-DRNN, DRNN and FNN, respectively.

Table 5 Test performance comparison of various identifiers with external input
Full size table

Remark 3

For OLS based DRNN the computational complexity mainly comes from the term \(mN{k^2}\) in (43), and for FRA-based DRNN mainly comes from the term mNk in (44). Therefore, the computational effort of OLS increases exponentially as the selected recurrent neurons increase, while FRA increases slowly.

Fig. 8
figure 8

Testing MSE versus the growth of hidden neurons number for various identifiers

Full size image
Fig. 9
figure 9

HPCMHS identification result of various models

Full size image
Table 6 Performance comparison of various identifiers in HPCMHS (NHY=Number of nodes in hidden layer, TCP=Total count of parameters to be tuned)
Full size table

Figure 9 further illustrates the temperature prediction performance among five identifiers. It can be easily seen that both models have produced a good performance in smooth area. However, the LSGD-based identifiers outperform others in sharp area, where the medium temperature is significantly raised and hot spot may occur. The foreknowledge of hot spots is vital, as in microwave drying process, the quality of grain presented in the hot spot regions may be severely affected due to the elevated temperature. Knowledge of temperature rise when microwaves are used for heating grains would be helpful to select suitable power level and treatment time while developing microwave heating systems for various applications. The various details of five identifiers are listed in Table 6. As with the previous examples, the proposed identifier yields the minimum MSE and MAE value and requires least number of parameters to be tuned as compared with DRNN and FNN. The practical implication is that FRA-LSGD-DRNN can be trained with shorter time and still achieve a proper accuracy level. This suggests a more impact and accurate network model can be used for HPCMHS identification.

6 Conclusion

This paper presents a two-stage method for DRNN based identification. The new design FRA-LSGD-DRNN integrates benefits from the high efficiency of FRA in constructing DRNN model and the improved optimization capability of LSGD in tuning network parameters. Its performance is tested on a simulated nonlinear system with two different inputs and a real-life data measured from microwave rice drying experiment. Results show that (1) compared with OLS and grid-search method, the proposed FRA based constructive method is more computational efficient and the resulting model has fewer parameters to be tuned; (2) superior identification accuracy and robustness can be achieved by using LSGD that is compared with conventional GD principle; (3) DRNN outperforms FNN in nonlinear dynamic system identification due to its information storage ability; (4) the proposed model can better learns the sharp temperature variation of heated medium during microwave heating process, namely it has the ability to predict the emergence of hot spot in advance, which helps to improve safety in actual microwave applications.

Since the HPCMHS is able to continuously generate new experimental data under different operational conditions, its dynamics is changing over time. The adaptive online identification of such non-stationary systems is particularly challenging since our proposed DRNN construction is off-line and fails to be sufficient to track the data change. Further works will focus on the adaptive tunable DRNN for nonlinear and non-stationary systems identification, more powerful and fast structure and parameter optimization approach will be employed to address this issue. Besides, with the advances of sensor technology, multi-point temperature can be measured during microwave heating or drying process, which will also considered in our further work, thus an identification model with higher fidelity and lower complexity will be established for HPCMHS.