On the Structural Optimization of a Neural Network Model Predictive

Article pubs.acs.org/IECR

On the Structural Optimization of a Neural Network Model Predictive Controller Mahsa Sadeghassadi,* Alireza Fatehi, Ali Khaki Sedigh, and SeyedMehrdad Hosseini Advanced Process Automation and Control (APAC) Research Group, Industrial Control Center of Excellence, Faculty of Electrical & Computer Engineering, K. N. Toosi University of Technology, P.O. Box 16315-1355 Tehran, Iran ABSTRACT: This paper provides a new algorithm for tuning the two most effective parameters in nonlinear model predictive control (NMPC). Tuning is performed in two steps. First, a new method based on Barron’s formula, the bicoherence nonlinearity test, and inphase-quadrature demodulation is proposed to determine the number of hidden layer neurons in a two-layer neural network. In the second step, a fuzzy algorithm is introduced to tune the input weight matrix in the objective function to make the tuning problem more practical and precise. To show the effectiveness of the proposed method, several examples are discussed including a simple flow process and a more complex pH neutralization problem. The method is also evaluated in the laboratory scale pressure and level processes. It is shown that the proposed method leads to tuning the number of neurons and the weight matrix with an acceptable performance. strategies available in the literature. Shridher et al.8 proposed a tuning strategy by employing the first order plus delay time (FOPDT) models. Neshasteriz et al.9 proposed a predictive control strategy for the plants described by second order plus delay time (SOPDT) models and an automatic method for selection of input weight factor to ensure stability and performance. Iglesias et al.10 used the analysis of variance (ANOVA) on the parameters of an estimated FOPDT model to determine the input weighting vector as a function of these approximated parameters. In all of the above methods, a linear model is used for tuning the parameters. Recently, some researchers have paid attention to tuning methods which directly use a nonlinear model of the process to tune the MPC parameters. For example, an off-line tuning method is presented by Ali et al.11 in which a multiobjective optimization problem is solved for a nonlinear process model. An online fuzzy based mechanism is also presented by Ali et al.12 which uses time domain performance specifications for developing fuzzy rules. Ghazzawi et al.13 presented an online adaptation strategy by deploying the sensitivity analysis of the closed loop response to the tuning parameters. Kawai et al.14 used the PSO algorithm to tune NMPC parameters. This paper deals with the above two issues in the NMPC. In the first step, a new algorithm based on Barron’s formula is presented to determine the number of hidden layer neurons. In the next step, a fuzzy inference system is used to determine the input weighting factor in the NMPC objective function. The paper is organized as follows: In section II, some mathematical backgrounds are discussed. In section III, the method for determining the number of hidden layer neurons using Barron’s method is described. Tuning of MPC objective function parameters is proposed in section IV. In section V, simulation results

I. INTRODUCTION Model predictive control (MPC) is a model based controller which determines the control signal by optimizing a predefined objective function. The objective function includes the control signal and the process model output deviation from the reference trajectory over a future time horizon. Thus, a suitable model is used to predict the future behavior of the process, and tuning the weights of the prediction error and the control signal affects the overall controller performance. Linear prediction models are widely used due to their online optimization simplicity, as reviewed in Camacho et al.1 However, in the presence of process inherent nonlinearities, linear prediction models may fail and in such cases a nonlinear model in MPC formulation becomes necessary. Different nonlinear models have been developed in the literature of nonlinear model predictive control; see, for example, Allgower et al.2 Among them, neural networks are wellknown and popular due to their global approximation property. However, selecting an optimal structure of neural networks for a specific problem is an unsolved task. For a two-layer feedforward neural network, an important parameter is the number of neurons in the hidden layer. Exhaustive search is favored among the neural networks by practitioners; however, it faces the combinatorial complexity problem. In this method, neural models of different complexity are estimated and the resulting errors on the test data set are compared.3 This method is computationally demanding due to the complexity in the estimation of each nonlinear model. Wu et al.4 presented some empirical rules to obtain the upper and lower bounds on the number of hidden neurons. Barron5 proposed a method of optimizing the number of hidden layer neurons using complexity regularization and statistical risk. Random search optimization approaches can also be used to select the hidden neuron numbers in order to minimize a predefined objective function. For example, Suraweera et al.6 used the particle swarm optimization (PSO) algorithm to search for the optimum neuron numbers and Son et al.7 used the genetic algorithm (GA). Weighting parameters of the objective function also strongly affect the controller performance. There are specific tuning © 2012 American Chemical Society

Received: Revised: Accepted: Published: 394

August 15, 2011 October 8, 2012 November 26, 2012 November 26, 2012 dx.doi.org/10.1021/ie201815r | Ind. Eng. Chem. Res. 2013, 52, 394−407

Industrial & Engineering Chemistry Research

Article

Table 1. Fuzzy Rules for Computation of γ

Figure 1. Fuzzy sets of ISE on a normalized universe of discourse.

Figure 2. Fuzzy sets of oscillation measure.

Figure 4. pH Process: Estimation of titration curve.

