An Electric Oven Identification Using Step Response Based on the Prony Approach

1. Introduction

Since the dawn of civilization, mankind has sought ways to simplify life by modifying its environment. The first step in accomplishing this is to understand the environment. For humans, it is crucial first to recognize the environment. This paper aims to identify the model of an oven for the design of its temperature controller. System identification methods can generally be categorized as linear and nonlinear. Within these categories, identification can be performed using gray box[1] or black box[2]1 methods. In the black box identification method, only input-output signal samples are utilized without using any information about the system's internal variables where most of the time, it is difficult or impossible to access this information. A lack of information about the internal components of a system leads to an extensive list of nonlinear identification approaches[3]. Given these explanations, it makes sense to start with a simple structure and then move on to more complex methods if the expected results are not achieved. In[2] and[4] an electric oven is modelled with a first-order linear TF while in[5], an Artificial Neural Network (ANN) approach was used for modelling and temperature control of an electric oven. However, due to the computational complexity of ANNs, this approach may not be suitable for online or practical applications. The authors in[6-7] used energy equations to model an electric oven and control its temperature, focusing on energy efficiency.

In the black - box category of identification methods, the first and simplest methods are identification using impulse and step responses[8]. The Prony is an old method in the literature on systems identification[9-10]. However, it is still one of the most powerful batch2 identification methods[11-12]. This method was examined with two other methods in the review article by[12] using four different case studies. In the review, Prony's method is praised for its performance in order reduction and noise resilience. Properties that can be achieved in the Matrix Pencil (MP) and Eigensystem Realization Algorithm (ERA) methods through the use of auxiliary tools. In the literature on system identification, Prony's method is particularly noteworthy in the identification of power system modes[13-14] and their small signal stability analysis[15], since most of its control methods are based on models. Other applications of Prony's method include the localization of short-circuit faults in DC systems[16], signal recovery[17], and the detection of wave distortion[18] in the field of control systems.

In this paper, the Prony method was chosen because of its ability to identify an oven system. The identification has to be done in the form of a black box due to the lack of information described in the next section, and the Prony method supports this approach. However, implementing the Prony method in practice can be difficult. Prony uses the samples of the impulse response of the system in two stages[19], which cannot be extracted in many practical systems, including the desired oven system. In this category of systems, high-order dynamics and restrictions in the input signal amplitude result in no response at the output. Therefore, the use of step input is recommended for three reasons. First, it does not have the problem mentioned for the impulse input. Secondly, its response can easily be converted into an impulse response suitable for identification. Thirdly, due to the high-frequency range, it fulfills the persistent excitation condition of the input.

In general, the main contribution of this paper is to reconcile Prony's theory with the practical realities of implementation. This can be summarized as follows: first, sampling the output of an industrial system using common tools; second, utilizing step and pulse (rectangular) inputs in two different stages of the Prony method instead of an impulse input; third, updating the Prony method relations to parameterize the step response as original data and using a rectangular signal as secondary data; Fourth, convert the parameterized step response into discrete and continuous transfer functions of the oven system; and finally, using a valid criterion to select a model with the appropriate order in terms of fitting and parsimony.

2. The oven system under study

The oven system, under study in this paper, located in the robotic laboratory, is shown in Figure 1. The system is powered by household current (220V AC and 50Hz) and the main box temperature is the system output. Household current is supplied to the heater through a Solid-State Relay (SSR) that can be controlled by a controller. Because the heater is in the intermediate chamber outside the main box, the oven system experiences dead time and high-order dynamics.

Specifically, the heater is positioned in the middle layer on the bottom, left, and back sides. It is insulated from the back with fiberglass and is insulated (from the ambient temperature) on the bottom and left side with one and two layers, respectively. The top and right sides have no heater and are connected to the outside (the ambient temperature) with two layers, and it is the same for the oven door on the front side with a different conduction coefficient. As the transfer coefficients are unknown and the heater in the middle layer does not provide a uniform temperature, the black box identification method is a suitable approach to determine the mathematical model of the oven.

