BACKGROUND ART There is a number of known non-linear control methods that were designed with the aim of achieving more efficient control of non-linear processes than the conventional linear methods such as the well-known and industrially proven PID controller (Bequette 1997, Henson & Seborg 1997, Qin & Badgwell 1998). All of those methods are based on some form of modelling of the process non-linearity and using this knowledge for improving the control strategy by applying non-linear controller action. Some of them are based on detailed theoretical (first-principles) modelling, while other use experimental modelling based on general non-linear models such as artificial neural networks, fuzzy models, wavelets, multiple models etc. (Sjoberg & al. 1995, Murray-Smith & Johansen 1997) where model parameters are extracted from measurement data using a form of model identification (Ljung 1987), or controller parameters are calculated from measurements directly.

Self-tuning or adaptive control methods (Astrom & Wittenmark 1995) assume time-varying process dynamics and therefore calculate model parameters, or update the controller directly, during controller operation. With adaptive controllers this usually occurs continuously while with self-tuning controllers it is triggered on demand. In both cases, injection of perturbation signals may be used. While most adaptive methods assume that the process is linear, this concept is also applicable to non-linear processes.

An overview of existing patents in related fields includes: WO 01/79945 to Hess et al. describes a non-linear controller-based on an adaptive-linear model-- predictive controller for which the model is obtained from a first-principles non-linear : process model.

US5335164 to Gough, Jr. et al. describes a method and apparatus for non-linear adaptive control based on Laguerre series network model and predictive control.

US6055524 to Chen describes a direct adaptive non-linear controller based on an artificial neural network.

US5740033 to McCroskey et al. describes a control system based on a non-linear model predictive controller and an interactive modeller.

US5682309 to Fontaine et al. describes a non-linear controller based on a non-linear state-space model, state estimation and on-line optimisation using sequential quadratic programming.

WO 01/92974 to Calise et al. describes an adaptive non-linear controller based on a linear controller and model inversion of an adaptive artificial neural network.

US5394322 to Hansen describes a model-based method for self-tuning of a linear PID controller.

EP0571080 to Molnar describes a method for self-tuning of a linear P1D controller.

EP0360206 to Yokokawa et al. describes a method of linear PID controller self-tuning based on three performance indexes.

US4407013 to Arcara et al. describes a method of linear P1D controller self-tuning from a discrete-time model.

While advanced in theory, there are several drawbacks to known methods and controllers in practice. Due to the complexity of implementation and computational demand, the hardware requirements are high. The complexity of tuning makes them unattractive to non-specialised engineers. This is particularly the case with methods based on detailed theoretical modelling, while reliability of experimental modelling is often questioned. In any case, control is subject to imperfect modelling.

The presented controller attempts to overcome these problems and provide an efficient and user- friendly tool for control. of a certain practically very important class of non-linear processes that may be represented by a fuzzy set of low-order local linear models that are selected using a scheduling variable. The model is obtained by means of experimental modelling using a special on-line learning procedure combining model identification with pre-and post-identification steps that provide reliable operation. There is a choice of several control algorithms suitable for different processes that may be used tor control and whose parameters are automatically tuned from the model. The controller monitors the resulting control performance and may react to detected irregularities. All the algorithms are adapted for implementation on standard low-cost industrial hardware platforms such as programmable logic or open controllers, accompanied by software that simplifies the initial configuration from a personal computer.

SUMMARY OF THE INVENTION The self-tuning non-linear controller according to the present invention is intended for control of a class of non-linear processes that may be represented by a fuzzy set of low-order local linear models that are selected using a scheduling variable. The model is obtained by means of experimental modelling using a special on-line learning procedure combining model identification with pre-and post-identification steps that provide reliable operation. There is a choice of several control algorithms suitable for different processes that may be used for control and whose parameters are automatically tuned from the model. The controller monitors the resulting control performance and may react to detected irregularities. The controller is suitable for implementation on hardware platforms such as programmable logic or open controllers.

The self-tuning non-linear controller is based on a multi-faceted model comprising a set of first- order discrete-time linear local models with time delay and offset and second-order discrete-time linear local models with time delay and offset which are placed at different positions of a scheduling variable s, that can be as a whole interpreted as the said Takagi-Sugeno fuzzy model.

The self-tuning non-linear controller consists of a Signal Preprocessing Agent, an Online Learning Agent, a Model Information Agent, a Control Algorithm Agent, a Control Performance Monitor, and an Operation Supervisor, where the Signal Preprocessing Agent prepares the process signals for real-time control and stores them to a buffer for, batch-wise post-processing, 'the Control Algorithm Agent performs non-linear control based on the multi-faceted model in real time and automatically retunes the controller from the multi-faceted model, * the Online Learning Agent periodically inspects the buffer of recent real-time signals and extracts new model information when appropriate, the Model Information Agent maintains the active multi-faceted model, the Control Performance Monitor periodically inspects the buffer of recent-real-time-signals- and extracts features of closed-loop performance of detected events, that may be used to trigger modification of Operation Supervisor configuration, 'the Operation Supervisor coordinates controller operation.

DISCLOSURE OF THE INVENTION The invention will be described with reference to the accompanying drawings, illustrating: Fig. 1: Self-tuning non-linear controller overview Fig. 2: Fuzzy membership functions of local models in the Multi-Faceted Model Fig. 3: Online learning procedure Fig. 4: Fuzzy Gain-Scheduling Controller Control Algorithm Agent overview Fig. 5: Dead-Time Compensation Controller Control Algorithm Agent overview Fig. 6: Rule-Based Neural Controller Control Algorithm Agent overview Fig. 7: Control Performance Monitor Fig. 8: Fuzzy logic membership functions for Performance Index calculation in CPM Fig. 9: Fuzzy logic rule table for Performance Index calculation in CPM Fig. 10: The sample neutralisation process Fig. 11: pH control performance comparison between FGSC and linear PI controller 1. Overview The proposed controller is built in the form of a multi-agent system where several independent agents (modules) interact with each other. The system includes the following agents, as shown in Fig. 1: Signal Preprocessing Agent (SPA) (1), Online Learning Agent (OLA) (2), Model Information Agent (MIA) (3), 'Control Algorithm Agent (CAA) (4), 'Control Performance Monitor (CPM) (5), Operation Supervisor (OS) (6).

2. Multi-faceted Model The model-based controller is founded on a multi-faceted model (MFM) (7) that includes several-model forms required by the online-learning mechanism and control algorithm agents.

