Adaptive identification method for simulation and control of glass melting process

Adaptive identification method based on modulating function approach is described in this paper. The real glass melting process is modelled by a linear system, which parameters values can be changed on-line during the identification procedure. The overall model is composed of several submodels. The models, which inputs can be controlled, are approximated by the Strejc transfer function. This approach allows to reduce the number of model parameters. Moreover, PID tuning methods dedicated to the Strejc models can be applied. The developed procedure was applied for data collected from a real glass containers production installation.


I. INTRODUCTION
The paper presents an on-line system identification method applied for a glass melting forehearth temperature control system. Maintaining an appropriate temperature profile in the forehearth is crucial for final products quality. Two main problems are often considered for glass containers production process. The first one concerns finding optimal temperatures for each segment (zone) of the forehearth depending on the current glass mass flow (pull rate) in a steady state. For this purpose nonlinear models based on partial differential equations are usually applied. These models can be very useful for simulation purposes, but often not convenient for real-time control. The problem of their synthesis was discussed in [5] and [7]. An example of Finite Element Method (FEM) based model synthesis for glass melting process was presented in [6]. The second issue, addressed in the paper, involves finding an optimal control strategy of transitions between two steady states. Authors in both [5] and [7] performed simplified models synthesis based on previously obtained results. The motivation of this paper is similar, but in contrast to the mentioned papers, on-line process identification idea is assumed. The process dynamics for a single forehearth section is approximated with a linear model near a selected operating point. The model parameters can be obtained with the use of Modulating Function Method (MFM). Different variants of this method can be found in works: [2], [3], [4].
The identification approach applied in the article can be divided into two stages -passive and active, which results from the fact that performing a step response experiment is often inadmissible for real industrial installations. Firstly, the linear model identified with the use of MFM, is used for simulating the system step response. Subsequently, the Strejc model is obtained based on this response and formulas for optimal PID controller parameters can be utilized.
The paper outlines the following sections. Section I briefly introduces the analysed problem. Sections II and III discuss theoretical basis of the implemented identification methods. Developed algorithm is presented in details in Section IV. Section V shortly describes the glass containers production installation. The results of performed experiments and their discussion are given in Sections VI and VII.

II. MODULATING FUNCTIONS METHOD
For Single Input-Single Output (SISO) Linear Time-Invariant (LTI) model (1) n + 1 parameters a i and m + 1 parameters b i should be identified. It is assumed that signals u(t) and y(t) are measured directly. Their derivatives can be calculated numerically, but the results can be inaccurate especially for high order signal derivatives in case of noisy measurements.
MFM assumes the use of function inner product with the modulating function φ(t). In the approach described in [3] and applied in this paper, single modulating function φ(t) is convoluted with each term of the differential equation: It is assumed that functions φ(t) and φ (i) (t) values are given on the interval [0, h]. Function φ(t) has the compact support and below conditions must be met: Obtained continuous functions y i (t) and u i (t) are defined for the interval [h, T ID ], where T ID is the identification interval length. Differential equation (1) can be transformed into the algebraic one with the same parameters a, b. Identification results significantly depend on the chosen modulating function.
In the described case, the Loeb and Cahen function: was chosen. In the developed algorithm, the MFM is applied for the MISO system with K inputs and a single output (5).
After performing the described convolution transformation for (5), algebraic equation (6) can be obtained: The term (t) represents a difference resulting from signal noises and modelling errors. It can be treated as a performance index of the identification procedure.
Linear constraint vector (12) is introduced to avoid the trivial solution: Hence, the minimisation task is solved with the use of Lagrange multiplier technique: Solution of (13) has the form: Problem of the optimal constraint vector θ choice was discussed in [2]. It was proven that the minimal value of the performance index J can be obtained for the eigenvector of the matrix G corresponding with its minimal eigenvalue.