Fig. 1 Oven system located in the robotics Lab, University of Mohaghegh Ardabili

The oven system under study is equipped with a PT100 sensor located inside the main housing, which records samples via the Arduino board. However, due to physical limitations, step response samples were used instead of impulse responses. To protect the system's circuitry from possible damage, only 20% of the maximum input value (denoted by 0.2u) is applied instead of 100% (full step input). The recorded step response for the oven system under study is shown in Figure 2.

Fig. 2 Step response of the oven system to 0.2u

In Figure 2, approximately 7199 samples of the output (oven temperature) were recorded with a sampling time of one second. Additionally, in order to initiate the step response from zero, an initial temperature of 22.7°C was deducted from the recorded data. According to this figure, the final temperature value is 73.8°C.

3. Prony's method principles and its adaptation to the desired system

It is assumed that the TF of the system under study can be represented as a rational fractional expression with two polynomials with [TeX:] b_{i}$b_i$ real coefficients in the numerator and [TeX:] a_{i}$a_i$ in the denominator. In other words, the TF of the system under study is assumed to be as shown in Eq. (1).

[TeX:] G\left( s \right) = \frac{b_{0} + b_{1}s + ... + b_{m}s^{m}}{a_{0} + a_{1}s + ... + a_{n}s^{n}} = \sum_{i = 1}^{k}{\frac{A_{i}}{s + \alpha_{i}} + \sum_{i = k + 1}^{\frac{n + k}{2}}\frac{D_{i}\omega_{\text{ni}}^{2}}{s^{2} + 2\xi_{i}\omega_{\text{ni}}s + \omega_{\text{ni}}^{2}}}$G\left(s\right)=\frac{b_0+b_{1}s+...+b_m{s}^m}{a_0+a_{1}s+...+a_n{s}^n}={\sum }_{i=1}^k\frac{A_i}{s+{\alpha }_i}+{\sum }_{i=k+1}^{\frac{n+k}{2}}\frac{D_i{\omega }_{\mathrm{\text{ni}}}^2}{s^2+2{\xi }_i{\omega }_{\mathrm{\text{ni}}}s+{\omega }_{\mathrm{\text{ni}}}^2}$ (1)

Where [TeX:] \alpha_{i}${\alpha }_i$s are k real poles separated from complex conjugate ones and [TeX:] A_{i}$A_i$ and [TeX:] D_{i}$D_i$ are the values of the residues in the corresponding poles. Here it is assumed that the system has no repeated (multiple/double) poles and that it is separated in the form of distinct real and complex-conjugate pairs as in Eq. (1). Combining the two summations given in Eq. (1), a more unified form of this TF is given in Eq. (2).

[TeX:] G\left( s \right) = \frac{b_{0} + b_{1}s + ... + b_{m}s^{m}}{a_{0} + a_{1}s + ... + a_{n}s^{n}} = \sum_{i = 1}^{n}\frac{A_{i}}{s + \left( \sigma_{i} + j\omega_{d_{i}} \right)}$G\left(s\right)=\frac{b_0+b_{1}s+...+b_m{s}^m}{a_0+a_{1}s+...+a_n{s}^n}={\sum }_{i=1}^n\frac{A_i}{s+({\sigma }_i+j{\omega }_{d_i})}$ (2)

when compared with Eq. (1), for the first k terms, [TeX:] \omega_{d_{i}}${\omega }_{d_i}$ is zero and [TeX:] \sigma_{i} = \alpha_{i}${\sigma }_i={\alpha }_i$.

By applying the step input and using the inverse Laplace transform, the step response of the desired system can be written as Eq. (3).

