Comparison of statistical metrics and a new fuzzy method for

This work compares statistical techniques for model validation to quantify the quality of identified models used in multivariable MPC controllers. Add...
2 downloads 11 Views 6MB Size
Download PDF
Article Cite This: Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

pubs.acs.org/IECR

Comparison of Statistical Metrics and a New Fuzzy Method for Validating Linear Models Used in Model Predictive Control Controllers Christiam Segundo Morales Alvarado* and Claudio Garcia† †

Laboratory of Automation and Control, Department of Telecommunications and Control Engineering, Escola Politécnica of the University of São Paulo, 05508-900, São Paulo, Brazil ABSTRACT: In the advanced process control area, model predictive control (MPC) implementations have been successful in many industrial applications. Despite being an optimization-based control technique, sometimes problems occur with the control algorithm when the dynamic model is not adequate. This work compares statistical techniques for model validation to quantify the quality of identified models used in multivariable MPC controllers. Additionally, a fuzzy validation system is proposed, showing the consistency between the model validation and the predictive controller performance. Multivariable identification, model validation, and predictive controller implementation are performed in an industrial-scale pH neutralization pilot plant.

1. INTRODUCTION System modeling is an important area of scientific research. In the last 60 years, modeling techniques have been widely used in many fields, such as aeronautics, power systems, industrial processes, and bioengineering, among others, aiming to describe the static and dynamic behaviors of complex systems, by a set of equations. The models are composed of a finite number of parameters and variables, which allow recreating the same features and behavior of the real system.1 In the process control area, dynamic models allow performing predictions and simulations, by describing process behavior in the time and frequency domains. Due to the demands of industrial processes and their degree of complexity, the aim of mathematical models is to develop control systems, such as MPC (model-based predictive controllers), by using mathematical theory of optimal control, seeking more efficiency and cost minimization, keeping the product quality. MPC controllers reduce the variability of the process, enabling to operate closer to the constraints of quality of the products, thereby maximizing the profit of plants. However, the most time-consuming task in MPC applications is to obtain the process model.2 To obtain dynamic models of industrial processes, system identification techniques are applied using the data collected under certain operating conditions of the process.3 In the last 20 years, system identification techniques have provided reliable alternatives to obtain mathematical models, applied to different chemical processes.4−7 Thus, system identification allows building model-based control systems, considering that the identified model is a reliable description of the plant behavior around the operating point of the process. © XXXX American Chemical Society

In spite of the distinct existing system identification techniques, quantifying the quality of the model is an open problem and, consequently, it is difficult to choose the most appropriate model for the control system design. The aim of this paper is (i) to develop a metric based on fuzzy logic to validate identified models, composed of statistical metrics, dynamic metrics in the time and frequency domain, and controller robustness analysis, (ii) to compare the validation results between the proposed validation index and statistical indexes used by Sprague and Geers,8 Hvala et al.,9 and Ljung,10 and (iii) to evaluate the consistency between the validation procedures and the predictive controller performance. The identified models are obtained by applying prediction error methods (PEM) proposed by Ljung,10 the asymptotic method (ASYM) proposed by Zhu,11 and techniques that optimize the predictions k-steps-ahead, such as MRI (MPC Relevant Identification), proposed by Huang and Wang12 and Potts, Romano, and Garcia.13 All the identified models are tested in the quadratic dynamic matrix control (QDMC) multivariable controller, quantifying their performance by using the PI index proposed, which is based on the process variable variability and the time response. This last is obtained when a new set point is selected. This methodology is applied to an industrial-scale pH neutralization plant, located in the Industrial Process Control Laboratory, IPCL, of the Escola Politécnica of the University of São Paulo, Brazil. Received: Revised: Accepted: Published: A

September 27, 2017 January 21, 2018 February 26, 2018 February 26, 2018 DOI: 10.1021/acs.iecr.7b04044 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

among the main statistical metrics used in the model validation are performed in subsection 2.1. In subsection 2.2, the fuzzy validation system proposed is presented. In subsection 2.3 is presented the PI index to evaluate the predictive controller performance. 2.1. Comparison of Statistical Metrics. Statistical metrics allow quantifying the agreement level between two sample vectors, i.e., real values and values predicted by the model. In addition, these metrics show the typical dispersion between the experimental result and the simulation estimated by the model. Usually, the statistical metrics used to quantify the quality of the model generated by the system identification are fit10 (eq 1), Theil’s inequality coefficient (TIC)17 (eq 2), root mean square error (RMSE)9,18 (eq 3), relative error (REL)9 (eq 4), and CS&G8 (eq 5).