III. THE STREJC MODEL IDENTIFICATION
The Strejc model identification method was the first time proposed in [10]. This relatively simple method of model identification based on its step response is widely used. Different ways of obtaining Strejc model can be found in papers: [11], [12], [13]. Strejc model transfer function with delay has the form: New methodology of the Strejc model identification was developed in [1]. It is based on the assumption that a normalised step response of an unknown system should be equal to a normalised step response of the identified Strejc model in two chosen points of these two characteristics. The first one is the inflection point t 0 and the second point is arbitrary chosen. In contrast to the commonly used graphic procedures of finding the inflection point for the system step response, analytical formulas for any n ≥ 1 were developed. In practice, it may occur that the chosen system rank is to high for the system and valid parameters values cannot be obtained for the specified n. Experiments for noisy measurements were performed in [1] and confirmed the method effectiveness and its filtering properties.
Transfer function for the Strejc model without delay can be expressed as: The step response for this model is given by the formula: In Section II the letter h was used for denoting the support of the function φ(t). From now, the function h(t) will indicate the Strejc model step response. This should not cause a misunderstanding. The inflection point of function h(t) can be obtained after calculating the second derivative as below: The obtained non-trivial solution has the form: Substituting (19) into (17), the formula for the inflection point value h(T p ) can be obtained: The inflation points values h n for different ranks of the Strejc model, assuming normalized step response characteristics, where h n (T p ) = h(Tp(n)) k are presented in Table I. During the identification procedure, the experimental object characteristics h 0 is used. It is assumed that this characteristics is also normalised to 1. Developed identification procedure is presented as Algorithm 1.
The equation (21) was solved for n = 1, . . . , 5 in [1]. Values of T n and τ n are presented in Table I. Identification experimental results are presented in the following part of this paper. Algorithm 1 Strejc model identification procedure.
Step 1. Select the desired Strejc model rank.
Step 2. Determine the value of the time t 0n for which the value h 0 (t 0n ) = h n .
Step 3. Determine the value of the time T 90 for which the value h 0 (T 90 ) = 0.9.
Step 4. Obtain the values T n and τ n assuming τ n = t 0n − T p and solving the below equation:

IV. PROPOSED IDENTIFICATION METHOD
As it was pointed out in Introduction, the aim of the described identification procedure is to obtain the linear MISO model, which is able to describe the process dynamics near a selected operating point. The overall model structure is presented in Figure 1. Model parameters are identified for the determined intervals and can be changed during the identification procedure as a result of changing working conditions. Simulated output of the whole MISO system can be obtained as a sum of SISO model outputs. Single SISO LTI system dynamics for k-th input can be modelled with the use of the transfer function (22).
Individual SISO models are identified in two ways during the identification procedure: • As a part of the MISO system according to the procedure described in Section II. All identified submodels have then the same denominator.
• As a separate SISO system using the appropriate system input and simulated SISO system output. State-space system representation was introduced to perform simulation of the models with updated parameters starting with nonzero initial conditions. State equation for the k-th SISO submodel is presented as (23). The matrix C k has full rank, thanks to which the whole simulated system state can be observed and utilized as an initial condition for the system simulation in successive time intervals.
In order to compare the identified model response with the real system output, the additional performance index (24), except (8), was introduced. It enables to define the difference between the real system output y(t) and the simulated MISO model output y s (t). In the proposed method, parameters of the selected SISO or the whole MISO model are updated, if the model with less performance index value can be obtained in the last identification interval. Detailed method description is presented as Algorithm 2. Simulated k-th SISO system output is denoted as y sk . For Algorithm 2 t 0 denotes the plant operating point, in contrast to Algorithm 1, where this indication was reserved for the transfer function inflection point.
V. GLASS FOREHEARTH CONTROL SYSTEM Described identification algorithm was tested for data collected from a real glass containers production installation. Forehearth is a part of the installation, where molten glass temperature is being stabilised. This process has a big impact on final products quality. Typical forehearth is divided into few zones, as shown in Figure 2. Each zone, except spout, has its own temperature controller. Calculated control signal is an input for installed actuators, such as mixture pressure valves, cooling air valves and cooling dampers. Temperature in each zone is measured with thermocouples in range of 900-1400 o C. Temperature values depends on produced container type. During the production they should be stabilised with accuracy up to 1 o C. The data used in described experiments were collected for the last forehearth zone before spout. This zone is equipped with a single PID temperature controller adjusting mixture pressure. Maintaining the desired temperature set point in this section is particularly important, because it is the closest zone to the containers forming machine. Algorithm 2 MISO model identification procedure.
Step 1. Find operating point and perform system linearization in this point.
Step 2. Perform MISO model identification starting from t 0 to t 0 + T start .
Step 3. Perform system simulation starting from t 0 to t 0 + T start + T . Save state variables values x(t 0 + T start + T ) as the initial condition for the next interval. Set the interval counter p = 2.
Step 4. Calculate the performance index value in the current interval E cur according to (24).
Step 5. Check if E curr < tr, where tr is the threshold value. If the condition is met, do not change the MISO model parameters and go Step 8. Otherwise, go to Step 6.
Step 6. For each k = 1, . . . , K perform the SISO model identification, for the input signal u k (t) and the output signal y sk (t) = y(t) − j=1,...,K:j =k y sj (t) in the last interval. Calculate the performance index E k for each SISO model. Perform also MISO model identification in the last interval and calculate the performance index value E K+1 for this model.
update the model parameters.
Step 8. Increment the interval counter p. Perform forward simulation for the current interval with the initial condition x 0 .
Step 9. Update the model initial condition for the next interval x 0 = x(t 0 + T start + T · (p − 1)).