[TeX:] S\left( s \right) = \frac{1}{s}G\left( s \right) = \frac{1}{s}\frac{b_{0} + b_{1}s + ... + b_{m}s^{m}}{a_{0} + a_{1}s + ... + a_{n}s^{n}} = \sum_{i = 1}^{n + 1}\frac{B_{i}}{s + \left( \sigma_{i} + j\omega_{d_{i}} \right)}$S\left(s\right)=\frac{1}{s}G\left(s\right)=\frac{1}{s}\frac{b_0+b_{1}s+...+b_m{s}^m}{a_0+a_{1}s+...+a_n{s}^n}={\sum }_{i=1}^{n+1}\frac{B_i}{s+({\sigma }_i+j{\omega }_{d_i})}$ (3)

[TeX:] \Rightarrow s\left( t \right) = \sum_{i = 1}^{n + 1}{B_{i}e^{- \left( \sigma_{i} + \omega_{d_{i}} \right)t}}$\Rightarrow s\left(t\right)={\sum }_{i=1}^{n+1}{B}_i{e}^{-({\sigma }_i+{\omega }_{d_i})t}$

In which [TeX:] \sigma_{n + 1}${\sigma }_{n+1}$ and [TeX:] \omega_{d_{n + 1}}${\omega }_{d_{n+1}}$ are formed due to the step input and both are equal to zero. Assuming that the TF given in Eq. (1) can appropriately model the oven system, the recorded samples in Figure 2 must be fitted in the time-domain equation given in Eq. (3). However, since the sampling rate is equal to a constant value such as T, the sample's exact values can be represented as Eq. (4).

[TeX:] s_{k} = \sum_{i = 1}^{n + 1}B_{i}e^{- \left( \sigma_{i} + j\omega_{d_{i}} \right)\text{kT}} = \sum_{i = 1}^{n + 1}B_{i}z_{i}^{k}\ \ \ k = 0,1,2,...,N - 1$s_k={\sum }_{i=1}^{n+1}{B}_i{e}^{-({\sigma }_i+j{\omega }_{d_i})\mathrm{\text{kT}}}={\sum }_{i=1}^{n+1}{B}_i{z}_i^{k}k=0,1,2,...,N-1$ (4)

where [TeX:] z_{i} = e^{- \left( \sigma_{i} + j\omega_{d_{i}} \right)T}$z_i=e^{-({\sigma }_i+j{\omega }_{d_i})T}$, k is the sample number, and N is the total number of samples.

In the Prony method, the goal is to find [TeX:] z_{i}$z_i$s and [TeX:] B_{i}$B_i$s in such a way that [TeX:] s_{k}$s_k$ matches the step response samples. Expanding Eq. (4) for different samples leads to Eq. (5).

[TeX:] \begin{bmatrix} B_{1}z_{1}^{0} + ... + B_{n + 1}z_{n + 1}^{0} \\ B_{1}z_{1}^{1} + ... + B_{n + 1}z_{n + 1}^{1} \\ \vdots \\ B_{1}z_{1}^{N - 1} + ... + B_{n + 1}z_{n + 1}^{N - 1} \\ \end{bmatrix} = \begin{bmatrix} z_{1}^{0} & z_{2}^{0} & \cdots & z_{n + 1}^{0} \\ z_{1}^{1} & z_{2}^{1} & \cdots & z_{n + 1}^{1} \\ \vdots & \vdots & \vdots & \vdots \\ z_{1}^{N - 1} & z_{2}^{N - 1} & \cdots & z_{n}^{N - 1} \\ \end{bmatrix}\begin{bmatrix} B_{1} \\ B_{2} \\ \vdots \\ B_{n + 1} \\ \end{bmatrix} = \begin{bmatrix} s_{0} \\ s_{1} \\ \vdots \\ s_{N - 1} \\ \end{bmatrix}$\left[\begin{array}{c}\hfill B_1{z}_1^0+...+B_{n+1}{z}_{n+1}^0\hfill \\ \hfill B_1{z}_1^1+...+B_{n+1}{z}_{n+1}^1\hfill \\ \hfill ⋮\hfill \\ \hfill B_1{z}_1^{N-1}+...+B_{n+1}{z}_{n+1}^{N-1}\hfill \end{array}\right]=\left[\begin{array}{cccc}\hfill z_1^0\hfill & \hfill z_2^0\hfill & \hfill \cdots \hfill & \hfill z_{n+1}^0\hfill \\ \hfill z_1^1\hfill & \hfill z_2^1\hfill & \hfill \cdots \hfill & \hfill z_{n+1}^1\hfill \\ \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill ⋮\hfill \\ \hfill z_1^{N-1}\hfill & \hfill z_2^{N-1}\hfill & \hfill \cdots \hfill & \hfill z_n^{N-1}\hfill \end{array}\right]\left[\begin{array}{c}\hfill B_1\hfill \\ \hfill B_2\hfill \\ \hfill ⋮\hfill \\ \hfill B_{n+1}\hfill \end{array}\right]=\left[\begin{array}{c}\hfill s_0\hfill \\ \hfill s_1\hfill \\ \hfill ⋮\hfill \\ \hfill s_{N-1}\hfill \end{array}\right]$ (5)