This paper is organized as follows. Section 2 presents a short overview about model validation, statistical metrics, and the fuzzy validation system proposed. In section 3, a system identification procedure is applied to the pH neutralization pilot plant, using GBN signals to excite the pH process and prediction error methods, the asymptotic method, and MRI algorithms to estimate the model parameters. Section 4 describes the predictive controller used in the pH process. The experimental results of the validation procedure and set point tracking tests are shown in section 5. In section 6, the conclusions are drawn.

2. MODEL VALIDATION Validating the identified models is an important step in the system identification area, as it enables the evaluation of quality of the model according to the established objectives. The validation methods begin by comparing results of two time series (the real data collected and the simulated data of the model identified) employing visual or statistical techniques.14 The most common way to validate a model is by visual analysis, even though it is not often the best option, because it can be subjective and prone to errors. The fit index is commonly used to evaluate model performance, by comparing point-by-point the real and the simulated data. This index informs only how close the simulated data and the real data are, lacking information about how to proceed if the fit index is lower than expected. Other mathematical criteria to analyze and compare data in time series, such as statistical metrics (Theil’s inequality coefficient, root mean square error, relative error) are used, as proposed by Sprague and Geers8 and Hvala et al.9 to validate or to invalidate mathematical models. However, it is not common to relate the validation procedure and the predictive controller performance, when different identification techniques to obtain models are used. According to Sargent,15 the main goal of the model validation is to check if the model achieves its objectives, such as applicability domain and a satisfactory range of accuracy, consistent with the real system. Based on these targets, it is impossible to get a pure and absolute validation; however, there are particular established techniques that can validate models in a quantitative form, for a specific purpose and within a range of values. The model is considered valid if its accuracy is within the range of the acceptable values required for its intended purpose. A test set to determine whether the model is valid is described as follows:16 • The statistical test allows us to validate mathematical models by analyzing two data sets: the prediction data of the model identified and the real process collected data. This test requires a confidence level. • Visual comparisons intuitively allow verifying whether the predicted output of the model and the real process data are close, showing that the model reproduces the real system. For this, a set of empirical data is necessary, regardless of the data used to estimate the model (crossvalidation). • Spectral analysis allows evaluating the confidence level of the model through the frequency response, by using a particular frequency band. In spite of having other techniques to validate mathematical models, each validation test can present uncertainties, showing incompatibility among their results. Therefore, comparisons

⎛ y(t ) − y (̂ t ) ⎞ fit (%) = 100 ·⎜1− ⎟ ⎝ y(t ) − y ̅ (t ) ⎠

(1)

N

∑t = 1 (y(t ) − y (̂ t ))2

TIC =

N

∑t = 1 (y(t ))2 +

(2)

∑i (yi − yi ̂ )2

RMSE =

REL =

N

∑t = 1 (y (̂ t ))2

(3)

N

∑i ((yi − yi ̂ )2 /yi 2 ) N

(4)

CS&G = ⎛ ⎛ ⎞2 ⎛ N N 2 ∑i = 1 yi ̂ (t ) yi (t ) ⎜ ∑i = 1 yi ̂ (t ) − 1⎟ + ⎜ 1 cos−1⎜ ⎜ ⎜ ⎜ ∑N y (t )2 ⎟ ⎜ ∑N y ̂ (t )2 ∑N y (t )2 ⎜π i=1 i ⎝ ⎠ ⎝ i=1 i i=1 i ⎝

2 ⎞⎞ ⎟⎟ ⎟⎟⎟⎟ ⎠⎠

(5)

