Example III: Van der Vusse Reactor

The process considered in this section is a third order exothermic van der Vusse reaction

A →B →C (3.21)

2A →D

placed in a cooled CSTR. It is a strongly nonlinear process with a nonminimum-phase behavior and input multiplicity. The model of the system is given in the following form:

where x₁ [mol/l] is the concentration of A, x₂ [mol/l] is the concentration of B, x₃ [K] is the temperature in the reactor, u₁ [1/s] is the dilution rate of the reactor, u₂ [J/s] is the heat exchanged between the CSTR and the environment, x₁₀ is the concentration of A in the inlet stream and x₃₀ is the temperature of the inlet stream. From the application point of view, u₁ is chosen as the system input whileu₂ is kept constant through the experiments [61].

To estimate the model structure, 960 data points were used (see Fig. 3.6).

In this example, the same four methods as those in the previous example were used. Likewise, 10 experiments were performed for the third and fourth methods. The input and output order were constrained to four: T = {u(k− 1),· · ·, u(k−4), y(k−1),· · ·, y(k−4)}.

In the first method, the model consisted of 45 polynomial terms (m= 8 and d = 2). Similar to the previous experience, this model was very accurate for the one-step ahead prediction, but it was unstable in the free-run simulation.

In the second method, the err’s were calculated for the 45 polynomial terms, and the terms that have very small values (below 0.01) were eliminated.

After that, seven terms remained:

3.4 Application Examples 51

0 100 200 300 400 500 600 700 800 900 1000

0 0.2 0.4 0.6 0.8 1 1.2 1.4

y(k)

0 100 200 300 400 500 600 700 800 900 1000

0 20 40 60 80 100

Time [sec]

u(k)

Fig. 3.6. Collected input–output data of the van der Vusse reactor Table 3.5. Results of Example III.

Free-run MSE^a One-step-ahead MSE^a min mean max min mean max

Method 1 - ∞ - - 0.002

-Method 2 - 20.3 - - 0.31

-Method 3 5.30 9.65 ∞ 0.088 0.24 0.62 Method 4 2.44 11.5 ∞ 0.15 0.53 0.99

aMSE in 10⁻³

y(k−2), y(k−3), u(k−1), u(k−2), y(k−4), u(k−3), y(k−1).

All of the bilinear terms were eliminated by OLS. This model is a linear model, it was stable in free-run simulation, but it was not accurate.

In contrast to the previous example, the third method result in free-run unstable models in five out of ten experiments. It is due to that the van der Vusse process has a complex nonlinear behavior, consequently this example is more difficult than the previous. The model which had the smallest MSE value in one-step ahead prediction, was unstable in free-run prediction, too.

The best model (free-run MSE) consisted of the terms

u(k−2), u(k−1), u(k−1)u(k−1)u(k−4), y(k−4), u(k−4)u(k−3), u(k−4), u(k−3).

The simplest model consisted of the terms

u(k−1), u(k−2), u(k−1)u(k−1), u(k−3)y(k−1), y(k−4).

In the fourth method, three out of the ten models were unstable. The most accurate model contained the terms

y(k−3), y(k−1)u(k−1), y(k−2)

u(k−3), (y(k−3) +u(k−2))(y(k−4) +y(k−1)), y(k−2).

The simplest model contained the terms

y(k−1), y(k−2), u(k−1)/u(k−2), u(k−2), u(k−1).

0 20 40 60 80 100

0.2 0.4 0.6 0.8 1 1.2

t

y

Fig. 3.7. Free-run simulation of resulted models. Solid line: original output, Dotted line: estimated output by Method 2., Dashed line: estimated output by Method 4.

(the best model)

As these results show, the model that gives the best free run modeling per-formance is given by the GP-OLS algorithm (see Fig. 3.7), and the structure of this model is quite reasonable compared to the original state-space model of the system.

3.5 Summary

This chapter was devoted to possible improvements in structure identification of nonlinear dynamical models. The main problem of identification of nonlin-ear dynamical models is that it is difficult to determine a model structure that is simultaneously accurate and simple. This chapter presented a new method for solving this problem for linear-in-parameters models.

One of the most preferred data structure identification method is Genetic Programming (GP). Genetic Programming selects potential solutions from a

3.5 Summary 53 given search space of possible structures. Although GP is usually able to find relatively good solutions, it tends to generate overparameterized models. This chapter proposed a method in which superfluous terms of linear-in-parameters models are eliminated during the identification procedure with the Orthogonal Least Squares method. The goal of the application of OLS for GP is to obtain more robust and interpretable models and improve the efficiency of GP.

The first part of the chapter presented the class of linear-in-parameters models, and the Orthogonal Least Squares (OLS) method. The second part of the chapter outlined the proposed approach, and the implementation of this approach in MATLAB. Finally, the proposed method was analyzed through model order identification problems of simulated chemical reactors. The re-sults of these examples justify that the GP-OLS method provides an efficient way for model order selection and structure identification of nonlinear input-output models.

4

Spline Approximation with the Use of Prior Knowledge

In many practical situations, laboratory and industrial experiments are expen-sive and time consuming and accurate measurements cannot be made. This problem results in a small number of data points that are often noisy and ob-tained at irregular time intervals. Hence, data smoothening and re-sampling are often required to handle such data sets. A good method for this purpose is the cubic spline approximation [62, 63] that provides acceptable results and has practical relevance [64], e.g. it was used for the estimation of reaction kinetic models [25, 26].

The identification of reaction rates is an important problem, as chemical process models usually contain reaction networks, and the kinetic parameters of these models often cannot be identified on the basis of a priori knowledge alone. Therefore, process modeling is generally supported by experiments to identify the kinetic parameters.

Many methods have been suggested to obtain reasonable estimates for these rate coefficients from experimental data [65]. For example, an inter-esting method for estimating rate constants of complex kinetic models from isothermal, batch or plug flow reactor data was described in the well-known paper of Himmelblau at al. [66]. This method, similarly to other approaches, is applicable when large numbers of data points along the concentration trajec-tories are available. If there are few data points, it is worth fitting individual splines for each measured variable and re-sampling the data based on these splines [25,26,64]. Unfortunately, the low number of the measured data points and the measurement noise also affect the quality of a spline approximation.

Hence, there is a need for a new approach that can handle this problem.

The main goal of this chapter is to develope a multivariate spline ap-proximation method that employs a priori knowledge in the approximation.

It will be shown that a priori knowledge, e.g. assumed material balance or the visual inspection of the data, can be easily incorporated into the spline approximation. It should be noted that, though, this idea is used alone for the cubic spline approximation in this chapter, it can be applied to other approximation techniques too.

This chapter is organized as follows. In Sect. 4.1, the applied spline model is presented. Sect. 4.2 describes howa priori knowledge can be implemented as constraints on the parameters of the splines. Sect. 4.3 presents the application examples.