[TeX:] \mathbb{Z}_{N \times \left( n + 1 \right)}\mathbb{B}_{\left( n + 1 \right) \times 1} = \mathbb{S}_{N \times 1}${\mathbb{Z}}_{N\times (n+1)}{𝔹}_{(n+1)\times 1}={𝕊}_{N\times 1}$

There are two challenges when determining the unknown parameters of Eq. (5).

To overcome the first challenge, more information (perhaps more recorded samples) is required. According to Eq. (4), the components of [TeX:] \mathbb{Z}$\mathbb{Z}$ are the poles of the desired system. These poles can be determined by finding the roots of the denominator polynomial of the TF in Eq. (1). For this purpose, it is necessary to put additional samples in the different equations of the system given in Eq. (6) and obtain [TeX:] a_{1}$a_1$ to [TeX:] a_{n}$a_n$ coefficients. The relation Eq. (6) can be directly derived from Eq. (1) where [TeX:] x_{i}$x_i$s and [TeX:] y_{i}$y_i$s are the samples of the input and output signals, respectively.

[TeX:] G\left( z \right) = \frac{Y\left( z \right)}{X\left( z \right)} = \frac{b_{0} + b_{1}z^{- 1} + ... + b_{m}z^{- m}}{1 + a_{1}z^{- 1} + .... + a_{n}z^{- n}}$G\left(z\right)=\frac{Y\left(z\right)}{X\left(z\right)}=\frac{b_0+b_1{z}^{-1}+...+b_m{z}^{-m}}{1+a_1{z}^{-1}+....+a_n{z}^{-n}}$ (6)

[TeX:] y_{i} + a_{1}y_{i - 1} + ... + a_{n}y_{i - n} = b_{0}x_{i} + b_{1}x_{i - 1} + ... + b_{m}x_{i - m}$y_i+a_1{y}_{i-1}+...+a_n{y}_{i-n}=b_0{x}_i+b_1{x}_{i-1}+...+b_m{x}_{i-m}$

This can be useful when an input is used so that the right side of the difference equation given in Eq. (6) becomes equal to zero. In Prony’s method, the appropriate impulse signal is one and for [TeX:] i = n$i=n$ to N-1 the samples are placed in Eq. (6)[20]. However, as mentioned above, applying the ideal impulse input to the oven system does not significantly affect its output (temperature), i.e., cannot be measured with a sensor. To solve this challenge, we used a pseudo-impulse signal that equals 1 from 0 to 500 seconds it is equal to 1 (100% duty cycle) and then becomes zero (a rectangular pulse with a width of 500 seconds). Since the sampling time is 1 second, Eq. (6) can be rewritten for [TeX:] i = 500 + n$i=500+n$ to [TeX:] i = N - 1$i=N-1$ as follow:

[TeX:] \begin{bmatrix} y_{500 + n} \\ y_{501 + n} \\ \vdots \\ \vdots \\ y_{N - 1} \\ \end{bmatrix} = \begin{bmatrix} y_{499 + n} & y_{498 + n} & \cdots & \cdots & y_{500} \\ y_{500 + n} & y_{499 + n} & \cdots & \cdots & y_{501} \\ \vdots & \ & \ddots & \ & \ \\ \vdots & \ & \ & \ddots & \ \\ y_{N - 2} & \ & \ & \ & y_{N - 1 - n} \\ \end{bmatrix}\begin{bmatrix} - a_{1} \\ - a_{2} \\ \vdots \\ \vdots \\ - a_{n} \\ \end{bmatrix},\ \ \ N > 501 + n$\left[\begin{array}{c}\hfill y_{500+n}\hfill \\ \hfill y_{501+n}\hfill \\ \hfill ⋮\hfill \\ \hfill ⋮\hfill \\ \hfill y_{N-1}\hfill \end{array}\right]=\left[\begin{array}{ccccc}\hfill y_{499+n}\hfill & \hfill y_{498+n}\hfill & \hfill \cdots \hfill & \hfill \cdots \hfill & \hfill y_{500}\hfill \\ \hfill y_{500+n}\hfill & \hfill y_{499+n}\hfill & \hfill \cdots \hfill & \hfill \cdots \hfill & \hfill y_{501}\hfill \\ \hfill ⋮\hfill & \hfill \ddots \hfill \\ \hfill ⋮\hfill & \hfill \ddots \hfill \\ \hfill y_{N-2}\hfill & \hfill y_{N-1-n}\hfill \end{array}\right]\left[\begin{array}{c}\hfill -a_1\hfill \\ \hfill -a_2\hfill \\ \hfill ⋮\hfill \\ \hfill ⋮\hfill \\ \hfill -a_n\hfill \end{array}\right],N>501+n$ (7)

In Eq. (7), the second challenge of the matrix relation of Eq. (5) appears again. This challenge appears since the number of equations is greater than the number of unknown parameters, and therefore the equations cannot have a unique solution. To solve this, it is common to use the pseudo-inverse formula in which a suitable unknown parameters vector that holds for almost all components is obtained in Eq. (7)[21]. For this:

[TeX:] \mathbb{\text{D.}}\mathbb{= Y}$\mathrm{\text{𝔻.a ̂}}=𝕐$; [TeX:] [TeX:] \mathbb{D}^{T}\mathbb{\text{D.}} = \mathbb{D}^{T}\mathbb{Y}${𝔻}^T\mathrm{\text{𝔻.a ̂}}={𝔻}^T𝕐$; [TeX:] [TeX:] = \left( \mathbb{D}^{T}\mathbb{D} \right)^{- 1}\mathbb{D}^{T}\mathbb{Y}$a\; ̂={\left({𝔻}^T𝔻\right)}^{-1}{𝔻}^T𝕐$ (8)

Where [TeX:] \mathbb{D}$𝔻$ is the output sample matrix in Eq. (7) and [TeX:] \mathbb{Y}$𝕐$are samples (5500 samples) in Eq.(8), as illustrated in Figure 3. It should be noted that the initial value of temperature is subtracted from this signal. These samples present the response of the oven system to the pseudo-impulse input (the rectangular signal in the interval time between 0 to 500 seconds).

Fig. 3 The system response to the pseudo-impulse signal

As [TeX:] \widehat{\mathbb{Z}}$\hat{\mathbb{Z}}$ is determined using Eq. (5), the unknown vector of [TeX:] \widehat{\mathbb{B}}$\widehat{𝔹}$ can be obtained using a pseudo-inverse formula as given in Eq. (9). It should be noted that a pole at [TeX:] z = - 1$z=-1$ is added to the overall poles due to the applied step input signal.