where y(t) represents the real outputs, ŷ(t) are the estimated outputs, y(t) ̅ is the mean of the real outputs, and N is the vector size. Each metric is understood as follows: • High values of the fit index indicate that the real and the predicted values are close. • For the TIC, RMSE, REL, and CS&G indexes, values close to 0 indicate that the identified model can predict values close to the real system. 2.2. Fuzzy Validation System. Herein, a fuzzy validation system with three main components is proposed: statistical tests, dynamic model evaluation, and predictive controller robustness analysis. 2.2.1. Statistical Tests. Statistical tests are used to accept or to reject the model, by imposing a confidence level. In this work, the Wilcoxon nonparametric test is used, which neither requires assumptions about the variable distribution nor is conditioned by any probability distribution and statistical parameters.19 This test is described by eq 6: Z=

T−

n(n + 1) 4

n(n + 1)(2n + 1) 24

(6)

where Z is the Wilcoxon statistical index for a confidence interval of 95%, vector sample size is denoted as n, and T represents the lower number of T+ and T−. These last B

DOI: 10.1021/acs.iecr.7b04044 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 1. Block diagram of the proposed fuzzy validation system.

parameters are the sum of the positive and negative posts of the collected data difference, respectively. 2.2.2. Model Dynamics Assessment. Model dynamics is assessed in the time and frequency domain. In both cases, the system step response of the identified models is used to quantify the validation procedure. In the time domain, the integral of time absolute error (ITAE) index is used to quantify the error e(i), represented by the difference between real and predicted data, as shown in eq 7.

σR k 2 = σξ 2 σIk 2 = σξ 2

(7)

N−1 0

(8)

k

μ I = μξ k

1 N

1 N

n=0

n=0

(13)

(14)

N

(9)

φ=

∑ |χq,α 2 i

i=1

N−1

∑ sin(Ω0kn)

(12)

n=0

Thereby, q and α represent the mathematical model parameter number and significance level, respectively. This last is considered as the probability of rejecting the null hypothesis in statistical tests when it is true. A contribution of this paper is to consider the degree of freedom as the number of parameters of the identified linear models, and as a result, in eq 15, the φ parameter allows quantifying the tolerance level between the χ2 distribution and Mnk value.

N−1

∑ cos(Ω0kn)

∑ sin 2(Ω0kn)

H1: M nk ∉ χq , α 2

In eq 8, Rk and Ik are the real and imaginary components of the error, respectively, Ω0 is the fundamental frequency, with value Ω0 = 2π/N, N is the vector size, and k represents positive integers k ϵ {0, 1, 2, ..., N − 1}. Assuming that Rk and Ik are normally distributed random variables, the means μRk and μIk and the variances σRk2 and σIk2, of the real and imaginary parts, are shown in eqs 9−12, respectively. μ R = μξ

N−1

H0: M nk ∈ χq , α 2

∑ ξ(n)e−jΩ kn n=0

(11)

n=0

The square root is extracted to obtain a greater penalty for Mnk values. In addition, Mnk values present characteristics of a χ2 distribution. To validate or to invalidate linear models, a hypothesis test is proposed, on the basis of the following expression:

In the frequency domain, the error ξk, between the real and the predicted data is analyzed using the discrete Fourier series,20 as shown in eq 8, by comparing the spectral residual with additive white Gaussian noise: ξk = R k − jIk =

∑ cos2(Ω0kn)

⎛ R k − μ R ⎞2 ⎛ Ik − μ I ⎞2 k k ⎟⎟ ⎟⎟ + ⎜⎜ ⎜⎜ ⎝ σIk ⎠ ⎝ σR k ⎠

M nk =

∑ i·|e(i)| i=1

1 N2

N−1

Then, the normalized gain Mnk is used to quantify the residual in the frequency domain, as shown in eq 13.

N

ITAE =

1 N2

− M nki|

(15)

If φ is zero, the identified model and the real process responses are close and the error presents Gaussian noise characteristics.

(10) C

DOI: 10.1021/acs.iecr.7b04044 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Low are defined as a function of the white Gaussian noise properties generated by the error between the real and the predicted data. A High variable value is above 90% when the obtained error is close to the white Gaussian noise. A Medium variable is between 20 and 60% and the Low variable is below 30% when the obtained error is far from the white Gaussian noise. (ii) For the dynamic metrics input variable (Figure 3), if

