Ind. Eng. Chem. Res. 1994,33, 2140-2150
2140
Nonlinear Rule Based Model Predictive Control of Chemical Processes Ching-Yu Peng and Shi-Shang Jang’ Department of Chemical Engineering, National Tsing-Hua University, Hsin-Chu, Taiwan, Republic of China
Several difficulties are still encountered in the direct use of a nonlinear model in the area of model based control. A time series rule based model is employed in this work to perform nonlinear control in cases where linear approaches have failed. T h e rule based model is basically comprised of a set of rules which are related to the time series of input and output data. The proposed control approach filtered out the high-frequency disturbances using possibility theory. An on-line identification phase is required if persistent changes of some parameters frequently occur. The identification algorithm maximizes the membership of the disturbance parameter in the immediate past. The control objective minimizes the square errors of the output and set point in a time horizon projected into the immediate future. Physical examples are simulated t o demonstrate the implementation of this approach. 1. Introduction
Linear control theory has found extensive application in the chemical process industry. Although chemical processes are generally nonlinear, a simple linear conventional controller can control the system adequately in most cases. Advanced control techniques are required in about 10% of all chemical processes, with its economic benefits being relatively profound. Model predictive control theory has been previously summarized by Garcia and Morari (1982). Most processeswhich require advanced control techniques are highly nonlinear processes. Approximate, or “nominal”, linearized models obtained via various empirical techniques (Richalet et al., 1978; Culter and Ramaker, 1979)have been conventionally employed; in addition, designed controllers must be detuned for uncertainty of such models. Alternatively, a trustworthy nonlinear model can be “globally linearized” (e.g., Hensen and Seborg, 1991), yet robustness consideration due to uncertainty in the nonlinear modelis stillnecessary. Direct on-line usage of a physical nonlinear model was proposed by Jang et al. (1987);however, this approach requires an on-line identification phase to compensate for model inaccuracies. Enormous computation power is required for performing both on-line identification and on-line optimization via a nonlinear model. A first principle model for a chemical plant is usually highly nonlinear, difficult to solve. Further, the model is always relatively inaccurate, containing several undetermined parameters. Such models are generally impractical for the purpose of on-line control. Several investigators have proposed the application of various simplified nonlinear modeling forms, e.g., neural network models (Hernandez and Arkun, 1992), or polynomial models (Patwarhan and Madhavan, 1993;Hernandez and Arkun, 1993) to perform predictive control. A set of rules which are related to the time series of input and output of the system is used in this work. The idea of rule based, or expert system, control has been originally proposed by Mamdani and co-workers (e.g. Procyk and Mamdani, 1979). Their algorithm is basically not related to process model, and rules are updated online by control performance rating. Their approach is adaptive, rather than “model predictive” in nature. Lu and Holt (1990) have suggested the employment of a time series rule based model, but they did not formulate a
“predictive control- methodology. In this work, a rule based model predictive control formalism complete with identification and filtering is derived. Simplified nonlinear models must be obtained either empirically or semiempirically in light of the fact that an accurate first principle model is currenty unavailable. The advantage of using a rule based model is that acquisition and validation of rule based knowledge has been a quite active area of research (e.g., Shieh and Joseph, 1992, Saraiva and Stephanopoulos, 1992). Moreover, uncertainty can be handled on the sound basis of fuzzy mathematics. A systematic strategy of determining the structure of an empirical nonlinear model is presented in this work using knowledge from the first principle model or extracting knowledge from process dynamics in a rule based form of current plant data. 2. Determination of Orders of a Time Series Rule
Based Model Given a process with the following lumped system dynamics:
wherex = ( x ~ , x ~ , . . . , x=~(y1,y2,...,y” ~,Y IT, u = (ul,uz,...,ul’)T, f = (fl,f2,...,f’)T, and g = (g1,g2,...g”’)T. If we assume that the controller is digital and the sampling time is constant and that
Equation 1 can then be converted into an equivalent discrete system, Le.,
where 4 = (@l,@2,...@9Tand y = (y1,y2,...,yi”IT. The above system dynamics is described in this work by an input/ output discrete system:
* Author to whom correspondence should be addressed. o s a s - ~ a a ~ ~ ~ ~ ~ ~ s ~ ~1994 - ~American ~ ~ o $ Chemical o ~ . ~ oSociety ~o
Ind. Eng. Chem. Res., Vol. 33, No. 9, 1994 2141
Definition: If the dynamics of a system can be described by eq 4, the integers n and m can then be defined as the output and input orders of the system, respectively. Analytical forms of f and g in (1) are known as first principle models and are typically rather difficult to obtain and solve. Likewise, analytical forms of 6 and 7 in eq 3 as well as F in eq 4 are also quite difficult to be defined explicitly. The F in eq 4 is represented here as a set of N time series rules in the following form:
F
0.015
yk 0.005
.."
* * q *, : , 11
13
15
17
19
k
Figure 1. Tracking performance of a full order rule based model, where the response with symbols of is the plant output and the response with symbols of + is the model predicted output.
The first order of business in constructing such a rule set clearly involves determining the input and output orders n and m. Selecting these orders subsequently determines how these rules should be written as well as the size of memory required to store them. Given some rough knowledgeof the first principle models (i.e., eq 3), the maximum input/output orders can be determined using the following fact: - - Definition: The characteristic matrices &, 1, and b of the discrete system (3) can be defined as:
following example. Example 1: The following discrete nonlinear system is considered. = 0.5~:
r,
i= [ajj] aij=
t
ny nonzero real number
if 4' is a function of xi if 4i is not a function of X' (5a)
B = [Bij]
t
B, =
ny nonzero real number
ift$ is a function of uj if $i is not a function of uj (5b)
Yk
The characteristic matrices for the above system are demonstrated as being:
&=[, O f d ] Yk+l
(5c)
-
b = [D,I
t
D, =
ny nonzero real number
y ( q ) = 1Fcq.I- i>-'B
+b]u(q)
a]
According to this above fact, the implicit model for this nonlinear system can be given as = f&ktYk-l,Uk)
(7)
A rule based model (Rl) with input/output orders of 2,l can predict the system behavior using a local interpolation algorithm given in the next section. Figure 1 shows the output trajectory and the model predicted trajectory subjected to an initial condition. The predicted trajectory sufficiently correlates with the real plant output. Example 2: A more complicated case is considered, = 0.5(~:)~ - 0.2x:
if y' is a function of uj if y' is not a function of uj
(5d) Fact: The input/output model (4) corresponding to (3) has its maximum input output orders as the following linear system, if 4 and y are invertiable to all Xk:
F=[O
= fyk + deyk-l + abuk
Yk+l
if y i is a function of X' if yiis not a function ofxi
E=[;]
The corresponding linear system therefore has an input order of 1 and an output order of 2:
F = [rij] rij= ny nonzero real number
2
= xk
+ o.luk
+
= 0 . 5 ~ : 0.4~:
In this case,
(6)
where q is the shift operator, such that qyk = Yk+l, 7 is an identity matrix. Proof: See Appendix. In the case that the system is with time delays, it is easy to perform the above analysis by adding some extended state variables (see, e.g. Astrom and Wittenmark, 1990). Application of the above fact is demonstrated in the
The corresponding linear system can then be expressed in the following form: Yk+l
= cf -I-c)Yk + d e h 1 + abuk - abfuk-1
This equation indicates that an input/output model (Rl) with orders of 2,2 can accurately predict the system
2142 Ind. Eng. Chem. Res., Vol. 33, No. 9, 1994
-0.0003 -0.OooL
-00006 -0.0007,
+ * + +
3
5
7
9
13
11
15
+ +
* +
17
19
k
Figure 2. Performances of output tracking using a full order ( R = 2, m = 2) rule baaed model and a rule baaed model with R = 2, m = 1, where the response with symbols of 0 is the plant output, the response with symbols of 0 is the full order model predicted output, and the response with symbols of + is the non-full order model predicted output.
behavior. This result is provided in Figure 2. However, if a model with input/output order of 1,2 is constructed, the model cannot predict the trajectory of output as shown in Figure 2. Given this basic structure of the rules, a flexible framework for representing system dynamics is made available. Only an approximate knowledge of the cause and effect relation between system changes and state variables has notably been used. Neither quantitative information concerning system dynamics nor the exact functional dependence of the cause and effect relation is necessary. In many cases, the effective input-output orders of a system are less than the maximum possible order. Actual plant data and statistical tools can be employed in further reducing the number of variables for the input/output model and fill in the details of the rule set. The burden of constructing an accurate first principle model can be substantially alleviated, despite the fact that a substantial amount of work is still required in terms of validating and simplifying this rule based knowledge. An adequate rule based model is assumed here to be able to be obtained using the above analysis and plant data. Model predictions can be performed via the method described in the next section. 3. Prediction of System Dynamics of a Nonlinear System by a Time-Series Rule Based Model and Possibility Theory
Given a set of time series rules of the form (Rl),and the current and immediate past trajectory, what is the model prediction of the output in the immediate future? Finding a rule that “precisely” matches the trajectory is usually not possible in light of the measuring error and the availability of only a finite set of discretized rules. It is next assumed here that the membership of each discretized state of either the observed or the manipulated variable Y (or U) is triangular as shown in Figure 3, i.e.,
if y < Y - E
if
y L Y + t
where E is an adequate tolerance depending on the system dynamics. Given noise-free data Yk-i, i = 1, ..., n, the
y
Y-E
Y
Y*E
Figure 3. Membership functions of a discretized state and a noisefree measurement.
Figure 4. Membership functions of a discretized state and a noisecontaminated measurement.
“possibility” of validity of the hypothesis of a given rule
Rj in rule set (Rl) is given by PRjyk+l(yk+l)
= min[Py,(Yk),***, Psr,+l(
Yk-n+l),Clu,( uk)”.Puh+l(Uk-m+l)
1 (9)
However, measurements would be noisy in most cases. According to many researchers (e.g., Civanlar and Trussell, 1986;Hohle, 19831,is noises are Gaussian, the uncertainties of a measurement g can be approximated by a near trapezoidal shape fuzzy membership (Figure 4), Le., PLYW=
f 0 y -9
if y < g - b x a
+ ( b - c)a
( b - c)a
1: 9
+
if j ~ - b X a l y < y - c X a if 9 - c X a l y < g + c X a
(b-c)a-y
( b - c)a
if 9 - c x a 5 y 5 9 + b x a if y L g + b X a
(10) with a being the average noise standard deviation of measurement error, and b is an appropriate constant so that the confidence level >0.999 and c, confidence level < 1/21/2 (Civanlar and Trussell, 1986) ( b = 2.5, c = 0.4 here). It should be noted that eq 10 is not unique for Gaussian density function. This is a very active area for fuzzy mathematics researchers. Non-Gaussian noise can be treated in a similar way. Therefore, the “possibility” of the measured variable 9 in state Y is obtained using the extension principle in fuzzy mathematics (Zimmerman, 1991) by (as shown in Figure 4)
PJ,O = ~ a x ( m i n [ ~ & ) , ~ ~ W l l Y
(11)
Consequently, the validity of the hypothesis of a given
Ind. Eng. Chem. Res., Vol. 33, No. 9, 1994 2143 o observations x forcasted behavior
I
With the possibility membership of the future state calculated in eq 12, a defuzzified prediction OfQk+lis given by
' tk-p
-Hl
1
immediate post
Equation 13 is an optimal estimator in the sense that the memberships of the noisy past n + 1 observations are maximized for all of the rules (see Figure 4). 4. On-Line Identification If any slow varying disturbance inputs occur that can be detected off-line or measurable by slow on-line sensors, e.g., a gas chromatograph, so that the data cannot be used by model predictive controllers, the system equation could be rewritten as
i
;
tk-n.,tk-n*2
~
i
tk.3tk-2tk-l
l
F'
-PT
i
l
i
tk*l tk+2 tk*3
present time
I
tk+Q
H,=QTl immediate future
Figure 5. Scenario of the horizon approach.
at i step before is denoted as di. Using the prediction procedure described in the previous section, it is possible to obtain a prediction of the current state 9; denoted by
96(d') = F(Stk-l,".,9k-n,Uk-l,...,uk-mtdi)
(17)
An optimal estimator of the current state, as well as the disturbance vector, maximizes the membership of each estimate state 9; by the following extension principle:
~ ~ ~ (=9max(min[~~ki(u),~yk(u)l~ ;) (18) Y
A filtered value of the current state should be given by P
The rule based model, in this case, should include this low frequency disturbance input using the following form: If Y k - n is in state Yk-,,, and Yk-,,+l is in state Yk-,,+l, ..., and Y k - 1 is in state Yk-1, and if Uk-,, is in state Uk-m, and Uk-,,+l is in state Uk-m+l, ...,and uk-1 is in state Uk-1, and if d is in state D,then yk is in state Yk. (R2) Given this rule based model, the disturbance can be identified using on-line data. The d is assumed here to be discretized into R states, for a specific discretized state of disturbance Dr, and a given set of measurement 3k-1, 3k-2, ...,Qk-n, U k - 1 , Uk-2, ...,Uk-m, the predicted current state 9; corresponding to Dr using the prediction procedure described in the revious section. A defuzzified value of the disturbance k is calculated as
:
R
However, d and 9 k cannot be determined by a single trajectory since the measurements are noisy. A possible method of solving this phenomenon involves employing a horizon approach, as shown in Figure 5. Consider a sufficient long identification horizon (i = 1, ...,P ) in the immediate past. The estimation of the disturbance vector
Given f k , one can use eq 16 to find a filtered
ak
5. Optimal Control Figure 5 considers a Q-step time horizon in the immediate future which is termed as the optimization horizon. The most general model based control strategy can be formulated as an optimization problem:
2144 Ind. Eng. Chem. Res., Vol. 33, No. 9, 1994
The above optimization problem can be solved by the followingprocedure. Consider a filtered trajectory $'k-n+l, $'k-n+2, ..., $'k in the immediate past together with the estimated disturbance values 2 and forecasted immediate future system behavior Yk+l,Y k + 2 , ...,Yk+Qgiven uk,uk+l, ..., uk+Q-l. The exhaustive search of all feasible control policy actions can be performed by forward branching search dynamic programming. This forecasted behavior can notably be obtained by the methodology described in section 2. However, fine-tuning control actions is still necessary since uk, uk+l,..., uk+Q-1 are discretized. The following is proposed here: let pi be the membership of a rule with an objective function of *i 1-
*i( yk+l, yk+2,**',
yk+Q, uk,uk+l,**',
+
I
I
I .1
uk+Q-l)
*max
if
h(Yk+l,Yk+2,"',Yk+Q,Uk,Uk+l""tUk+gl) 2 0
0 if h(Yk+l,Yk+2,"',Yk+Q,Uk,Uk+lt"',Uk+gl) 0.999 c = standard deviation of measurement error;confidence level
< 1/21/2
d = disturbances
D = characteristic matrix 7 = identity matrix
x:+1
= $i(Zk,Xf)
Xj+l
= +j(Zk,Xf)
(AI)
then there exists a function 4 such that 4+2
provided
$i,
$j
= #(zk+l,zk)
(A21
are invertible to xf, where Zk = [ x i , -
4,...,xpl,xp!...,Xk"1T. xi
hence zjk,+l = $j(zk,$;l(zk,x6+1)), = (P(zk,xi+I) or xk+l = v-'(zk,x:+1). $j(zk+l~f+1) = $ j ( z k + l t d ' ( z k , ~ ~ + l )= ) 4But, ~:,+2 (zk+l,zk,d+l)= '?(zk+l,zk). 2. Note: If = $i(~k,xf), but x $ + ~= $ j ( Z k ) , then there is no need to perform the above elimination; one has xjk+2 = $j(zk+l)* 3. Consider the following system,
Proof:
=
I);~(Z~,Z&~),
Xk+l
= $(Xk)
(A31
+
where = [$1,$2, ...,$ J T , a variable xf can be eliminated from this nonlinear state space equation if $i,$j are invertible for xf v j such that $j = $.(xk,zf), but the order of the system is increased by 1. +here exists function vector 4 such that state space system (2) can be rewritten as zk+2
= 4(zk+I,zk)
(A4)
provided that a t least one of the elements of $ contains x;, where 1 2
zk = [xk,xk,...,
,...,xkln
zi-l xi+l
k
3
k
T
T
4 = [41,42,"',4n-11
4. The above elimination can be repeated until only one variable xi remains; it yields
2150 Ind. Eng. Chem. Res., Vol. 33, No. 9, 1994 i
= f(xk+n-1,x;+,-2,...,Xk)
i
(A51
Note that some terms may be missing in light of the fact that some variables that to be eliminated from the variable set are not of functions of all other variables. In that case, the order of (A51 is reduced as shown in (2). 5. Property 2. Consider a corresponding linear system of (A3) xk+l
=
(A61
then there exists a linear system equivalent to (A6) zk+2 = 4'2,
= Alzk+l -k
A2Zk
= (qA1 + AZ)zk (A71
or (q21- QA1- A2)zk = 0
(-48)
where A1 = [A:]
Then
A!.'I =
e e
A; =
A, = [A;]
, if di is not a function of 2;1k+1 any real number, if di is a function of Z J ~ + ~ , if di is not a function of zi
any real number, if r # ~ is ~ a function of zi
i = 1, ...,n - 1
ProoE The same procedure as with the proof of property 1 is used. 6. The above elimination can hence be repeated until only one variable x i remains. This yields (qn - q"-lA1 - qn-zx,
- ...-;.,,A
=0
(A91
Note that the above procedure is actually the same as Gauss-Jordan elimination of a variable from a set of linear system. The following is a well-known fact in linear system theory. 7. Lemma 1. Consider a linear system (A6) corresponding to (A3); if ( q l - A) is invertible, then C ( q l A)-1 xf = 0 is equivalent to (A9), where C = [O,...,u,...,01, a is any real number and is located at the ith position of the vector. All real numbers are assigned in a way that cannot be eliminated by subtraction and addition. 8. Lemma 2. Consider the following nonlinear system: xk+l
(A101
= $(xk&k)
yk = Y(xL,uk)
where y and u are scalars, and corresponding linear system can be expressed as
+ BUk yk = cxk + D U k
xk+l
(All)
= axk
Then there exists an equivalent input/output system with maximum orders equivalent to the followinglinear system: Y ( Q ) / u (= ~ )C(q1- a)-' B
provided that (q1- a) is invertible.
+D
(A121
Proof: (A12) can be easily obtained using lemma 1 and the previous procedures of Gauss-Jordan elimination and nonlinear variable elimination. 9. Lemma 2 can be easily extended to the more general case in the Fact.
Literature Cited AstrBm; Wittenmark. Computer Controlled Systems, Theory and Design, 2nd ed.; Prentice-Hak Englewood Cliffs, NJ, 1990. Civanlar, M. R.; Trussell, H. J. Constructing Membership Functions Using Statistical Data, Fuzzy Seta and Systems. 1986,18,1-13. Culter, C. R.; Ramaker, B. L. Dynamic Matrix Control: A Computer Control Algorithm. Presented at the Eighty-Sixth National Meeting of the AIChE, Houston, April 1979. Economou, C. G.; Morari, M.; Palsson, B. 0.Internal Model Control. 5. Extension to Nonlinear Systems. Znd. Eng. Chem. Process Des. Dev. 1986,25,403-411. Garcia, C. E.; Morari, M. Internal Model Control 1. A Unifying Review and Some New Resulta. Znd. Eng. Chem. Process Des. Dev. 1982,21 (2),308-323. Henson, M. A.; Seborg, D. E. An Internal Model Control Strategy for Nonlinear Systems. AZChE J. 1991,37 (7),1065-1081. Hernandez, E.; Arkun, Y. Study of the Control-Relevant Properties of Back Propagation Neural Network Modela of Nonlinear Dynamical Systems. Comput. Chem. Eng. 1992,16(4),227-240. Hernandez, E.; Arkun, Y. Control of Nonlinear Systems Using Polynomial ARMA Models. AZChE J. 1993,39 (3),446-460. Hohle, U. Fuzzy Measures as Extensions of Stochastic Measures. J. Math. Anal. Appl. 1983,92,372-380. Jang, S. S.; Chen, L. Q. Nonlinear Model Predictive Control Based on an Impliciit Model. Chin. Znst. Chem. Eng. 1992,23(4),241251. Jang, S. S.; Joseph, B.; Mukai, H. On-Line Optimization of Constrained Multivariable Chemical Processes. AZChE J. 19878, 33 (l),26-42. Jang, S. S.; Joseph, B.; Mukai, H. Control of Constrained Multivariable Nonlinear Process Using a Two Phase Approach. Znd. Eng. Chem. Process Des. Dev. 1987b,26,2106-2114. Jang, S. S.; Wong, D. S. S.; Wong, S. J. Optimal Robust Linear Controller Design for Chemical Processes Using An Extended Regional Mapping Approach. Chem.Eng. Sci. 1991,47(a),20572068. Joseph B.; Wang, F. H.; Shieh, D. S. S. Exploratory Data Analysis: A Comparison of Statistical Methods with Artificial Neural Networks. Comput. Chem. Eng. 1992,16 (4),413-423. Lu, Z.H.; Holt, B. R. Nonlinear Robust Control: Table Look-Up Controller Design. Proc. Am. Control Conf. 1990,2758-2763. Patwarhan, S. C.; Madhavan, K. P. Nonlinear Model Predictive Control Using Second-Order Model Approximation. Znd. Eng. Chem. Res. 1993,32,334-344. Procyk, T.J.; Mamdani, E. H. A Linguistic Self-organizing Process Controller. Automatica 1979,15,15. Richalet, J.; Rault, A.; Testud, L.; Papon, J. Model Predictive Heuristic Control: Application to Industrial Process. Automatica 1978,14,413. Saraiva, P. M.; Stephanopoulos, G. Continuous Process Improvement Through Inductive and Analogical Learning. AZChE J. 1992,B (2),161-183. Shieh, D. S. S. Automated Knowledge Acquisition Form Routine Data For Process and Quality Control. Ph.D. Thesis, Washington University, St. Louis, 1992. Zimmermann, H. J. Fuzzy Set Theory and Its Application; Kluwer Academic Publishers: London, 1991.
Received for review January 4, 1994 Revised manuscript received May 19, 1994 Accepted June 27, 1994. 0 Abstract published in Advance ACS Abstracts, August 15, 1994.