[TeX:] {\widehat{\mathbb{B}}}_{\left( n + 1 \right) \times 1} = \left( {\widehat{\mathbb{Z}}}_{\left( n + 1 \right) \times N}^{T}{\widehat{\mathbb{Z}}}_{N \times \left( n + 1 \right)} \right)^{- 1}{\widehat{\mathbb{Z}}}_{\left( n + 1 \right) \times N}^{T}\mathbb{S}_{N \times 1}${\widehat{𝔹}}_{(n+1)\times 1}={\left({\hat{\mathbb{Z}}}_{(n+1)\times N}^T{\hat{\mathbb{Z}}}_{N\times (n+1)}\right)}^{-1}{\hat{\mathbb{Z}}}_{(n+1)\times N}^T{𝕊}_{N\times 1}$ (9)

Finally, it can be concluded that:

[TeX:] {\widehat{s}}_{k} = \sum_{i = 1}^{n + 1}{\widehat{B}}_{i}{\widehat{z}}_{i}^{k}\ \ \ k = 0,1,2,...,7198${\hat{s}}_k={\sum }_{i=1}^{n+1}{\hat{B}}_i{\hat{z}}_i^{k}k=0,1,2,...,7198$ (10)

Regarding Eq. (10), the system order (n) has not been determined yet. This parameter is available to the designer.

Applying Z-transform to Eq. (10), the discrete TF given in Eq. (6) can be obtained as follows:

[TeX:] G\left( z \right) = \left( \sum_{i = 1}^{n + 1}\frac{{\widehat{B}}_{i}}{1 - {\widehat{z}}_{i}z^{- 1}} \right)\frac{z - 1}{0.2z}$G\left(z\right)=\left({\sum }_{i=1}^{n+1}\frac{{\hat{B}}_i}{1-{\hat{z}}_i{z}^{-1}}\right)\frac{z-1}{0.2z}$ (11)

According to the details considered in Eq. (4) and using [TeX:] s_{i} = \sigma_{i} + j\omega_{i} = - \frac{\ln z_{i}}{T}$s_i={\sigma }_i+j{\omega }_i=-\frac{ln{z}_i}{T}$, the continuous TF is obtained by converting discrete TF poles to continuous TF poles. Here we get Eq. (12):

[TeX:] G\left( s \right) = \left( \sum_{i = 1}^{n + 1}\frac{{\widehat{B}}_{i}}{s - \frac{{\widehat{\ln z}}_{i}}{T}} \right)\frac{s}{0.2}$G\left(s\right)=\left({\sum }_{i=1}^{n+1}\frac{{\hat{B}}_i}{s-\frac{{\widehat{lnz}}_i}{T}}\right)\frac{s}{0.2}$ (12)

For the system order (n) to be determined, the following mean squared error (MSE) criterion is utilized. For this purpose, the obtained model is evaluated for different values of the n, and the order that results in the lowest value of MSE is chosen as the final system model order.

[TeX:] MSE = \frac{1}{N}\left( \sum_{k = 0}^{N - 1}\left( s_{k} - {\widehat{s}}_{k} \right)^{2} \right)$MSE=\frac{1}{N}\left({\sum }_{k=0}^{N-1}{(s_k-{\hat{s}}_k)}^2\right)$ (13)

4. Order finding and simulation

4.1 Implementation of Prony's algorithm

For the models with different orders to be evaluated, the proposed identification method was tested with different inputs and the Mean Squared Error (MSE) results are summarized in Table 1. The identification of the oven system was performed with a 0.2u input signal, but the MSE can also be used to evaluate a model at different operating points. For this purpose, 0.4u and 0.6u inputs were applied to the obtained models, and the MSE was calculated. The results of all evaluations are shown in Table 1. To get a better understanding of the modelling accuracy, the responses of different order models to a 0.2u input are shown in Figure 4.

Table 1 Mean Square Errors of different order models

Fig. 4 Step response of the identified oven system with different orders (identified using 0.2u input)