Otherwise, there are components in the frequency domain whose errors are not considered as Gaussian noise. 2.2.3. Robust Control Using Uncertain Model. In an identification procedure, there are some uncertainties in the gain or dynamics of the identified model that do not allow true representation of the real process in one operating point. These problems can be transferred to the control structure, which is not able to operate according to what was planned. Even though the uncertainties are unknown, the most common is the additive uncertainty (ΔA), as shown in eq 16. G(z) = G0(z) + ΔA

(16)

In eq 16, G(z) and G0(z) represent the real and the identified plant, respectively. In this equation, the transfer function of the additive uncertainty MA(z) is the relationship between the control signal u(z) and the uncertain signal v(z), as shown in eq 17. MA (z) =

u(z) = −K (z)(I + K (z) G0(z))−1 v (z )

(17)

Figure 3. Universe of discourse for input dynamic metrics.

where K(z) is the controller transfer function used. The small gain theorem is used to analyze the additive uncertainty of linear models, expressed by eq 18.21 ΔA

∞

= σ ̅(M(jw)A ) < 1

the model presents a good response in time and frequency domains, the linguistic variable High must be above 70%, the Medium variable is between 20% and 60%, and the Low variable is below 30%. (iii) To represent the input robustness variable (Figure 4), it can be considered that the predictive

(18)

In eq 18, ΔA is the additive uncertainty, ||·||∞ is the infinity norm, and σ̅[·] are the maximum singular values, which are used to quantify the controller robustness. 2.2.4. Proposed Method for Validating Linear Models. The block diagram of the proposed fuzzy validation system is shown in Figure 1. All values used by the fuzzy validation system must be normalized (f(u)) from 0 to 10, using eq 19: Xs = Ls + (X − Lx )

(Us − Ls) (Ux − Lx )

(19)

where X is the value obtained within interval [Lx, Ux], Lx and Ux being the minimum and the maximum values of the scale, respectively. Xs is the new normalized value within interval [Ls, Us] with the minimum Ls and maximum Us values, respectively. The components of the fuzzy system are described as follows. The fuzzification transforms scalar values to fuzzy values. For this, three linguistic variables (Low, Medium, High) and two pertinence functions (trapezoidal and triangular) are defined and adapted to the heuristic knowledge of the validation procedure. To develop the pertinence functions, we considered the following requirements: (i) for the statistical metrics input variable (Figure 2), the linguistic variables High, Medium, and

controller using the identified model can present an additive uncertainty below 1% when the identified model is close to the real process. Thus, the High variable is above 40%, the Medium variable is between 10% and 50%, and the Low variable is below 30%. (iv) For the validation value output variable (Figure 5), we considered that the High variable presents a value above 50% when the identified model is adequate to be used in the

Figure 2. Universe of discourse for input statistical metrics.

Figure 5. Universe of discourse for output validation value.

Figure 4. Universe of discourse for input robustness.

D

DOI: 10.1021/acs.iecr.7b04044 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Finally, the defuzzification transforms the fuzzy output into a real value. The centroid method is used to convert fuzzy to scalar, which determines the center of the area of the combined membership functions. The centroid method is computed by the following expression:

MPC controller, a Medium variable must be between 15 and 65%, and the Low variable must be below 35%. In the fuzzy inference machine, fuzzy rules are defined to take a decision. This stage is computed using an inference engine, knowledge base, and membership functions. The mathematical translation of each fuzzy rule is converted into an output value of each fuzzy input. We used here the fuzzy rules proposed by Mamdani,22 known as the Max-Min method. Twenty-seven rules were implemented in the fuzzy algorithm using IF-THEN structures. Figures 6, 7, and 8 show the surface

n

G (B ) =

∑i = 1 uiφB(ui) n

∑i = 1 φB(ui)

(20)

where G(B) is the defuzzified value (centroid), n is the number of elements, and φB(ui) represents the membership value of element ui in fuzzy set B. 2.3. Proposed Index To Evaluate the Performance of Control Loops. To quantitatively assess each model in the QDMC controller, two performance indexes are used: rise time (tr) and variability index (VI), proposed by Shunta.23 Index tr is the mean time needed for the process variable to cross the new set point for the first time. The variability index is calculated as seen in eq 21. ⎛ σ + s⎞ VI (%) = 100⎜1 − fbc ⎟ σtot + s ⎠ ⎝

(21)