where the notation ∑Π means summing over all distinct ways of partitioning x(k + τ1)x(k + τ2), ..., x(k + τ2N) into pairs. Signals such as sine wave and pseudo random binary signal (PRBS) are in accord with this class of inputs.16 If a Gaussian input passes through a nonlinear block, this theorem can be used as a nonlinearity measure, but for non-Gaussian input, there is no analogous theorem.17 This is a key assumption for bicoherence nonlinearity measure which is mentioned by Choudhury et al.18 and is used in this paper. The cumulant-generating function Cx(t) is defined as

Cx(t ) ≡ ln(Mx(t ))

(3) tx

where Mx(t) is the expectation of the transformation e , which is called the moment-generating function. Cumulants can be derived from the Taylor’s series expansion of the cumulant generating function. The nth order cumulant can be written for any random series x(k) as

Figure 3. Fuzzy sets of γ on a normalized universe of discourse.

are presented. Experimental results are proposed in section VI. Finally, section VII provides the conclusion and summaries.

cn(τ1 , τ2 , ..., τn − 1)

II. MATHEMATICAL BACKGROUND a. Nonlinearity Measure. There are several measures to estimate the nonlinearity extent of a given data series x(k). Many of the nonlinearity measures can be performed if the input signal x(k) is jointly Gaussian random, Isserlis;15 i.e., the input signal series x(k + τ1)x(k + τ2), ..., x(k + τ2N+1) (N = 1, 2, ...) are normalized, which means E[x(k + τi)] = 0, E[x2(k + τi)] = 1, and also E[x(k + τ1)x(k + τ2), ..., x(k + τ2N + 1)] = 0

≡ cum[x(k)x(k + τ1)x(k + τ2), ..., x(k + τn − 1)] =

∑ ∏ E[x(k + τi)x(k + τj)]

t=τ

(4)

Bicoherence is formally defined as the double discrete Fourier transform (DDFT) of the third order cumulant of the time series. This becomes

(1)

B(f1 , f2 ) ≡ DDFT[c3(τ1 , τ2)] ≡ E[X (f1 )X (f2 )X *(f1 + f2 )]

(2)

where X(f) is the discrete Fourier transform coefficient of the time series x(k) at frequency f. The asterisk is the complex conjugate operator, and E[·] is the expected value operator.

(5)

E[x(k + τ1)x(k + τ2), ..., x(k + τ2N )] =

∂ nCx(t ) ∂t n

395

dx.doi.org/10.1021/ie201815r | Ind. Eng. Chem. Res. 2013, 52, 394−407


Article

Figure 5. pH Process: Nonlinearity indices for 10 points among the whole search area: (a) in first iteration; (b) in last iteration. Dashed lines are related to the γ determination stage.

Note that, as the second order cumulant is a measure of the correlation between the two frequencies, bicoherence is a measure to show the skewness or the extent of phase interaction in a signal. As a result, it can be used to quantify the certain type of nonlinearity according to the quadratic phase coupling problem. Quadratic phase coupling is the presence of additional harmonics (f1 + f 2) in the filtered signal which are phase coupled to the original frequencies f1, f 2. By normalizing the bispectrum over the range [0 1], bicoherence is defined as bic 2(f1 , f2 ) =

|B(f1 , f2 )|2 E[|X(f1 )X(f2 )|2 ]E[|X(f1 + f2 )|2 ]

Figure 6. pH Process: Spectral representation vs frequency and operating condition.

(6)

According to Choudhury et al.,18 a nonlinearity index (NLI) is defined as NLI = |bic max 2 − (bic 2 + 2σbic2)|

(7)

2

where bicmax is the maximum value of the squared bicoherence. bic 2 and σbic2 are the average and standard deviation of the squared bicoherence, respectively. This measure determines whether the system is linear or nonlinear and also the nonlinearity extent of the system. Theoretically, for linear signals, squared bicoherence should be a constant for any bifrequencies and bicoherence should be flat in a three-dimensional space. As a result, NLI should be zero for linear systems. However, practically there are confidence levels around zero where NLI can be assumed zero and the process can be assumed linear. For example, for a data length of 4096, if NLI is 0.01 or less, the system can be assumed linear. It is preferable to use a large number of data points for the bicoherence nonlinearity measure. However, if it is not possible to obtain a long data set, the threshold values used in this measure can be changed to obtain reliable results. Choudhury et al.18 proposed some threshold values according to the data length. The threshold values can also be obtained experimentally. In this paper, the plant is excited by sinusoidal inputs, which meet the Isserlis condition,15 to obtain the required data for the bicoherence test. The oscillation period of input signal is userdefined and chosen around the time constant of the process based on the experience of authors. To use a threshold value of 0.01, it is important to meet the assumptions underlying the suitable data length of 4096 to be used in the nonlinearity test. Thus, the sampling time of the nonlinearity test is considered small enough to have the data length of 4096 in a reasonable

Figure 7. pH Process: Determining the optimal number of neurons using the trial and error method.