Using Eqs. (11) and (12), the parameters obtained in Prony’s curve fitting method can be utilized to get the discrete and continuous transfer functions of the oven system. For example, for [TeX:] n \in \left\{ 3,5,10 \right\}$n\in \{3,5,10\}$ the continuous transfer functions are extracted as follows:

The third-order model:

[TeX:] G\left( s \right) = \frac{- 8.11s^{3} + 0.07s^{2} + 2.82s + 0.004}{0.2s^{3} + 0.16s^{2} + 0.03s + 0.000007} = \frac{- 75.2\left( s - 0.59 \right)\left( s + 0.58 \right)\left( s + 0.001 \right)}{\left( s + 0.5 \right)\left( s + 0.3 \right)\left( s + 0.0003 \right)}$G\left(s\right)=\frac{-8.11s^3+0.07s^2+2.82s+0.004}{0.2s^3+0.16s^2+0.03s+0.000007}=\frac{-75.2(s-0.59)(s+0.58)(s+0.001)}{(s+0.5)(s+0.3)(s+0.0003)}$ (14)

The fifth-order model:

[TeX:] G\left( s \right) = \frac{- 11.27\left( s + 1.42 \right)\left( s - 0.06 \right)\left( s + 0.001 \right)\left( s^{2} - 2.84s + 6.22 \right)}{\left( s + 0.33 \right)\left( s + 0.05 \right)\left( s + 0.0003 \right)\left( s^{2} + 0.6s + 3.62 \right)}$G\left(s\right)=\frac{-11.27(s+1.42)(s-0.06)(s+0.001)(s^2-2.84s+6.22)}{(s+0.33)(s+0.05)(s+0.0003)(s^2+0.6s+3.62)}$ (15)

The tenth-order model:

[TeX:] G\left( s \right) = \frac{- 11.52\left( s - 24.61 \right)\left( s - 1.04 \right)\left( s - 0.01 \right)\left( s + 0.003 \right)\left( s^{2} - 0.12s + 1.01 \right)\left( s^{2} - 0.01s + 3.16 \right)\left( s^{2} + 0.1s + 6.25 \right)}{\left( s + 0.0005 \right)\left( s + 0.009 \right)\left( s^{2} + 0.36s + 0.9 \right)\left( s^{2} + 0.42s + 2.6 \right)\left( s^{2} + 0.46s + 4.9 \right)\left( s^{2} + 0.48s + 8.1 \right)}$G\left(s\right)=\frac{-11.52(s-24.61)(s-1.04)(s-0.01)(s+0.003)(s^2-0.12s+1.01)(s^2-0.01s+3.16)(s^2+0.1s+6.25)}{(s+0.0005)(s+0.009)(s^2+0.36s+0.9)(s^2+0.42s+2.6)(s^2+0.46s+4.9)(s^2+0.48s+8.1)}$ (16)

4.2 Akaike information criterion

There are two important principles in determining the model order: the fitting principle and the parsimony principle. The fitting principle focuses exclusively on minimizing the MSE, which often results in a high-order model. While this type of model may be suitable for purposes such as fault detection, if the goal of the modelling is to design a controller for the system, a lower-order model is sufficient and simplifies the calculations. In contrast, according to the principle of parsimony, a lower order model is chosen even if it is weaker in terms of fitting. Among the practical methods for determining the model order determination, the variance method, the covariance matrix method[22], the Akaike information criterion (AIC) and the f - test method are frequently used[22]. In this paper, we use the AIC method due to its efficiency and simplicity in calculation[23].

According to the simulations (Table 1), for models with [TeX:] n > 27$n>27$, due to the unnecessary order, the MSE increases, and the model loses accuracy. Also, choosing a model with an order of 27 is solely based on the fitting principle which is considered a large order for a model. The AIC is used to determine the final oven system model considering both the principles of parsimony and fitting.