Specifically, it includes, a set of first-and second-order discrete-time linear local models with time delay and offset, which are placed at different positions of a scheduling variable s. The model equation for second order models is y (k) =-a1, jy (k-1)-a2 jy (k-2) + bl ju (k-1-du y) + +b2,ju(k-2-duj)+c1,jv(k-1-dvj)+c2,jv(k-2-dvj)+rj and for first order models simplified to y (k) =-a1, jy (k-1) +b1,ju(k-1-duj)+c1,jv(k-1-dvj)+rj where k is the discrete time index, - j is the number of the local model, y (k) is the process output signal (controlled variable, CV), u (k) is the process input signal (manipulated variable, MV), v (k) is the (optional) measured disturbance signal (MD), du is the delay in the MV-CV path, dv is the delay in the MD-CV path, ai, bi, ci and r are parameters of local models.

The MFM can be interpreted as a Takagi-Sugeno fuzzy model, for the second order written as where =/) is the degree of fulfilment of the j-th membership function (it is a function of the scheduling variable s). Normalised triangular membership functions are used as illustrated in Fig. 2, so that at any value of s.

The scheduling variable is calculated in each time instant as follows s (k) = krr(k) + kyy (k) + kuu (k) + kvv (k) where the coefficients kr, ky, ku, and kv are predefined constants.

3. Online Learning Agent OLA (2) analyses the buffer of recent real-time signals (8), prepared by SPA, and generates new local models (9) when conditions are appropriate. It is invoked periodically or upon demand by OS. Since it is a computationally intensive task, it-runs as a low-priority task in a multitasking- operating system.

The following procedure is carried out whenever OLA is invoked, as illustrated in Fig. 3: Excitation check (10), Copying of the current MFM from MIA (11), Selection of local models (12), Identification of selected local models (13), Verification/validation of selected local models (14).

Submission of modified local models to MIA (15).

Excitation check: if the variance of the signals r (k), y (k), u (k), and v (k) in the active buffer is below their specified thresholds, the execution is cancelled (16).

Selection of local models : the local models for which the sum of their corresponding membership functions ß, {s) over the active buffer normalised by the active buffer length exceeds 0.2 are selected. Further processing does not include other local models.

Identification of local models: the parameters of the selected local models are identified using the Fuzzy Instrumental Variables (FIV) identification method, an extension of linear instrumental variables identification procedure (Ljung 1987) for the specified MFM.

The identification is performed for each local model (and its corresponding fuzzy set, denoted by the indexj) separately (17). The most general form of the procedure where all parameters of a second-order model are estimated is described below.

Firstly, the vector of the estimates of the plant parameters 6j= [a1,j,a2,j,b1,j,b2,j,c1,j,c2,j,rj] is initialised (18) from the active MFM in MIA, and the covariance matrix Pj is initialised to 105 I (identity matrix).

Then, an approximation of Pj and #j is obtained using weighted least-squares (LS) identification (19) in the recursive form (to avoid matrix inversion) with a dead-zone preventing result degradation by noise : for k = (buffier_start + max (duj, dvj) +2) to buffeo eizd <BR> <BR> <BR> #j(k+1)=[-ßjy(k),-ßjy(k-1),ßju(k-duj),ßju(k-1-duj),ßjv( k-dvj),ßjv(k-1-dvj),ßj]T<BR> <BR> <BR> <BR> <BR> <BR> <BR> <BR> <BR> if abs(ßjy(k+1)-#jT(k+1)#j(k)) > noise_threshold<BR> <BR> <BR> <BR> <BR> <BR> <BR> <BR> <BR> <BR> #j(k+1)=#j(k)+Pj(k)+Pj(k+1)#j(k+1)(ßjy(k+1)-#jT(k+1)#j(k))& lt;BR> <BR> <BR> <BR> <BR> <BR> Pj (k)#j(k+1)#jT(k+1)Pj(k)<BR> <BR> <BR> <BR> Pj(k+1)=Pj(k)-1+#jT(k+1)Pj(k)#j(k+1) endif endfor where"jT iS the vector of measurements.