running time. In this paper, the sampling time is selected to have 128 samples in each oscillation period of the input signal. Therefore, the sufficient data length is obtained through 32 oscillation cycles. b. Describing Function Analysis. Assume that a stable system is excited by an input signal of the form u(t) = u0 + a sin(ω0t + φ0) which is a sinusoidal plus a constant u0, that is, the input operating point. If the system is linear, the output will be a sine wave of the same frequency but with a scaled amplitude and shifted phase. In the presence of nonlinearity, however, the output includes harmonics of the input signal frequency. These virtual harmonics show the nonlinear behavior of the system. The frequency response can be used to characterize linear systems. However, many of the nonlinear systems show quasi-linear behavior in the sense that the amplitude reduces 396



Article

Figure 8. pH Process: (a) Output training data, (b) input training data, (c) output test data and one step ahead prediction, and (d) prediction error.

quickly with each successive harmonic in the spectrum that represents the low-pass characteristic of the system. Such systems can be described by an approximately equivalent transfer function that is the ratio of the phasor representing the fundamental component of the output to the phasor representing a sinusoidal input signal. Because of the low-pass characteristic of the systems studied in this research, this condition is satisfied for them. There are many demodulation methods that can be used to separate the fundamental frequency. In this paper, the inphasequadrature (IQ) demodulation method is used for frequency separation: I1(ω) =

For target functions satisfying a given smooth property, the mean integrated squared error between the estimated network and the target function is bounded by O(Cf 2/n) + O(nd log N /N )

where d and n are the size of input and hidden layers, respectively. N is the size of training data. O(·) plays the balancer role between accuracy and neural network size. Cf is defined as the Fourier magnitude distribution of the target function f. The Fourier representation of function f is

Tss + k × T

∫T +(k−1)×T y(t ) exp(−j(ωt + θ)) dt ss

f (x ) = (8)

Cf =

eiωxF (̃ ω) dω

(12)

∫

|ω||F(̃ ω)| dω

(13)

where |·| stands for the Euclidean norm. The main problem on deploying this function is that the target function f is unknown in most practical cases. Therefore, Cf should be derived from input−output data of the process. Wechsler21 showed that Cf relates the system frequency domain properties to the network complexity. Thus, it can be estimated by analyzing the spectral characteristics of the target function. b. Neuron Number Determination. In this paper, a new approach is presented to determine Cf from the input−output data. Since the process is nonlinear, describing function approximation could be used to estimate the pseudo frequency response at the frequency vector ω around a certain operating point. However, each of these quasi-linear frequency response models are only valid around a certain operating point and have a different value for Cf. The greatest value of Cf which is for the operating point with the highest nonlinearity extent can be used to determine the neuron numbers. For this purpose,

(9)

This procedure should be repeated over a range of frequencies to obtain the sinusoidal input describing function model of the nonlinear plant.

III. PREDICTION MODEL STRUCTURE a. Problem Formulation. Barron5 proposed a method for approximating the number of hidden layer neurons in a two-layer feedforward neural network with a sigmoidal activation function. In Barron’s formula, neuron numbers are determined by minimizing a specific objective function. It is shown by using a neural network model with the following number of neurons: n ∼ Cf (N /(d log N ))1/2

∫

If ωF̃ (ω) is integrable, Cf can be defined as

where Tss and k are the steady state time and period index, respectively. y(t) is the system response to a sinusoidal with frequency ω plus a constant which drives the system to the operating point. According to Nassirharand et al.19 and Schwartz et al.,20 pseudo frequency response is calculated as follows: G1(ω) = 2I1(ω)/(aT )

(11)

(10) 397



Article

Figure 9. pH Process: Parameter convergence of γ and response obtained with NMPC in (a) region 1, (b) region 2, (c) region 3, and (d) region 4.

a bicoherence nonlinearity test is performed to find the operating point with the highest nonlinearity extent. In this point, the plant is more complex and, therefore, the model should be more flexible. For a neural network model with one hidden layer, it means that the number of neurons in the hidden layer has to be increased. Thus, finding the number of neurons for that operating point determines an upper bound which guarantees the sufficient neuron numbers for other operating points.

The procedure for determining the number of hidden layer neurons and training the network involves the following steps: (i) Specifying an approximation of the process time constant. (ii) Performing a bicoherence test to find the maximum NLI over the system operating range. (iii) Using the inphase-quadrature demodulation method to obtain G1(ω) in eq 9. (iv) Replacing F̃ (ω) by G1(ω) to find Cf in eq 13. 398



Article

Table 2. Estimation of the Closed Loop Response for Initial and Final Values of γ operating point

initial value of γ

final value of γ

original overshoot

trained overshoot

overshoot of reference model

initial ISE

final ISE

pH ∼ 6.5 pH ∼ 7.8 pH ∼ 8.3 pH ∼ 9.3

γ = 300 γ = 300 γ = 300 γ = 300

γ = 18.34 γ = 42.24 γ = 16.53 γ = 3.12

17.44% 18.02% 11.52% 6.66%

27.23% 25.58% 24.97% 25.44%

25.58% 25.58% 25.58% 25.58%

0.014 0.012 0.022 0.036

0.0019 0.0001 0.0001 0.0097