[TeX:] \left( \text{AIC} \right)_{n} = 10\log\left( \sum_{k = 0}^{N - 1}\left( s_{k} - {\widehat{s}}_{k} \right)^{2} \right) + p${\left(\mathrm{\text{AIC}}\right)}_n=10log\left({\sum }_{k=0}^{N-1}{(s_k-{\hat{s}}_k)}^2\right)+p$ (17)

where p=n+m+1. The value of AIC for [TeX:] n \in \left\{ 1,...,50 \right\}$n\in \{1,...,50\}$ is calculated and the results are illustrated in Figure 5.

Fig. 5 Akaike information function values in Eq. (17) for models with orders from 1 to 50

The AIC has its minimum value for n=17, which is 128.34. Based on this criterion and considering both the fitting and parsimony principles, the model with n=17 is therefore the most appropriate model for the oven system.

The overall process using the Akaike information function to obtain a linear model with an order of 17 can be summarized as shown in Figure 6. It should be noted that Figure 5 represents the evaluation result within the diamond shown in the flowchart of Figure 6. This figure shows that the decision-making process to determine the appropriate model is inductive and not straightforward.

Fig. 6 Flowchart of the model order finding process

5. Discussion

Due to the presence of noise in the input/output sampled data, the statistical approach to parameter identification has been investigated in several references[24-25]. Most of the analyses in this area refer to the Least Squares (LS) method, which uses linear regression equations to minimize the sum of squared errors and determine system parameters[26]. This method is commonly used to identify system parameters by identifying the coefficients of impulse responses or difference equations[19]. According to the BLUE theorem, unbiased estimation is achieved when the noise in the sampled data is white, making this method the best linear estimator. However, white noise is a theoretical concept and can only be approximated in practice, just like impulse signals. The justification for using the LS method in its ordinary format lies in the ignorable error caused by non-white (color) noise. In this paper, Prony's method is also selected for the same reasons. In the following, the initial steps taken to reduce the level and impact of noise to increase the success rate of the Prony method implementation are outlined.

One way to reduce the effect of noise is to filter the sampled data with low-pass filters, which is done by programming. Another measure to reduce the effects of noise is to use an RTD temperature sensor instead of a thermocouple. Thermocouples can have a very high noise level as they generate very low voltages.

Another important point is that in Prony's method (also in the LS method), the last equation to be solved (here Eq. (7) and then Eq. (5)) actually represents a system of equations with more equations than unknowns. Solving this system using the transpose matrix in Eqs. (8) and (9), the so-called pseudo-inverse, means that the fitting operation is performed twice in the Prony method, similar to the LS method.

The next point to discuss is the error analysis in Table 1. When a linear model is used to model a nonlinear system, the desired model is expected to be valid only at the point of operation where the samples were taken. Too much displacement will invalidate the model, but this does not mean that we should abandon the advantages of linear models. As can be seen in Figure 7, the main difference between the step responses obtained with 20%, 40%, and 60% inputs is one coefficient. In such cases, it is easy to use new linear models by means of the Gain scheduling technique[27], so that when the operation point of the system changes, the corresponding coefficient in the linear model changes.

Fig. 7 The response of the oven system to 20%, 40% and 60% of the step input

6. Conclusion

This paper modelled an industrial electric oven using the Prony-based curve-fitting method and AIC. The Prony method can utilize all the samples from the system's output in its calculations by employing pseudo-inversion. This method has been demonstrated to provide the best linear approximation for the data used. Applying step inputs with amplitudes of 20%, 40%, and 60% of the unit step to the models under evaluation from different orders in Table 1 and the results obtained in Figure 7 confirmed that the model calculated with the 0.2u input data can perform well for other step inputs simply by adjusting the DC gain of the TF. Using AIC designers can strike a balance between the principles of fitting and parsimony, and select the AIC coefficients based on the modelling purpose. Finally, it can be concluded that based on the desired results obtained in Figure 4, the need to use more complex models, such as non-linear ones, is unnecessary.