Finally, improved estimates of Pj and Q are calculated using weighted instrumental variables (IV) identification method (20) for k = (buffer start + max (duj, dvj) +2) to buffer end vj (k+1) [-pjy (k),-pjy (k-1), Pju (k-duj), Pju (k-1-duj), Pjv (k-dvj), Pjv (k-1-dvj), PjIT m y (k+i) _ ; T (k+1) e ; (k) j=l Xj (k+1) = [-ßjy (k),-ßjy (k-1), ßju (k-duj), ßju (k-1-duj), ßjv (k-dvj), ßjv (k-1-dvj), ßJ] T if abs (jy (k + 1) _ v jT (k +1) 0./(k)) > noise threshold ûj (k + 1) = Oj (k) F p, + i) x, +1) (/ ?, 3 +1)-/ ( + i) e, () P (k (k) x (k+1) yrT (k+1) P (k) 1+WjT (k+I) Pj (k) Zj (k+l) endif endfor where X ; (k+1)T is the vector of instrumental variables and y (k) is the simulated output.

For identification of first order models a very similar procedure is used where the following vectors are modified: #j = [a1,j,b1,j,c1,j,rj] #j(k+1)=[-ßjy (k), ßju(k-duj), ßjv(k-dvj), ßj]T xj(k+1) = ßj# (k), ßju(k-duj),ßjv(k-dvj), ßj]T - In case of lack of excitation in the model branch from v (MD) to y (CV), or when MD is not present at all, variants of the method with the following vectors are used for second order models <BR> <BR> <BR> #j=[a1,j,a2,j,b1,j,b2,j,rj]<BR> <BR> <BR> <BR> <BR> <BR> #j(k+1)=[-ßjy(k),-ßjy(k-1),ßju(k-duj),ßju(k-1-duj),ßj]T <BR> <BR> <BR> <BR> <BR> xj(k+1)=[-ßj#(k),-ßj#(k-1),ßju(k-duj),ßju(k-1-duj),ßj]T and for first order models, respectively, #j=[a1,j,b1,j,rj] #j(k+1)=[-ßjy (k), ßju(k-duj), ßj]T xj(k+1)=[-ßj#(k), ßju(k-duj),ßj]T In case of lack of excitation in the model branch from u (MV) to y (CV), variants of the method with the following vectors are used for second order models <BR> <BR> <BR> #j=[c1,j,c2,j,rj]<BR> <BR> <BR> <BR> <BR> #j(k+1)=[ßjv(k-dvj),ßjv(k-1-dvj),ßj]T xj(k+1)=[ßjv(k-dvj),ßjv(k-1-dvj),ßj]T and for first order models, respectively, #j = [c1,j,rj] #j(k+1) = [ßjv(k-dvj), ßj T XX (k+1) = [ßjv(k-dvj), ßj]T Verificationlvalidation of local models : This step is performed by comparing the output of a selected local model, denoted by the index j, with actual process output in the proximity of the local model position, specified by the corresponding membership functions Pj, and calculating the normalised sum of mean square errors (MSE) vaf. When performed on the same data that was used for parameter identification it is called verification, while with model parameters obtained on another data set it is referred to as validation.

The local model is described as #j = [a1,j,a2,j,b1,j,b2,j,c1,j,c2,j,rj]; in case of a first- order model or when measured disturbance v is not present, the coefficients a2,j,b2,j,c2,j are set to 0. The following procedure is used: six stars = max (duj + dvj ? + model_order ; ym_dd=y(sim_start-2) ; ym_d=y(sim_start-1) ; Safe = 0 ; sumbetaj = 0 ; for k buffer_start + sim_start to buffer end calculate membership function<BR> p function ßj = ßj(s(k)) sumbetaj = sumbetaj + ßy ; ym = [-ym_d-ym_dd u(k-duj-1) u(k-duj-2) z(kdvj-1) z(k-dvj-2) 1]T##j; ym_dd = ym_d, ym_d = ym ; vaf = vaf + (y (k)-ym) z * (3 ; end vaf = vaf I (sumbetaj + 10-lao) ; For each of the selected local models, this step is carried out with three sets of model parameters O : the initial model from MIA (21), the model estimate of the LS method (22) and the model estimate of the IV method (23). The model with the lowest vafj of the three is sent to MIA (15) as the result of online learning, along with its vaf that serves as a confidence index (cij) and a flag chgj indicating whether the model is new or not.

4. Model Information Agent (MIA) MIA maintains the active MFM (model-in-use) and its status information.

MIA accepts online learning results from OLA. If the model confidence index is below the specified threshold, the automatic mode may be disabled. When a new local model is received, it is accepted if it passes the stability test (and possibly other tests) and its confidence index is sufficient. If the local model is accepted, its agj flag is set, which means that it is ready for controller tuning by the CAA's tuning layer and not used yet. Each local model also has a scj flag indicating whethei the local model had been tuned since start-up or not.

MIA contains a mechanism for inserting additional local models.

MIA may also store the model-in-use to a local database (24) or recall a previously stored one from there.

At initial configuration, MIA is filled with local model that are not accurate but provide a reliable (although sluggish) performance of the controller, similar to the Safe mode. With online learning through experiments or even normal operation (when the conditions are appropriate), an accurate model of the plant is built by receiving identified local models from OLA.

5. Control Algorithm Agent The CAA (4) includes functionality of an advanced non-linear control algorithm and automatic tuning of its parameters from the MFM. Several CAAs with different operating principle and properties may be used in the controller and may be interchanged in the initial configuration phase.

The CAA supports several modes of operation : Manual 7node : open-loop operation Safe mode : a fixed linear controller, tuned conservatively Auto mode (or several auto modes with different tuning parameters): non-linear controller The CAAs share a common interface of interaction with other agents of the controller and a common modular internal structure, consisting of three layers: The con7-7 ol layer (25) includes the functionality of a local linear controller (or several local linear controllers at once), including all functionality required for industrial controllers, such as anti-windup protection, reasonable switching between open-loop and closed-loop operation, etc.

The scheduling layer (26) performs switching or scheduling (blending) of tuned local linear controllers, so that in conjunction with the control layer a fixed-parameter non-linear controller is realised.

The tuning layer (27) performs automatic tuning of controller parameters from the MFM whenever new models are generated (indicated by the ag flags) and auto-tuning is enabled.

When tuning is completed, the ag flags are cleared and the tuned flags are set. The tuning layer maintains the controller parameters, that are derived from the MFM in general but CAA-specific, and the controller parameters also include a copy of the MFM from which they were tuned. It replaces the parameters of the control and scheduling layer in such manner that their real-time operation is not disturbed.

Fuzzy Gain-scheduling Controller (FGSC) C, 4A An overview of FGSC is shown in Fig. 4.

The control layer of FGSC (28) consists of a single PID controller in the velocity form (in order to use velocity-based linearisation when blending local controllers), equipped with anti-windup protection, and is defined by the following procedure carried. out in each discrete time step k : % INPUTS : % Tsámp-sampling time of the PI control algorithm % r (k)-reference signal % y (k)-process output signal %limit_umax,_umin,_duup,_dudown-actuator limitation @ % controller controller mode- % u_manual-controller output for the manual mode % Kp, Ti, Td-Current PID controller parameters % beta-0 : disturbance rejection, 1 : reference following) % KTd-The filter time constant ratio, KTd = Td I Tf) % OUTPUTS: % u (k)-controller output % INTERNAL SIGNALS (previous values stored for the next sample): % uP_1, uD_2, uD-1, ucu_1, r_1, uu_1, u_1 % Re-calculation of all filter and gain parameters if Ti #= 0 Kpi = Tsamp / Ti ; else Kpi=0 ; endif Kpd = Td / Tsamp ; Tf = Td / KTd; Kf=1 - exp(-Tsamp / Tf) ; % Calculation of anti-windup parameters if (beta=1) Tmp = 1 + Kpi + Kpd ; else Tmp = Kpi; endif Tmpl = 1 / (Tmp + 1e-10); Daw = Tmpl / (Kp * Kf + 1e-10) ; if (beta == 1) Kawl=Kpi * Tmpl,'Kf Kaw2 = 1Kf-Kpd * Tmp1 ; else Kaw1 = 1 / Kf; K"2 = 1-Kf ; endif % Calculation of PID terms e (k) r (k)-y (k) ; uP = Kp * e (k) ; uI = Kpi * Kp * e (k) ; uD = Kp * Kpd * (e (k) + (beta-1) * r (k)) ; % Calculate unlimited unfiltered controller output uu (k) ucu 1 + uI + uD - 2*uD_1 + uD_2 + uP - uP_1 + Kaw1 * (u_1 - uu_1) - Kp * (1-beta) * (r (k)-r 1) ; uu (k) = Kf* ucu + Kaw2 * (u_1 - uu_1) + (1-Kf) * 1 ; % Calculate the new value of u (k) (limited controller output) if controller-mode Manual u (k) = u_manual ; else % Automatic control u (k) = uu (k) ; endif % Controller output limitation u (k) = max (min (u (k), limit_umax), limit_umin) ; duout= u (k)-u_1 ; if duout > Tsamp * limit_duup u (k) = u 1 + Tsanzp * limit_duup ; elseif duout <-Tsamp * limit dudown u (k) = u_1 - Tsamp*limit_dudown ; endif % Update internal signals for next sample uP_l=uP ; uD 2=uD_l ; uD_1 = uD; ucu_1 = ucu; r_1 = r(k); u_1 = u(k); uu_1 = uu(k) ; The scheduling layer of FGSC (29) is implemented as fuzzy blending of the controller parameters (in case of Ti, its inverse value) according to the scheduling variable s and the triangular membership functions ßj (30) of the local models in the MFM. The following procedure is used: % INPUTS : % CPM-controller parameters matrix CPM = [CPI CP2... CP",] % where m is the number of local controllers, sorted by scheduling variable % CP ; are local PID controller parameters, CPj = [Kpj Tij Tdj sj]T % s-the current scheduling variable % OUTPUTS: % CP = [Kp Ti Td s]T - scheduled controller parameters at the current scheduling variable s % Find two neighbouring controllers j = 1; s_higher = s1 ; while (s > s_higher) & (j <= m) j = j + 1 ; ifj <= m s higher = sj ; endif endwhile % Determine scheduled parameters if j = 1 % s is lower than s1 CP = CPI ; elseifj = m + 1 % s is higher than Sm CP = CPm ; else % s is between si and s,, actual fuzzy scheduling tmp = (s-sj-1)/(sj-sj-1) ; Kp = Khi-1 + tmp * (Kpj - Kpj-1); Td = Tdf l + tn7p * (Tdj - Tdj-1); if (Tij == 0) # (Tij-1 == 0) Ti = 0; % PI controller else Ti = 1/( 1/Tij-1 + tmp* (1/Tij - 1/Tij-1)) ; endif endif The tuning layer of FGSC (31) is based on the magnitude optimum (MO) criterion implemented using the multiple integration (MI) method (Vrancic & al. 2001). Discrete-time process parameters that are extracted from the MFM are converted it into continuous-time parameters. Then the so-called areas are calculated which are used for the calculation of the controller parameters, first the PI and then PID controller parameters. The following procedure is used for tuning of a single local PI controller from its corresponding local model (32) : % INPUT PARAMETERS : % Tsamp-sampling time % beta-1 : reference following, 0: disturbance rejection (boolean) % KTd-The ratio between Td and Tf % Kid-The ratio between Ti and Td in the serial PID structure % Kmax-maximum allowed open-loop gain (controller gain * steady-state process gain) % LMj = [au, a2j,b1j,b2j,duj,sj]- local model parameters, delay and scheduling variable % position, extracted from the MFM (second order models: a2j = b2j = 0) % OUTPUT PARAMETERS: % CPj-controller parameters, CPj = [Kpj Tij Tdj sj] % Flag_PI gain-The PI gain was modified from Kmax ; boolean % Transformation from discrete to continuous model by using bilinear transformation Kpr = (b1j + b2j) / (1 + a1j + a2j); b1c = - b2j * Tsamp / ( b1j + b2j) ; b2C = 0.25 * (b2j - b1j) * Tsamp2 / (b1j + b2j) ; aie = (1-a2y) * Tsamp/ (1 + au + azid ; a2, = 0.25 * (1 + a2j -a1j) * Tramp2/ (1 + au + a2y) ; Td = duj * Tsamp ; % Calculation of areas from continuous model for magnitude optimum calculation AO = Kpr ; A1 = Kpr *(a1c - b1c + Td) ; A2 = Kpr * (b2c - a2c - Td* b1c + 0. 5*Td2) + A1* a1c ; A3 = Kpr*(Td* b2c - 0.5*Td2* b1c + 0.166666*Td3) + A2* a1c - A1* a2c ; A4 = Kpr*(Td2 b2c *0. 5-0. 166666*Td3* b1c + 0. 04166666*Td) + A3* a1c - A2* a2c ; A5 = Kpr*(Td3* b2c *0.166666-0. 04166666*Td4* bic + 0. 008333333*Td5) + A4* a1c - A3* a2c ; % Calculation of PI parameters for magnitude optimum Flag_PI_gain = 0; Tmp1 =A1*A2-A3*A0 ; Tmp2 = 1 - beta2 ; Tmp3 = A02*A3 + A13 - 2*A0*A1*A2 ; if beta == 1 % Optimum controller-gain for reference-following if (Tmp1 #= 0) Kpj = A3 * 0.5 / Tmp1 ; else Kpj = Kmax/A0 ; % Use maximum gain Flag PI gain = 1; endif else % Optimum controller gain for disturbance rejection Sl = Tmp12 - A3 * Tinp2 * Tinp3 ; if Tmp3 ~= 0 Kpj = (Tmp1 - sqrt(S1)*sgn(Tmp1)) / (Tmp3*Tmp2) ; else Kpj = Kmax/A 0 ; % Use maximum gain Flag gain = 1; endif endif if (Kpj#A0 <= 0) # (abs (Kpj*A0) > Kmax) Kpj = Kmax / A0 ; % Gain correction Flag?PI-gain = 1 ; endif Tij = A 1 *Kp / (A0*Kp + 0.5 + 0. 5 *Kp2*A02Tmp2) ; Tdj = 0 ; Dead-time Compensation Controller (DTCC) CAA An overview of DTCC is shown in Fig. 5.