(v) Specifying the training data for identification which has the proper amplitude and frequency characteristics. (vi) Calculating the neuron numbers according to Barron’s formula (eq 10). It is important to note that three types of data series are used in the above algorithm. The first type is the data which is used for the bicoherence test. To perform this test, a sine wave is used as the input signal according to section II.a. Note that the absolute value of the nonlinearity index for an operating point is not as important as its relative value compared to other operating points. Thus, a comparison can be done for every data length (greater than the DFT length) if we use the same data length for all of the operating points. In this paper, a data length of 4096 is used for all of the simulations and implementations. However, the length of the data could be reduced by using a higher threshold or different types of nonlinearity measures. The second type is the signal which is used for performing inphase-quadrature demodulation and calculating parameter Cf. To obtain the necessary data, the system is excited in the most nonlinear operating point by a sine wave according to section II.b. The length of the signal is selected so that a complete period of output is accessible after the transient response. Thus, the length of data series for calculating Cf is in the order of 2 times the settling time of the plant. The third type is the training signal. An appropriate excitation signal for training the neural network nonlinear model is designed on the basis of Nelles.3 Neural networks are used to model processes with some extent of nonlinearity. Thus, it is necessary to provide a training signal which is sufficiently rich to identify a nonlinear process. An amplitude modulated pseudo random binary signal (APRBS) is used for this purpose. This signal can be designed to capture the limits of the process input range between umin and umax. Such a signal can drive the process through all the operating regions. A generally acceptable length of training data in a linear identification is around 500−1000 data points if a PRBS signal is used.22 In nonlinear identification, data needs to be gathered in different operating points. Therefore, the length of training signal N in eq 10 is chosen according to the operating area selected to train. The number of data points can be smaller if a narrower operating area is selected to train. The frequency range of the input signal is selected on the basis of the modeling purpose and limitations for the length of the training data. In this paper, the model is used in a control framework, so it is necessary to capture frequencies around the crossover frequency.3 Therefore, the frequency of input signal is selected in the range of [1/nτ n/τ ], where τ is the time constant of the process. n is chosen according to the frequency range suitable for identification. It is based on the experience of the user and the process properties. In this paper, n is chosen from the set of [3, 4, 5] and depends on the system. The sampling time for identification is also different with the sampling time of the bicoherence test. The sampling time for identification is selected as 1/m of the process time constant.

Figure 10. pH Process: Performance of nonlinear MPC: (a) set point and process output; (b) control signal.

m is chosen in the range of [5 15]. The sampling time for the bicoherence test is chosen according to the last paragraph of section II.a. Thus, it is much lower than the sampling time that is used for identification.

IV. OBJECTIVE FUNCTION TUNING The tuning parameters in the model predictive controller are the prediction horizon, control horizon, reference trajectory bandwidth, and diagonal matrices containing input and output coefficient weights. Prediction and control horizons are set to reasonable predefined values with respect to process time delay and computational limitations.8 The output weight matrix Λ is often the identity matrix. In most cases, the input weight matrix is set to Γ = γI, as I is the identity matrix. Our goal is to find the parameter γ in the NMPC. To determine γ, we follow its effect on the closed loop performance. The desired closed loop performance is denoted by a reference model. The reference model is chosen on the basis of criteria such as rising time, settling time, or overshoot. A fuzzy inference system (FIS) is deployed to determine γ. The FIS tunes γ so that the integral of squared error (ISE) between the closed loop output and the desired reference and the extent of oscillation in process output are reduced. Input and output fuzzy sets include seven membership functions (MFs). In this paper, all the membership functions are in the shape of symmetric triangles, except the two MFs at the extreme ends. Figure 1 depicts the membership functions for the first input. This input of the fuzzy inference system is scaled by a coefficient, gISE, in order to match the range [−1 1]. Selection of a suitable value for this scaling parameter is made on the basis of the knowledge about the process to be controlled and sometimes through trial and error to achieve the best possible one. 399



Article

The following measure is used for the oscillation detection:23 osccount = |∑ e| /∑ |e|

The oscillation measure is bounded between 0 and 1. The system responses with osccount = 1 have no oscillation, and systems with osccount ≤ 0.4 start a significant oscillation. Then, seven triangular membership functions are defined over the range [0 1], as illustrated in Figure 2. We also use seven triangular membership functions of Figure 3 for the output which is scaled by gγ. Selection of the suitable value for gγ is made to achieve the best possible control performance and sometimes based on the actuator limitations. The effect of γ on the stability and performance of the closed loop system should be studied to determine the rule base of the FIS. γ has to be decreased when the convergence of the process output to the reference trajectory needs to be sped up and increased when the output starts to oscillate. On the basis of this idea, γ must be decreased to improve the closed loop speed, before a significant oscillation starts in the process output. Table 1 illustrates the rule base used in the FIS. Mamdanie’s minimum implication and center of gravity defuzzifier are used in this paper.

(14)

Figure 11. Flow process: Nonlinearity indices for nine points among the whole search area.