In eq 21, σfbc and σtot are the theoretical minimum deviation and the total standard deviation of the output variable, respectively. The s parameter is the sensitivity factor of the variable scale, with a value of 0.1. With these two indexes, a performance index PI (%) was proposed to evaluate the QDMC controller using all the identified models, as shown in eq 22.

Figure 6. Surface generated by the rules of the inference machine between robustness and statistical metrics inputs.

PI (%) = α(100 − VI) + βtr

(22)

α and β are weighting factors equal to 0.7 and 0.3, respectively.

3. BLACK BOX SYSTEM IDENTIFICATION OF THE PH NEUTRALIZATION PROCESS pH neutralization can be found in different industrial processes, such as water treatment and biotechnological and chemical processes. In these processes, the pH is recognized as highly nonlinear, difficult to control, and very sensitive when the pH value is close to the neutral point. Much research has been related to pH control.24,25 These advanced control techniques allow us to solve the complex problem of pH control. One of the most used techniques is the model-based predictive control (MPC). To develop MPC controllers, the process model must be adequately obtained, through system identification techniques. System identification allows building mathematical models of a dynamic system, based on measured data and previous process knowledge. To create a process model, it is necessary to choose an appropriate structure with a finite number of parameters.10 In this application, the identification process was performed in an industrial-scale pH neutralization plant in closed loop, considering that the level and the pH PID controllers were kept in automatic mode, because it is not recommended to open the pH loop, due to the process high gain around the neutral pH. Because the pH is very sensitive to minimum variations in the flows of acid or base around its neutral value, if the process is left in an open loop, the pH of the solution tends to become acidic or basic over time. The pilot plant was built using industrial instrumentation such as sensors, control valves, a programmable logic controller

Figure 7. Surface generated by the rules of the inference machine between robustness and statistical dynamic metrics.

Figure 8. Surface generated by the rules of the inference machine between statistical metrics and dynamic metrics inputs.

generated by the fuzzy rules formulation, showing the relation between the inputs “robustness” and “statistical metrics”, “robustness” and “dynamic metrics”, and “statistical metrics” and “dynamic metrics”, respectively, with “validation value” output. E

DOI: 10.1021/acs.iecr.7b04044 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 9. P&ID of the pH neutralization plant.

the process. Therefore, random binary excitation signals known as GBN (generalized binary noise26) are used in the closed loop identification. This binary signal switches between two values −A and +A, with nonswitching probability pn‑sw varying between 0 and 1. The latter value is computed by using the dominant time constant τdom and the frequency ω of the GBN signal. In the pH neutralization plant, the dominant time constant of the level loop τlevel and the dominant time constant of the pH loop τpH are 124 and 83.5 s, respectively. The sampling times for the level and the pH loops are calculated on the basis of expressions 23 and 24, respectively. τ Δtlevel < level → Δtlevel ≅ 10 s (23) 10

(PLC), and a distributed control system (DCS), commonly used in industrial processes. According to Figure 9, the piping and instrumentation diagram (P&ID) of the pH neutralization plant is composed of the following main parts: • A primary acid tank (PAT), where the hydrochloric acid solution (HCl) is prepared with a concentration of 0.0056 mol/L and then stored in an intermediate tank (PATI) with controlled level, to maintain a constant flow entering the reactor; • A tank where the basic solution (BT) of sodium hydroxide (NaOH) is stored with a concentration of 0.0185 mol/L, used to neutralize the pH in the reactor tank. The flow to the reactor is manipulated by a dosing pump, and • A reaction tank (RT), where the pH neutralization occurs. To homogenize the mixture of reactants, a mechanical stirrer is used. This reactor has a pH sensor (AE-40) and a level transmitter (LIT-10). The level is controlled by a solenoid valve that manipulates the outlet flow. This pilot plant has two measured variables in the tank reactor: level and pH. The reactor level is controlled by a PWM (pulse width modulation) signal, which has a variation of 0− 100%, equivalent to 0−10 s, keeping the solenoid valve open or closed. For the pH control loop, a dosing pump is used to control the pH, varying its pulsation frequency. This signal has a variation of 0−100%, generating a pH variation between 2.35 and 12.21, approximately. 3.1. Excitation Signals Design for Multivariable Identification. The design of the excitation signals allows gathering data for the identification. The data must contain information regarding the time and the frequency domains of

Δt pH