The control layer of D rCC (37) consists of a single PID controller in the velocity form and up to 10 local predictive functional controllers (PFCs) equipped with anti-windup protection and MV constrains restriction. The control layer is implemented by using the following procedure that is carried out at each time instant (k): % INPUTS: % Tsamp-sampling time of the PI control algorithm % r (k)-reference signal % y (k)-process output signal % limit umax, umin, duup, dudown-actuator limitation % controller_mode - The controller mode % u rnanual-controller output for the manual mode % Kp, Ti, Td-Current PID controller parameters % Wht1... Wht10 - weight coefficients for the local controllers'outputs % data... dst10, hr... haro-tuning parameters for the local predictive functional controllers % klar... klmlo, k2m1... k2m10, k1d1...k1d10, k2d1...k2d10, am1m...am1m10, % am2m1...am2m10, dm1...dm10, dd1...dd10-local model parameters % delt-dead zone % pfc_enable % ff - is feed-forward enabled % OUTPUTS: % uct (k)-controller output % INTERNAL SIGNALS (previous values stored for the next sample): e_l, uD_1, u_1 % upfc=0 ; % Calculation of Pm controller outputs % Calculation of control error e (k) = r (k)-y (k) ; % dead zone restriction if (e(k) >= - delt& en <= delt) e = 0 ; endif; % Derivative block calculation uD = uD_1 + Td / (alfa*Td + Tsamp) * (e(k)-e_1) + Tsamp /(alfa*Td + Tsamp)*(e(k)- uD) ; % gain block calculation uP = Kp * (uD - uD_1); % integral block calculation uI-Kp*Tsampl/Ti*uD ; % output calculation, upid (k) = u_1 + uD + uI ; % antiwindup reset if (upid(k) >= limit_umax # upid(k) <= limit_umin) upid(k) = u_1 + uP; endif % past values updating e_I = e (k) ; uD_1 = uD ; u 1 = upid (k) ; % Calculation of PFC local controllers'outputs for j = 1 to num_local_models % Model update routine for i = DTCC vector lengtlz to 2 ylm_memoiJ = y1m_memoi-1j; y2m_memoij=y2m_memoi-1j; y1d_memoij=y1d_memoi-1j; y2d_memoij=y2d_memoi-1j; end y1m_memo1j=am1mj*y1m_memo1j+klmj*(1-am1mj)*vector_u(1) ; y2m_memo1j=am2mj*y2m_memo1j+k2mj*(1-am2mj)*vector_u(1); y1d_memo1j=am1dj*y1d_memo1j+k1dj*(1-am1dj)*vector_v(1) ; y2d_memo1j=am2dj*y2d_memo1j+k2dj*(1-am2dj)*vector_v(1) ; if (controller_mode >= 2)&(pfc_enable ==10 lambda = exp(-3*Tsamp / dstj) ; % Reference trajectory coefficient if dmj > ddj difj = dmj-ddj, keyj = 1 ; else difj=ddj-dmj; =0 ; endif ifff == 0 ddj=0 ; eej=(1-lambdaj^hj)*(y1m_memo(dm+1,i)+y2m_memo(dm+1,i)); aaj=vector_rf(1+dmj+hj)-lambdaj^hj*vector_rf(1+dm)-(1-lambda ^hj)*y(k); bbj=0 ; ccj= (lambdaj^hj)*(y1m_memo1j+y2m_memo1j); ff=0 ; ggj=-(am1mj^hj)*y1m_memo1j+(am2mj^hyj)*y2m_memoij); hhj= (k1mj* (1-am1mj^hj)+k2mj*(1-am2mj^hj)) ; else ddj=(1-lambdaj^hj) * (y1d_memo(1+ddj, i)+ y2d_memo(1+ddj, i)) ; ee= (1-lambdaj^hj)*(y1m_memodmj+1,j+y2m_memodmj+1,j); if keyj ==1 aaj =vector-rf(1+dmj+hy)-lambdaj^hj* vector_rf(1+dmj)-(1-lambdaj^ hj)*y(k) ; bdj=(lambdaj^ hj)*((amldj^difj)*yld_memolj+kldj*(1-amldj^difj) *v (k)+(am2dj^difj)*y2d_memolj+k2dj * (1-am2dj^difj)*vj); ccj = (lambdaj^hj)*(ylm_memolj+y2m_memolj); ffj=-(amldj^(hj+difj)*yld_memolj+kldj * (1-amldj^(hj+difj)) *v (k) +am2d ^ ( hj+difj)*y2d_memolj+k2d*j(1-am2dj^(hj+difj))*v(k)); ggj=-((amlmj^hj)*ylm_memolj+(am2mj^hj)*y2m_memolj); hhj=(klmj*(1-amlmj^hj)+k2mj*(1-am2mj^hj)) ; elseif keyj ==0 aaj =vector_rf(ddj + hj+1)-lambdaj^hj*vector_rf(ddj+1)-(1-lambdaj^hj)*y(k); bbj=(lambdaj^ hj)*((amlmj^difj)*ylm_memolj+(am2mj^difj) * y2m_memolj); ccj=(lambdaj^hj)*(yld_memolj+y2d_memolj); ff =-(amlmj^(hj+difj)*ylm_memolj+am2mj^(hj+difj)*y2m_memolj); ggj=-((amldj^hj)*yld_memolj+kldj*(1-amldj^hj)*v(k)+(am2dj^hj )* y2d_memolj + k2dj*(1-am2dj^hj)*v(k)); hhj =-(lambdaj^hj)*(klmj *(1-amlmj^difj)+k2mj *(1-am2mj^difj))+klmj *(1- amlmj^(hj+difj)+k2mj * (1-am2mj^(hj+difj)); endif endif % Manipulated variable computation upfclj= (aaj+bbj+ccj+ddj+eej+ffj+ggj)/hhj; % Controller dead zone if abs(vector_rf(1)-y(k))<delt upfclj = upfc_lj ; endif upfc = upfc + Whtj*upfclj; % Blending of the local PF controllers outputs endif endfor % Controller mode switch & final output if (controller_mode >=2)&(pfc_enable ==1) uact (k) = upfc ; elseif controller_mode>=1 uact (k) = upid ; elseif controller mode == 0 uact(k) = u_manual; endif if uact(k) >lim_umax uact (k)= lim_umax; endif if uact (k) < lim_umin uact(k) = lim_umin ; endif % Constrains checking routine if lim umax < (vector_u1 + lim_duup*Tsamp) umx = lim_umax; else umx = (vector_u1+ lim_duup*Tsamp) ; endif if lim_umin > (vector_u1-lim)dudown*Tsamp) umn = lim_umin ; else umn=(vector_u1-lim_dudown*Tsamp) ; endif if uact (k) > umx uact (k) = umx ; endif if uact < umn uact (k) = umn ; endif % Controller dead zone if (abs(vector_rf1-y(k)) <delt)&(controller_mode #= 0) uact(k) = vector_u(1); endif % Update of MV and disturbance vectors forj = DTCC_vector_length to 2 j = vector u. I j = vector v endfor vector_u1=uact(k) ; vector_v1= v(k) ; The scheduling layer (36) of DTCC is implemented as fuzzy blending of the PID controller parameters or local PFC controller outputs according to the scheduling variable s and the membership functions of the local models. The following algorithm is implemented: % INPUTS : % CPM controller parameters matrices CPM = [CPl CP2... ?,] % where m is the number of local controllers, sorted by scheduling variable % CPj are local PID controller parameters, CPj=[Kpj Tij Tdj Whtlj... Wht10j]T % s-the current scheduling variable % OUTPUTS: % CP = [Kp Ti Td Whtlj...Wht10j]T - scheduled controller parameters at the current scheduling variable s % Find two neighbouring controllers j = 1 ; higher = sl ; while (s > s_higher) & (j<=m) j=j+ 1 ; if j <= m S_higher = Sj ; endif endwhile % Determine scheduled parameters ifj == 1 % s is lower than s CP = CPI ; elseifj = m + 1 % s is higher than Sm CP = CPm; else % S is between S1 and Sm, actual fuzzy scheduling tmp=(S-Sj-1)/(Sj-Sj-1) ; Kp = Kpj-1 + tmp * (Kpj-Kpj-1); Td Tdj-l + tmp * (Tdj-Tdj-,) ; Ti=Tij-1+tmp*(Tij-Tij-1) ; Wht1=Whi1j-1+tmp*(Wht1,j-Wht1j-1) ; <BR> <BR> <BR> <BR> <BR> <BR> wht10=Whi10j-1+tmp*(wht10,j-Wht10j-1); endif The tuning layer of DTCC (33) is based on tuning tables for parameters of the predictive controller and the PID controller and classical approaches for PID controller tuning. There are two manners of obtaining PID settings. The first one is based on formulas approximating the results of a number of simulations with various stable plant models performed in advance. The second manner comprises the approach based on the internal model scheme. PFC tuning parameters are calculated based on formulas obtained with simulations (in the way as the first PID tuning manner). The following procedure is used for tuning of a single local PID and PFC controller (35) from their corresponding local model (34): % INPUT PARAMETERS : % Tsamp-sampling time % LMj=[a1j, a2j, b1j, b2j, duj, sj]- local model parameters, delay and scheduling variable % position, extracted from the MFM (second order models: a2j = b2j = 0) % speed factor-correction coefficient for PID controller tuning parameters % no_del_gain-PID controller gain when no delayed model is available %. tune manner-two different formulas for calculation of PID tuning parameters are provided % model_orders - odrer of local models % OUTPUT PARAMETERS : % CPj - PID controller parameters, Cl 'j = [Kpj Tij Tdj Wht1...Wht10] % CPaj-PFC controller parameters, CPaj = [dstj hj k1mj k2mj k1dj k2dj am1mj am2mj dmj ddj] % if (nzodel orders == 2) & (real roots), i = (alj + sqrt (a21j--4*a2j))/(-2) ; g= a1j - i; f= (b2j-g*b1j)/(i-g) ; h= (blj*i-b2j)/(i-g) ; f_f=(c2j-g*c1j)/(i-g); j *i-c2j) l (i-g) ;. dml = din *Tsamp; ddl dd*Tsamp ; am1mj=-i; am2mj=-g; k1mj=h/(1-am1mj); k2mj=f/(1-am2mj) ; k1dj=h_h/(1-am1mj) ; k2dj=f_f/(1-am2mj) ; T1mj=-Tsamp/log(am1mj); T2mj=-Tsamp/log(am2mj); Kj=k1mj+k2mj; Tj=T1mj+T2mj; elseif model_orders ==1 % first order model i = (all + sqrt(a21j-4*a2j))/(-2); g= a1j-i; f= (b2j-g*blj)/(i-g) ; h=(b1j*i-b2j)/(i-g); f_f=(c2j-g*c1j)/(i-g); h_h=(c1j*i-c2j)/(i-g) ; dm1j=dmj*Tsamp ; dd1j=ddj*Tsamp; am1mj=-i; am2mj=-g; k1mj=h/@ amlmj) ; k2mj=f/(1-am2mj) ; kldj=h_h/(1-am1mj) ; k2dj=f_f/(1-am2mj) ; T1mj=-Tsamp/log(am1mj); T2mj=-Tsamp/log(am2mj) ; Kj=kl mj+k27nj ; Tj=T1mj+T2mj; endif if dmlj=D tune_manner=1 ; endif kappa= dm1j/Tj; if kappa < = 0. 4 esp=0.9*dm1j; else esp=dnalj ; endif if drnlj==0, Kpj= speedfactor*no_del_gain/Kj ; if tune_manner==0 S0j=1. 0/(Kj*dm1j)*(0.287+0.21*Tj/dm1j+0.385*exp(-1.57*Tj/dm1j)) ; S1j=1/Kj*(-0.186+0.784*Tj/dm1j+0.435*exp(-1.4*Tj/dm1j)); S2j=dm1j/Kj*(0.076+0.207*Tj/dm1j+0.024*exp(-1.25*Tj/dm1j)) ; Tij=S1j/S0j; Tdj=S2jlSlj,- else tune-manner Tij=Tj+dm1j/2 ; Tdj=Til6,- endif else if tune-manner 0 S0j=1. 0/(Kj*dm1j)*(0.287+0.21*Tj/dm1j+0.385*exp(-1.57*Tj/dm1j)) ; S1j=1/Kj*(-0.186+0.784*Tj/dm1j+0.435*exp(-1.4*Tj/dm1j)) ; S2j=dm1j/Kj*(0.076+0.207*Tj/dm1j+0.024*exp(-1.25*Tj/dm1j)) ; Tij=S1j/S0j; Tdj=S2j/S1j; Kpj=0.43*speed_factor*S1j; elseif tune_manner ==1 Kpj=0.43*speed_factor*(Tj+dm1j/2)/Kj/(esp+dm1j/2); Tij=Tj+dm1j/2; Tdj=(Tj*dm1j/2)/(Tj+dm1j/2) ; endif endif dstj=0. 5*(1.2463*dm1j+0.18822*Tj) ; hj=ceil (5*0.2231*Tj/Tsamp) ; for i=l : 10, if i==j Whti=l ; <BR> <BR> <BR> <BR> else<BR> <BR> <BR> <BR> <BR> else<BR> Whti=0, endif endfor for i=1 to DTCC_vector_length ylm_memog=klmj*vector_u ; y2m_memoij=k2mj*vector_u; y1d_memoij=k1dj*vector_v; y2d_memoij=k2dj*vector-v; endfor Rule-based Neural Controller (RBNC) CAA An overview of RBNC is shown in Fig. 6.

The control layer of the RBNC (38) consists of a single PID controller in the velocity form equipped with anti-windup protection, and is specified by the following procedure that is executed at each discrete time step k : % INPUTS: % Ts-sampling time of the PI control algorithm % r (k)-reference signal % y (k)-process output signal % limit_umax, umin, slew_max, min - actutor limitation % cofatroller naode-The controller mode % u_manual - controller output for the manual mode % Kp, Ti, Td - Current PID controller parameters % alfa-tracking constant % beta-jump constant % OUTPUTS: % u (k)-controller output without constraints % uout (k)-controller output with constraints % INTERNAL SIGNALS (previous values stored for the next sample): % y_1, y_2, u_1, e_1, eu-1, r_1 % RBNC CONTROLLER ALGORITHM % Calculation of PID terms e (k) = r (k) -y (k); Tt = alfa*Ti ; % tracking time constant cO =Kp*(1+Ts/Ti) ; % e (k) multiplier cl =-Kp ; % e (k-1) multiplier f0 = -Kp*Td/(Ts^2); % y (k) multiplier fl =-2*f0 ; % y (k-1) multiplier f2 = fO ; % y (k-2) multiplier jump =Kp*(beta-1); % jump factor enters only when set point changed dole = c0*e + cl *e_1 ; % dP + dI terms using process error dDy =f0*y + f1*y_1 + f2*y_2 ; % dD term using process output du = dPIe + dDy + jurnp* (r-r_1) ; % total incremental control signal including jump term % Calculate unlimited controller output deltau = du + (Ts/Tt) *eu_1 ; % basic dP+dI+dD algorithm with anti-windup u(k)=u_1 + deltau; % controller output without constraints % Calculate the new value of u (k) (limited controller output) if controller-mode = = Manual u (k) = u_manual; deltau =u(k) - u_1 ; endif % Controller output limitation if deltau >= slew_max, % Controller slew rate limit du = slew_max ; % upper bound else if deltau <= slew_min du = slew_min ; % lower bound else dit = deltau; % in linear regime endif % Calculate the Total Controller output with constraints MM = u_1 + du; % Total controller output if uin >= limit umax % Controller output limit uout (k) = limit_umax; % upper bound else if uin <= limit_umin uou (k) t = limit umin ; % lower bound else uou (k) t = MM ; end % Update internal signals for next sample eu=uout-ui7l ; eu_1 = eu; u_1 = uout; e_1 = e; r_1 = r; y_2 = y_1; y_1=y; The scheduling layer of the RBNC (39) (i. e. the RBNC Supervisor) performs switching of controller parameters in accordance with linguistic rules and the values of the process output and target set-point. The kernel of the RBNC Supervisor involves two 2-4-1 artificial neural networks (ANN) (44). The neurons of all but the output neuron possess sigmoidal compression.

Each ANN is trained using a set of 11 coded linguistic rules. The object of the zone verification block (43) is to establish the indices of the current and target operating zones. An operating zone is defined as a region about each point where a local process model has been identified.

The boundaries of each zone are taken as the mid-points between contiguous points at which a local process model has been identified. These two indices form the inputs to the ANNs in the RBNC Supervisor. The RBNC Supervisor decides which candidate controller to use while the process output is in transition from one operating point to another.

% INPUTS : % Ts-sampling time of the PI control algorithm % controller mode-The controller mode % y (k)-process output signal % sp-target set point % ZLM-is a matrix of zone limits ZLM= [ZLj, Zl2, ZLnJ for each operating zone % where ZLi contains values that are pertinent to each operating zone % CPM-controller parameters matrix CPM=[CP1 CP2 ... CPm] % where m is the number of local controllers, sorted by scheduling variable % CPj are local PID controller parameters, CPj Kpj Tij Tdj sj] T % s-the current scheduling variable % LMj = [a1j, a2j, b1j, b2j, duj, sj]- local model parameters, delay and scheduling variable %-position, extracted from the : MFM---- (second order models: a2j=b2j=0) % model order-model order of the MFMs % OUTPUTS: % zone-zone index uniquely specifying the current operating zone and corresponding optimum % controller parameters % RBNC NEURAL RULE-BASED SUPERVISOR ALGORITHM Step 1: determine zone index of y Step 2: determine zone index of sp Step 3: determine whether to use increscent or decrescent switching logic that guarantees asymptotic stability between contiguous operating Zones. The procedure makes use of the Jury stability test.

The increscent switching logic has the following steps: 1) if the target set-point is in the same Zone as the current operating Zone then do not alter the Controller parameters 2) if the target set-point is in a Zone higher than the present operating Zone then increment the Zone index successively (and apply the corresponding Controller parameters) until the target Zone is reached.