V. SIMULATION RESULT To show the effectiveness of the proposed algorithm, two process simulated examples are used. These are the pH neutralization process as a complex plant and the flow control process as a simple plant. Example 1. The first example is a pH neutralization process which consists of a reactor tank where the pH of the content is controlled by sodium hydroxide flow rate as the manipulated variable and an unknown value for acetic acid (CH3COOH) rate which can be considered as the disturbance. Also, water flow is used to control the tank level in h = 10 cm. These streams are injected into the tank with dosing pumps. The objective is to produce a certain pH value which ranges from pH ∼ 6 to pH ∼ 9.5. The model of this plant is completely explained in Appendix A. pH control is a difficult task due to process high

Figure 12. Flow process: Determining the optimal number of neurons using the trial and error method.

Figure 13. Flow process: (a) Output training data, (b) input training data, (c) output test data and one step ahead prediction, and (d) prediction error. 400



Article

Figure 14. Flow process: (a) Parameter convergence of γ, (b) reference trajectory and the response obtained with NMPC using the initial and tuned parameter of γ. (c) Control signal obtained with NMPC using the initial and tuned parameter γ.

nonlinearity. The nonlinearity can be described by a titration curve, as shown in Figure 4. This curve shows the static nonlinearity which is the sensitivity of the steady state response to changes in the base flow rate. The slope of output change is very steep around pH 7.8. The steady state gain of process changes up to 20 times near this point compare to that around pH 6. The proposed method is now used to find the operating point with the highest nonlinearity extent. The process time constant is about 460 s. The input signal and the sampling time for bicoherence test are chosen to meet the conditions in section II.a. Thus, the sampling time for the bicoherence test is considered about 3.6 s. The input operating value u0 uniformly changes from 0.73 to 0.78 for 10 points. As a result, pH changes from pH 6.25 to pH 9.3. The nonlinearity indices in these points are shown in Figure 5a. The maximum value of NLI is between pH 8 and pH 9. However, the value of index changes rapidly in this operating range. Thus, 10 new points are selected in this operating range. The procedure continues until the difference between the maximum NLI and its neighbors is less than 0.01 which is the linear threshold in the bicoherence nonlinearity test. The nonlinearity indices for the last 10 points are shown in Figure 5b. The maximum NLI is for pH ∼ 8.315. Now we want to compare the result of the nonlinearity test with the titration curve by using spectral analysis, as shown in Figure 6. Figure 5a shows that the maximum values of nonlinearity indices are related to pH ∼ 7.2 and pH ∼ 8.3, which also have the maximum height of the second harmonic. Taking another look at the spectrum, we can see that the line with maximum height in the third harmonic is related to pH ∼ 7.8, which is the point of inflection of the titration curve. At this point, the value of the second harmonic is equal to zero, which coincides with the zero value of the nonlinearity index. By using the proposed strategy in section III, Fourier magnitude distribution is evaluated as Cf = 0.4973. Then, the

Figure 15. Flow process: Performance of nonlinear MPC: (a) set point and process output; (b) control signal.

number of neurons based on Barron’s formula for a data length of 23 480 and d = 5 as the input layer size is n = 16. To justify our proposed method, the result is compared with the best number of neurons when an exhaustive search method is used. Figure 7 shows that the result of our proposed method is sufficiently accurate. Although the error is reduced by increasing the number of neurons, it is important to note how much the fit is improved for the higher neuron number. Our aim is to avoid unnecessarily complex models. Choosing n = 16 as the number of neurons is a reasonable choice according to both fitness qualification and model complexity. 401



Article

about 30 s, which is chosen regarding the time constant of the process that is about 460 s. Thus, the sampling time is about 8 times larger than the sampling time used in the nonlinearity test. A long data set is required to excite the process through the desired operating region. The number of data points can be smaller if a narrower operating area is selected to train the network. The quality of modeling is shown in Figure 8c and 8d when validation data is introduced to the plant. The next step is to determine γ. Since the process behavior is much different in various operating points, it is necessary to choose different values of γ for these operating points. Therefore, the NLI is deployed to divide the operating points into four different regions which were shown in Figure 5 with dashed lines. The parameter γ changes from region to region but remains constant in each of them. The same reference model applies for all of these small regions. The response of the reference model is not necessarily better than the response of the process, but our aim is to make the system response as close as possible to the reference regardless of the fact that it becomes better or worse. The initial value for γ is set to 300 to be large enough to avoid instability during the tuning procedure. Figure 9 shows how this tuning parameter converges to different values in various regions in order to track the same desired closed loop reference model. It can be seen that the responses for γ = 300 are all sluggish and cannot track the reference trajectory. The closed loop responses for γ's at their final values are shown in the same figures and are summarized in Table 2. The specifications of the system response when gamma is at its initial value are named as original overshoot and initial ISE, while these parameters are named trained overshoot and final ISE when gamma converges to its final value. It is clear from Table 2 that the response specifications are very close to model reference specifications for the final value of γ. Thus, we can see that the algorithm is successful in generating suitable values for γ and the system outputs can accurately track the