VI. EXPERIMENTAL RESULTS
Historical data used in the experiments are presented in Figure 3. Mixture pressure in the controlled zone is treated as a system input, while molten glass temperature measured in the previous forehearth zone is treated as a measured disturbance.
In the first experiment two models were identified: • SISO 1 , where: u 1 -air-gas mixture pressure.
• SISO 2 , where: u 2 -molten glass temperature in the previous forehearth zone. Temperature measured in the current forehearth zone is treated as the real system output y. Applied models coefficients and parameters of the MFM are presented in Table II. Parameters specified for the identification procedure described in Section IV are presented in Table III.  In Figure 4 the experimental results are presented and compared with the historical data. Intervals in which changed models were applied for the first time for simulation are marked as below: • B -both models were changed.
• T -SISO 2 model parameters were changed. Model parameters identified during the simulation procedure are presented in Tables IV-V.   In the second experiment, identified models for SISO 1 were used in order to obtain the Strejc models according to the procedure described in Section III. Identified Strejc model parameters are presented in Table VI. Results of the identification procedure in both cases are compared in Table VII. Mean square error between the real system output and the model is assumed as a performance index. Based on the obtained Strejc models, PID controller parameters were adjusted according to the procedure presented by H. S. Preuss in [9]. The algorithm utilizes the system openloop step response. Strejc transfer function is assumed for the model. The method is based on the below criterion: where G z is the closed-loop system transfer function and ω r is the maximum considered angular frequency value for the system. This approach was implemented by Siemens in the series of SIPART temperature controllers. Formulas for PI and PID controllers can be found in [8]. In the experiment, formulas for PI controller were used, because there are no PID controller formulas for n = 2. The controller transfer function is given as R(s): (26) Parameters values for the defined simulation intervals are presented in Table VIII.
Results of the simulation with tuned PID controller parameters are presented in Figure 5. It can be seen that the controller was able to track the changes of the operating point, except the situation when the controller output reaches the lower limit. Analogously as in the real control system, controller output (mixture pressure) was limited in the range 0.6-6 kPa and antiwindup integrator clamping method was adopted.

VII. CONCLUSIONS
Conducted experiments proved that the presented advanced identification method can be successfully adopted for the glass forehearth installation control system. Both original and the Strejc models give satisfying approximation of the molten glass temperature changes. Also PID controller parameters obtained based on the Strejc model were appropriate for temperature set point tracking and stabilising system in the desired operating point.