3) if the target set-point is in a Zone lower than the present operating Zone then starting from the current Zone parameters decrement the Zone index successively (and apply the corresponding Controller parameters) until the target Zone is reached. while the decrescent switching logic has the following steps: 1) if the target set-point is in the same Zone as the current operating Zone then do not alter the., Controller parameters 2) if the target set-point is in a Zone higher than the present operating Zone then starting from the current Zone parameters, increment the Zone index successively (and apply the corresponding Controller parameters) until the target Zone is reached.

3) if the target set-point is in a Zone lower than the present operating Zone then decrement the Zone index successively (and apply the corresponding Controller parameters) until the target Zone is reached.

The tuning layer of the RBNC (40) uses a stochastic procedure, and in particular Simulated Annealing based on the Metropolis Algorithm with an embedded Jury stability test (45), to compute the globally optimum parameters of the PID controller corresponding-to each operating zone.

% INPUTS : % Ts-sampling time of the PI control algorithm % controller mode-The controller mode % LMj =[a1j, a2j, b1j, b2j, duj, sj]- local model parameters, delay and scheduling variable % position, extracted from the MFM (second order models: a2j = b2j = °) % OUTPUTS: % CP * _ (Kp Ti Td s JT-optimum controller parameters at the current scheduling variable s % beta--jump constant % RBNC TUNING ALGORITHM Step 1: determine the time to reach 63% of final value of step response Step 2: using the Persson and Astrom technique compute estimate of the parameters of the sub-optimum PID or PI controller Step 3: test stability of the closed system with the corresponding LMi using the Jury test Step 4: use the Simulated Annealing search algorithm with embedded Jury test to compute the globally optimum controller parameters 6. Control Performance Monitor (CPM) CPM analyses the buffer of recent real-time signals, prepared by SPA, and estimates features of closed-loop control response when conditions are appropriate. It is invoked periodically. Since it is a computationally intensive task, it runs as a low-priority task in a multitasking operating system. It consists of three submodules : the Buffer Founder (BF), the Situation. Classifier (SC) (57), and the Performance Estimator (PE) (58).

When CPM is invoked, BF forms new dynamic buffers from the current buffers according to the process conditions so that only the correct data are sent to SC. SC scans the buffer for the last recognisable event, ; such as a step change of the reference, a step change of the measured disturbance, an unmeasured disturbance, a steady state, or oscillation. If such an event is detected, the corresponding section of the buffer is sent to PE, otherwise the execution is terminated.

For the detected events, PE may extract the following features: overshoot, settling time, rise time, oscillation, tracking error, regulation error.

The following procedure is carried out whenever CPM is invoked, as illustrated in Fig. 7: Copy of real data buffer (47) Steady State check (48) Data Filtration (49) Reference and Measured Disturbance analysis (50) Oscillation check (51) Noise/Signal Estimation (52) Transient process settlement check (53) Steady State check (54) Features Estimation (55) Performance Index calculation (55) Send Results to Operator Supervisor and Human Machine Interface (56) Copy of real data buffer : This module calls functions, which copy a part of the real-time data buffer for background processing. The so-formed arrays r, y, u, and v are the main CPM input parameters.

Steady state check : By using peak-to-peak and linear fit criteria, this module checks whether all signals r, y, u, and v in the active buffer have reached steady state (SS). If so, there is nothing to evaluate and the processing is terminated; if not, the algorithm continues with the main part of CPM., Filtration : This module performs filtering of the input y and v buffers (if requested). Three filtering methods may be used: moving average, exponential, or fourth differences.

Reference and Measured Disturbance analysis : The main task of this module is checking the input r and v buffers for a change. If a step-wise change of either signal is recognized, the CPM continues. with oscillation checking routine. If there is an unrecognised. event, the CPM sends"UNRECOGNIZED EVENT"message to the operator panel and terminates any further computations.

Oscillation Check : This module checks the CV buffer for oscillation by comparing the decay ratio to its threshold.

The main task of this checking routine is to detect oscillation as soon as possible. This procedure is intended for avoiding accidents or undesired CV behaviour. If oscillation is detected, theCPM terminates and sends a flag to the Operation-Supervisor-for switching to the safe mode.

Noise/Signal Estimation : The main task of this module is calculating the Noise/Signal ratio of the controlled variable. An estimate of the noise variance is calculated on the base of a segment of the CV buffer before the registered change. The requested minimum data can be specified. The amplitude of the signal is calculated as an absolute value of the difference between the past and the new steady state value of MD (if there is a MD change) or as the absolute value of the reference change (if there is a SP change). The calculated ratio is compared to a configured threshold. 0 C Noise/Signal = (5 Y or I Lori ove If the threshold is not exceeded, the CPM processing of the buffer continues, otherwise it is terminated.

Transient process settlement check.- Check if the transient process has settled to a new steady state. If not, the processing is terminated.

Steady State check.

Check if the process was in steady state before the recognized event. If not, the processing is terminated.

Send Results to Operator Supervisor and Huma7i Machine Interface (HMI).

This module comprises functions that send messages to HMI ; and the flag to that indicates a request for switching to the safe mode to the Operation Supervisor.

Features Estimation : This procedure computes the transient process features like overshoot, settling time, rise time, oscillation, tracking error, and regulation error when the conditions are appropriate. <BR> <BR> <BR> <BR> <BR> <BR> <BR> appropriate_conditions = #oscillation_detected AND noise_signal < noise_signal_threshold AND transient processfinish AND steady_state before change ; if (reference or measured disturbazzce change AND appropriate_ conditions) Calculate overshoot, settling time, rise time, oscillation, tracking error, and regulation error; Performance Index calculation; endif Performance Index calculation Based on the computed process features (overshoot (59), settling time (60), decay (61) ) a generalized criterion for the quality of the transient process is computed by using fuzzy rules in Fig. 8. This criterion is called Performance Index (PI) (62). For each process feature a set of linguistic variables is introduced-poor, good, and excellent. A rule table based on expert knowledge is created (Fig. 9) which maps the PI values and all combinations of feature values.

PI has five linguistic variables-poor, satisfactory, good, very good, excellent, corresponding to the following values ; 0.0, 0.25, 0.5, 0. 75,1. 0. From the received value of PI using the defuzzyfication procedure the transient response quality is evaluated. PI lower than 0.25 triggers a request for switching to the Safe mode, while PI between 0.25 and 0.5 results in a sending warning message to the HMI BEST MODES FOR CARRYING OUT THE INVENTION The run-time module (RTM) of the proposed controller, including all described agents, is intended for implementation on a capable programmable logic controller (PLC) or an embedded controller. One of the test platforms consists of a Mitsubishi Al S series PLC with an INEA IDR SPAC20 coprocessor, based on the TMS320C32 digital signal processor of Texas Instruments, running at a clock frequency of 40 MHz and with 2MB of RAM. The initial configuration of the RTM is performed using a Configuration Tool (CT), running on a personal computer, that is used for transferring the RTM parameters to the PLC and includes a"wizard"that guides the engineer through the configuration of initial parameters.

The operation of the controller is illustrated with an application for pH control of a simulated neutralisation process in Fig. 10, described in detail in (Henson & Seborg 1994, Kavsek- Biasizzo & al. 1997, Gerksic & al. 2000). An acid stream (63), a buffer stream (64) and a base stream (65) are mixed in a tank. The acid and base streams are equipped with flow control valves. The pH of the mixture is measured with a sensor located downstream (66). The effluent pH is the controlled variable y, and the manipulated variable u is the flow of the base stream Q3.

The static characteristic of the process is highly non-linear and its open-loop gain changes by the factor 8, so it is very difficult for control with a conventional PID controller.

The non-linear model simulated on a personal computer using Matlab/SimulinkofThe MathWorks, Inc. , was used in the role of the real process, while the controller on the PLC was using a model obtained by online learning. The operating range between pH values 6 and 8 was covered with five local models placed, at positions 6,7, 7.15, 7.4, and 8, with respect to the known titration curve. s (k) = 0.3 y (k) + 0.7 r (k) was used as the scheduling variable. Each local model was trained with online learning through an open-loop experiment consisting of three step-changes of u of small amplitude about the operating point, and the local controllers were tuned automatically from the local models. In the upper part of Fig. 11, the process output (67) follows the reference signal (68) reasonably well under FGSC control, while the performance of a fixed PI controller in the lower part (69,70) is much worse. The PI controller, which is also used as the Safe mode, is tuned for stable performance in the high-gain operating area around pH 8, however its operation in the low-gain region below pH 7 is sluggish.

REFERENCES : Astrom, K. J., Wittenmark, B. , 1995, Adaptive control. Addison-Wesley Publishing Company, Inc. , Reading, Massachusetts.

Bequette, B. W. , 1991, Non-linear Control of Chemical Processes: A Review. Ind. Eng. Chem.

Res. , Vol. 30, pp. 1391-1413.

Gerksic, S. , Juricic, D., Strmcnik, S. , Matko, D. , 2000, Wiener model based non-linear predictive control. International Journal of Systems Science, Vol. 31, No. 2, pp. 189-202.

Henson, M. A., Seborg, D. E. , 1994, Adaptive non-linear control of a pH neutralisation process IEEE Transactions on Control Systems Technology, Vol. 2, pp. 169-182.

Henson, M. A. , Seborg, D. E. , 1997, Non-linear Process Control. Prentice-Hall PTR, Upper Saddle River, NJ..

Kavsek-Biasizzo, K. , Skrjanc, I., Matko, D. , 1997, Fuzzy predictive control of highly non-linear pH process. Computers & Chemical Engineering, Vol. 21, pp. S613-S618.

Ljung, L. , 1987, System Identification. Prentice Hall, Englewood Cliffs, NJ.

Murray-Smith, R. , Johansen, T. A. , 1997, Multiple Model Approaches to Modelling and Control. Taylor & Francis, London.

Qin, S. J. , Badgwell, T. J. , 1998, An Overview of Non-linear Model Predictive Control Applications. Proc. Non-linear Model Predictive Control Workshop Ancona'98, Ancona, Birkhäuser Verlag AG, Basel (2000).

Sjoberg, J. , Zhang, Q. , Ljung, L. , Benveniste, A.,. Deylon, B., Glorennec, P. -Y. , Hjalmarsson,<BR> H. , Juditsky, A., l995j Non-linear Black-Box Modeling in System Identification: a Unified_Oyerview. Automatica 31 (12), pp. 1691-1724.

Vrancic, D. , Strmcnik, S. , Juricic, D., 2001, A magnitude optimum multiple integration tuning method for filtered PID controller. Automatica, vol. 37, pp. 1473-1479.