In what follows, comparisons of results are provided. Kocijan et al.24 used the neural network model with 10 neurons in the hidden layer to capture the dynamics of the pH neutralization process. Yu et al.25 used eight hidden neurons to provide an appropriate model of the pH subsystem for a nonlinear predictive controller. Wior et al.26 found a nonlinear autoregressive exogenous (NARX) model for the experimental neutralization plant. The result was a neural network which has 11 neurons in the hidden layer. Note that our proposed method determines the number of neurons for the most nonlinear operating point. Thus, it is obvious that the number of neurons determined by this algorithm is slightly more than the number of neurons in the above references. The next step is the estimation of a NARX model by training a two-layer MLP neural network with n = 16 sigmoid neurons in the hidden layer. The training procedure with the Levenberg− Marquardt adaptation method is implemented by applying the APRBS input data shown in Figure 8b and the process output shown in Figure 8a. Note that the identification sampling time is

Figure 16. Pressure pilot plant: Nonlinearity indices for nine points among the whole search area.

Figure 17. Pressure pilot plant: (a) Output training data, (b) input training data, (c) output test data and one step ahead prediction, and (d) prediction error. 402



Article

Figure 18. Pressure pilot plant: (a) Parameter convergence of γ, (b) reference trajectory and the response obtained with NMPC using the initial and tuned parameter of γ. (c) Control signal obtained with NMPC using the initial and tuned parameter of γ.

Figure 19. Pressure pilot plant: Performance of nonlinear MPC around different operating points.

Example 2. This example considers the model of a flow pilot given in Maghoul.27 This pilot is used to maintain a constant predetermined water flow in a collection of pipes. The plant main parts are a source tank, a control valve, a pump, and some pipes. The operating region is limited between ymin ∼ 60 dcm3/h and ymax ∼ 140 dcm3/h. The nonlinearity test is in accord with the conditions mentioned in section II.a. The process time constant

reference model. As a result, the value of γ increases as the steady state gain increases, which coincides with what is expected experimentally. Figure 10 shows the results of applying the final NMPC controller to the process in different operating points. As depicted from this figure, the behavior of the closed loop control is more or less the same in all the operating regions. 403



Article

tuning algorithm tried to speed up the process response by reducing the weighting parameter which is initially chosen too large. The tuning procedure is stopped after a finite number of iterations at a value which gives the best trade-off between the error correction and the oscillation limitations. In order to use this value of γ in the NMPC structure, it is necessary to test the quality of process response all over its operating condition. The output performance in eight operating points is shown in Figure 15. This figure shows the acceptable tracking quality in all of the set point changes.

is about 50 s. Thus, the sampling time for nonlinearity test is considered about 0.4 s. Figure 11 is obtained by using nine points in this operating range and calculating the bicoherence measure and nonlinearity indices. NLI values over plant operating conditions are less than 0.01. Thus, the process can be considered linear. We get Cf = 0.0257. For N = 38027 and d = 5, the number of hidden layer neurons is n = 1 which emphasizes that the process is more linear. Note that, if n is obtained as less than two, it will be considered as n = 2. Figure 12 shows the result of neuron number determination using the exhaustive search method. It is clear that n = 2 is sufficiently appropriate. Figure 13 illustrates the input and output training series. As mentioned, the identification sampling time is selected in view of the process time constant that is about 50 s. Thus, the identification sampling time is 7 s. Estimation error in Figure 13d shows that the prediction model is sufficiently accurate. The tuning algorithm modified the weighting parameter in order to reduce the tracking error at the most complex operating point which is at 60 dcm3/h. It is obvious from Figure 14 that the

VI. EXPERIMENTAL RESULT The proposed method is implemented in laboratory scale pilot plants. The first step is to select the model structure. The idea is to select the regressors based on insights from the linear system identification and then determine the best possible network structure with the given regressors as inputs. Thus, the output order, input order, and delay time are set to 3, 2, and 1, respectively, and the NARX architecture is used as the nonlinear model in the following plants. Case 1. The first plant is a pressure pilot.28 The system main parts are a control valve, a compressor, and a tank. The compressor takes in the air at atmospheric pressure, raises the pressure of air up to 400 KPa, and then lets the compressed air through pipes into the plant. A valve is used to regulate the tank pressure. The output of the system is the difference between the absolute pressure of the tank and the atmospheric pressure. Hence, the output operating range is from ymin ∼ 0 KPa up to ymax ∼ 200 KPa. The bicoherence nonlinearity measure is performed to determine the most nonlinear operating point of the plant. The sampling time for bicoherence test is chosen according to section II.a. Thus, a data length of 4096 is gathered in each nine operating points uniformly spanned in the operating range. Results of the nonlinearity test are shown in Figure 16. It is obvious from Figure 16 that the NLI value increases as the

Figure 20. Level pilot plant: Nonlinearity indices for six points among the whole search area.

Figure 21. Level pilot plant: (a) Output training data, (b) input training data, (c) output test data and one step ahead prediction, and (d) prediction error. 404



Article

Figure 22. Level pilot plant: (a) Parameter convergence of γ, (b) reference trajectory and the response obtained with NMPC using the initial and tuned parameter γ. (c) Control signal obtained with NMPC using the initial and tuned parameter γ.

Figure 23. Level pilot plant: Performance of nonlinear MPC around different operating points.

The nonlinear model is generated by using the training series shown in Figure 17a and 17b. The accuracy of this model is tested by using a series of test input−output data, as shown in Figure 17c and 17d). The sampling time for identification is 1 s, so this data length has been obtained in less than 3 h. Since the NLI does not drastically change in the operating regions, γ is tuned for the operating point y = 191 KPa, which has

pressure value increases. It also shows that the system has a nonlinear behavior in high pressure as the nonlinearity indices change from 0.04 in low pressure to more than 0.15 in high pressure. By using the proposed strategy in section III, the Fourier magnitude distribution is evaluated as Cf = 0.0998. Then, the efficient neuron number for a data length of 7826 and d = 5 is determined as n ≈ 3. 405



Article

the maximum value of the NLI. The initial value of γ is set to 400; however, after the training, it is tuned to 1313.7. Figure 18 shows the result of tuning γ together with its implementation in the NMPC. Figure 19 shows that this value of γ is also suitable in generating an appropriate response in various operating points. Case 2. The second plant is a laboratory size level pilot.29 The process consists of a cylindrical water tank that is approximately 60 cm high, a pump, a control valve, a source tank, and several water pipes. The control valve is used to regulate the level of tank. The water level is considered in the range of [10 60 ] cm. The result of the bicoherence test is shown in Figure 20. The maximum NLI is less than 0.1, so the process is not highly nonlinear. It is obvious from Figure 20 that the NLI does not drastically change in the operating region. However, the proposed algorithm is used for the operating point of 32, which has the maximum NLI. We get Cf = 0.0207. Then, the efficient neuron number is obtained as n = 2 for a data length of 7826 and d = 5 as the network input dimension. The identification result with the NARX model and n = 2 is shown in Figure 21. The result of tuning γ is shown in Figure 22. Figure 23 shows that, by using γ = 745.8, the tracking procedure is done successfully.

Kb =

K w = [H+][OH−]

x1 = [CH3COOH] + [CH3COO−]

(A.5)

x 2 = [Na +]

(A.6)

[H+] + [Na +] = [OH−] + [CH3COO−]

(A.7)

Equation A.7 can be rearranged to give a cubic polynomial in [H+]. The resulting equation is [H+]3 + [H+]2 (ka + x 2) + [H+](ka(x 2 − x1) − k w ) − kak w = 0

(A.8)

The reaction invariants x1 and x2 are found from the following mass balances: Ah

dx1 = fa Ca − (fa + fb + fc )x1 dt

(A.9)

Ah

dx 2 = fb C b − (fa + fb + fc )x 2 dt

(A.10)

where fc is the rate of water flows into the tank, A is the tank crosssectional area and h is the level of the tank. An overall mass balance on the tank yields Ah

dh = fa + fb + fc − fo dt

(A.11)

The exit flow rate fo is modeled as fo = Cv√h, where Cv is the constant valve coefficient. Finally, the pH is calculated using pH = −log[H+]

(A.12)

The simulated model comprises eqs A.8−A.12. Equations A.9 and A.10 are bilinear differential equations, while eqs A.8 and A.12 are highly nonlinear equations. Hence, this plant is a highly nonlinear process.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected], [email protected].

APPENDIX A The mathematical model of the pH process is derived using chemical reaction invariants, conservation law, and equilibrium relations according to Sedigh.30 In this system, acetic acid of concentration Ca flows in at a rate of fa and reacts with sodium hydroxide of concentration Cb which flows into the tank at a rate of f b. These chemical reactions can be expressed as

Notes

The authors declare no competing financial interest.

■

REFERENCES

(1) Camacho, E. F.; Bordons, C. Model predictive control, 2nd ed.; Springer-Verlag: London, 2004. (2) Allgower, F.; Zheng, A. Z. Nonlinear model predictive control: Assessment and future directions for research; Birkhäuser Basel: Germany, 2000. (3) Nelles, O. Nonlinear system identification: from classical approaches to neural networks and fuzzy models; Springer: Germany, 2008. (4) Wu, Q.; Zhang, Q.; Zong, Ch.; Cheng, G. Vibration control of block forming machine based on an artificial neural network. In Proceedings of the International Symposium on Neural Networks, China, 2007. (5) Barron, A. R. Universal approximation bounds for superpositions of a sigmoidal function. IEEE Trans. Inf. Theory 1993, AC-39, 930−945.

CH3COOH ↔ CH3COO− + H+ NaOH ↔ Na + + OH− (A.1)

The equilibrium constants of acid, base, and water are [CH3COO−][H+] CH3COOH

(A.4)

In addition to the chemical reaction for the acid and base, the solution must remain electrically neutral at all times, giving

■

Ka =

(A.3)

Two reaction invariants for these chemical reactions are the total ionic concentration of the acid and base:

VII. CONCLUSION In this paper, an algorithm is presented to overcome the difficulties in using the neural network model predictive controllers both in tuning problem and determining the neural network complexity. In the first step, a method based on Barron’s formula, bicoherence nonlinearity measure, and IQ frequency separation is provided to consider frequency characteristics of the process for calculating the number of neurons needed for plant identification. Then, a fuzzy rule based method is used to estimate the input weight parameter γ in the objective function of the NMPC controller. The efficiency and accuracy of the proposed method have been investigated using illustrative examples. The efficiency of the proposed algorithm is that having some simple prior information about the process behavior leads to a nonlinear predictive controller which can track a desired trajectory. It is only necessary to perform simple preliminary tests to obtain a few prior knowledge of the process. Designing training and test data, determining the number of hidden layer neurons, training the network, and tuning the weight parameters in the NMPC objective function will be automatically performed.

H 2O ↔ H+ + OH−

[OH−][Na +] NaOH

(A.2) 406



Article

(6) Suraweera, N. P.; Ranasinghe, D. N. Adaptive structural optimization of neural networks. ICTer 2008, 1, 33−41. (7) Son, J. S.; Lee, D. M.; Kim, I. S.; Choi, S. K. A study on genetic algorithm to select architecture of a optimal neural network in the hot rolling process. J. Mater. Process. Technol. 2004, 153−154, 643−648. (8) Shridhar, R.; Cooper, D. J. A tuning strategy for unconstrained siso model predictive control. Ind. Eng. Chem. Res. 1997, 36 (8), 729−746. (9) Neshasteriz, A. R.; Khaki Sedigh, A.; Sadjadian, H. Generalized predictive control and tuning of industrial processes with second order plus dead time models. J. Process Control 2010, 20, 63−72. (10) Iglesias, E. J.; Sanjuan, M. E.; Smith, C. A. Tuning equation for dynamic matrix control in siso loops. Ing. Desarrollo 2006, 19, 88−100. (11) Ali, E.; Zafiriou, E. Optimization-based tuning of nonlinear model predictive control with state estimation. J. Process Control 1993, 3 (2), 97−107. (12) Ali, E. Heuristic on-line tuning for nonlinear model predictive controllers using fuzzy logic. J. Process Control 2003, l3 (3), 383−396. (13) Ghazzawi, A. A.; Ali, E.; Nouh, A.; Zafiriou, E. On-line tuning strategy for model predictive controllers. J. Process Control 2001, 11 (4), 265−284. (14) Kawai, F.; Ito, H.; Nakazawa, Ch.; Matsui, D.; Fukuyama, Y. Automatic tuning for model predictive control: can particle swarm optimization find a better parameter? IEEE International Symposium on Intelligent Control, Singapore, 2007. (15) Isserlis, L. On a formula for the product-moment co-efficient of any order of a normal frequency distribution in any number of variables. Biometrika 1995, 82, 351−368. (16) Billings, S.; Voon, W. Correlation based model validity tests for nonlinear models. Int. J. Control 1986, 44, 235−244. (17) Nicholas, J.; Oslon, C.; Michalowicz, J.; Bucholtz, f. the bispectrum and bicoherence for quadratically nonlinear systems subject to non-gaussian inputs. IEEE Trans. Signal Process. 2009, 57, 3879− 3890. (18) Choudhury, M. A. A. S.; Shah, S. L.; Thornhill, N. F. Diagnosis of process nonlinearities and valve stiction: data driven approaches; Springer: Germany, 2008. (19) Nassirharand, A.; Karimi, H. Input/output characterization of highly nonlinear multivariable systems. Adv. Eng. Software 2002, 17 (10), 825−830. (20) Schwartz, C.; Gran, R. Describing function analysis using MATLAB and SIMULINK. IEEE Control Syst. Mag. 2001, 7 (6), 19−26. (21) Wechsler, H. Optimal sampling neural network performance and pattern recognition, Proceedings of the International Neural Network Society Annual Meeting, Washington, DC, 1995. (22) Zhu, Y. Multivariable system identification for process control; Elsevier Science Ltd: Oxford, 2001. (23) Astrom, K. J.; Hagglund, T. Revisiting the Ziegler-Nichols step response method for PID control. J. Process Control 2004, 14, 635−650. (24) Kocijan, J.; Banko, B.; Likar, B.; Girard, A.; Murray-Smith, R.; Rasmussen, C. E. A case based comparison of identification with neural network and Gaussian process models. IFAC International Conference on Intelligent Control Systems and Signal Processing, Portugal, 2003. (25) Yu, D. L.; Gomm, J. B.; Williams, D. Model predictive control of a chemical process using neural networks. UKACC International Conference on Control, 1998. (26) Wior, I.; Boonto, S.; Abbas, H. S.; Werner, H. Modeling and control of an experimental pH neutralization plant using neural networks based approximate predictive control; Aalborg University: Aalborg, Denmark, 2010. (27) Maghoul, P. Physical identification MIMO system & design & apply MIMO controller. M.Sc. thesis, K. N. Toosi University of Technology, 2005. (28) RT 532; Pressure control trainer datasheet, G.U.N.T. Gerätebau GmbH. (29) RT 512; Level control trainer datasheet, G.U.N.T. Gerätebau GmbH. (30) Sedigh, A. K. Universal control systems; Technical report; K. N. Toosi University of Technology: Tehran, Iran, 2007.

407


On the Structural Optimization of a Neural Network Model Predictive

Recommend Documents