- Email: [email protected]
A survey of grinding circuit control methods: from decentralized PID controllers to multivariable predictive controllers
Powder Technology 108 Ž2000. 103–115 www.elsevier.comrlocaterpowtec
A survey of grinding circuit control methods: from decentralized PID controllers to multivariable predictive controllers Andre´ Pomerleau a b
a,)
´ Gagnon , Daniel Hodouin b, Andre´ Desbiens a , Eric
b
Department of Electrical Engineering, LaÕal UniÕersity, Quebec, Canada G1K 7P4 Department of Mining and Metallurgy, LaÕal UniÕersity, Quebec, Canada G1K 7P4 Accepted 20 September 1999
Abstract A conventional grinding circuit consisting of one open-loop rod mill and one closed-loop ball mill is essentially a two-input= twooutput system, assuming that the classifier pump box level is controlled by a local loop. The inputs are the ore and water feed rates and the outputs are the product fineness and the circulating load. The design problem is to find a control algorithm and a tuning procedure which satisfy specified servo and regulatory robust performances. A first approach is to use decentralized PID controllers and systematic tuning methods which take into account loop interactions. Another technique consists of adding decouplers or pseudo-decouplers to the decentralized controllers. Finally, the design of a fully multivariable controller is a possible option. To face the problem of performance robustness related to change of process dynamics, two options are studied. A design criterion involving the minimization of a penalized quadratic function on a future trajectory can be used. A second alternative is to track process dynamics changes using adaptive process modelling. The paper will present a comparison of these various strategies, for a simulated grinding circuit. A benchmark test, involving a sequence of disturbances Žgrindability, feed size distribution, change of cyclone number . . . . and setpoint changes, is used to compare the performances of the controllers. q 2000 Elsevier Science S.A. All rights reserved. Keywords: Grinding circuit; Distributed PID; Internal model control; Multivariable predictive control; Distributed adaptive predictive control
1. Introduction The objective of this paper is to study the performances of various control algorithms acting on a grinding circuit in face of different disturbances and setpoint changes. Even though numerous articles w1–3x have studied this control problem, it is hard to compare the results since the benchmarks are different. In this article, the characteristics of the process, from a control point of view are first founded. From there, various control algorithms will be applied and general conclusions on the performances will be drawn. The algorithms studied are: –
Decentralized PID controllers which are characterized by this fixed structure. Algebraic internal model with explicit decouplers. They are characterized by a model-based structure. Full multivariable predictive controllers. Distributed adaptive predictive controllers. These are
– – – )
Corresponding author.
characterized by variable parameters as function of the process model variation.
2. Process description The grinding process is a size reduction operation used in the mineral industry to liberate the valuable minerals from the discardable gangue. Fig. 1 shows the flowsheet of the plant: it is a standard two-stage grinding circuit which consists of a rod mill in open circuit and a ball mill in closed circuit with a hydrocyclone classifier. Two local controllers are assumed to be part of the process. Sump level is maintained at a constant value using a variable speed pump controlled by a PID. The percentage of solids in the feed is kept around a desired set-point using a feed water addition rate proportional to the fresh ore feed rate. Grinding typically represents around 50% of the total expenditure of the concentrator operation. This step prepares the ore for the concentration step such as flotation or magnetic separation. Therefore, it determines the feed characteristics to these subsequent separation stages and
0032-5910r00r$ - see front matter q 2000 Elsevier Science S.A. All rights reserved. PII: S 0 0 3 2 - 5 9 1 0 Ž 9 9 . 0 0 2 0 7 - 7
104
A. Pomerleau et al.r Powder Technology 108 (2000) 103–115
Fig. 1. Grinding circuit plant flowsheet.
influences the efficiency of the valuable minerals recovery process. The true objective of grinding is to obtain a proper liberation of the minerals coupled with a proper particle size distribution in order to maximize recovery and concentrate value without overgrinding, while keeping specific comminution energy and reagent consumption as low as possible, taking into account the prevailing economic conditions. These objectives being highly complex, simpler ones must be formulated. Generally, as the liberation degree cannot be measured on line, the production objective is to maintain at a set-point value the percentage of the product particles smaller than the reference size Žfor instance 70% smaller than 74 mm..
In parallel, the circulating load Žtons per hour of solids entering the ball mill. is to be maximized. High recirculation at a fixed product size means lower energy consumption. Excessive circulating load can result in classifier, mill or pump overload conditions. The selected controlled variables are therefore the size distribution and the circulating load. The fresh ore feed rate and the water addition to the sump box are the available manipulated variables. Since industrial design and testing are costly, the dynamic grinding circuit simulator DYNAFRAG is used w4x. It simulates the mineral processing plant with appropriate dynamic phenomenological models Žset of mass balance equations for the pulp, the water and the particles of different hardness.. The detailed mathematical formulation of the models involved in the simulator is given elsewhere w5x. From 100 to 150 simultaneous differential equations are to be solved. With DYNAFRAG, it is also possible to simulate disturbances such as sensor and actuator noises, natural variations in ore characteristics Žparticularly hardness changes., changes in circuit performance and mechanical upsets Žblockages in pumps and classifiers, etc... The simulator operates in a multi-computer environment. The ‘‘process’’ and the ‘‘control room’’ are physically separated. They are both simulated on their own micro-computer and are interconnected through digital-to-analog and analog-to-digital interfaces. This multi-computer architecture resembles industrial reality.
Fig. 2. Open-loop step response.
A. Pomerleau et al.r Powder Technology 108 (2000) 103–115
3. Process model Fig. 2 illustrates the open-loop responses of DYNAFRAG Žsampling period: 0.5 min.. Fig. 2a and b shows the size distribution ŽSD: percentage passing 74 mm. and the circulating load ŽCL. following an 8.8 m3rh step on the pump box water ŽPBW. at time t s 10 min. The responses to a 6 trh step on the rod mill feed ŽRMF. at time t s 10 min are given in Fig. 2c and d. Fig. 3 shows the same responses with white noise added to the measured variables. The same level of noise is kept during the closed-loop simulations. The linear process model identified at the operating point is: 0.6925ey5 s
CL s RMF CL s PBW
1 q 18 s 0.475 Ž 1 y 3.5s .
Ž 1 q 4.5s . Ž 1 q s .
SD s
y1.63 Ž 1 y 2.5s . ey8 s
RMF
Ž 1 q 28 s . Ž 1 q s . SD 0.393 Ž 1 q 92 s . eys s PBW Ž 1 q 21 s . Ž 1 q 2 s . The time constants are expressed in minutes. From these transfer functions, one concludes that a grinding
105
circuit is a highly coupled process where PBW addition has a very fast short-term effect on the size distribution, but where the RMF has the dominating long-term effect.
4. Selected control algorithms There are two main approaches to select the control loop architecture for multivariable processes: – decentralized Single-Input Single-Output ŽSISO. loops – multivariable controllers, where the latter can also be subdivided in two depending if the decoupler is explicitly expressed or not. Fig. 4 shows a multivariable process controlled by distributed SISO controllers. In this structure, the proper inputroutput pairing has to be chosen and the tuning of one controller depends on the tuning of the second controller. In a structure where the decoupler is explicitly expressed as shown in Fig. 5, a decoupler structure has to be chosen and the facility of decoupling depends on the choice of the coupling used. The controller tuning of the SISO controllers is done on the combined system including the process and the decoupler. Finally, in a fully multivariable controller as shown in Fig. 6, the designer does not need to choose the pairing, but the specifications have to be feasible.
Fig. 3. Open-loop step response Žnoise environment..
106
A. Pomerleau et al.r Powder Technology 108 (2000) 103–115
Fig. 4. Distributed SISO controllers.
The first studied algorithm is the standard decentralized PID controller which is usually the reference to compare different algorithms. PID controllers present some drawbacks. First, their performances are limited by the properties of the fixed structure. Second, the inputroutput pairing problem has to be solved for multivariable systems. Finally, an efficient tuning method has to be used to obtain the desired dynamics of the controlled loops. The second studied algorithm consists of distributed internal model controllers applied to the process preceded by a partial decoupler. The decoupler partially solves the pairing problem and facilitates the tuning. The internal model control structure which is model-based eliminates the problems of a fixed structure. The third studied algorithm is a full multivariable optimal predictive controller. It totally solves the pairing problem and takes into account the compromise between the performances and the manipulated variables energy through a quadratic criterion. Finally, a distributed adaptive predictive controller is used. Compared to the other algorithms, the adaptive part permits to take into account the nonlinearities of the
process and the time evolution of the parameters. It also solves the tuning problems of the distributed control algorithm. The time specification on the circulating load has been arbitrarily fixed to 25 min. This corresponds approximatively to a dynamic similar to the open-loop function between the CL and RMF. The specification on the size distribution has been fixed to 5 min. Here, since the open-loop transfer functions relating the size distribution to the RMF or to the PBW have zeroes, it is more difficult to relate the specification to the open-loop time constants.
5. Pairing of controlled and manipulated variables Distributed control consists of using n SISO to control an n = n process as shown in Fig. 4. The main difficulties of distributed control lie in the choice of the proper input–output pairing and in the fact that the tuning of each SISO controller depends on the tuning parameters of the second SISO controller. When manipulated variables are not properly selected, interactions between controlled and
Fig. 5. Distributed SISO controllers with decoupler.
A. Pomerleau et al.r Powder Technology 108 (2000) 103–115
107
Fig. 6. Full multivariable controller.
manipulated variables can result in undesirable control loop interactions, leading to poor control performances. The well-known Bristol’s relative gain array ŽRGA. w6x being based only on steady-state information does not consider the process dynamics and is inappropriate for many processes. This is especially the case when dynamics in transversal transfer functions are fast compared to direct transfer functions and present smaller static gain than
direct transfer functions. The presence of zeros in the transfer functions also degrades the performance of RGA. Consequently, a dynamic method must be employed. The basic and simple idea of the RDG ŽRelative Dynamic Gain. method is to extend the calculation of the Bristol’s relative gain for all frequencies. The method used here is the generalized dynamic relative gain w7x. For the process under study, the RDG frequency responses are
Fig. 7. RDG frequency response.
108
A. Pomerleau et al.r Powder Technology 108 (2000) 103–115
given in Fig. 7. For the dynamics required by the specification, the output Y1 , is paired with U1.
The system can be drawn as depicted in Fig. 8. For each closed loop, the control objectives are: Ž1. no permanent error and Ž2. a second-order set point tracking response with specified dynamics. Thus, the corresponding closedloop transfer functions can be selected as follows:
6. Algorithm 1: Decentralized extended PIDs Y1 Ž s . To tune the controllers, a frequency-based method, in which the open-loop characteristics are deduced from the closed-loop dynamics is used w8x. First, for tuning a controller, it is imperative to know the dynamics of the process to be controlled. With monovariable processes, the method consists of opening the loop and evaluating the transfer function seen by the controller. For multivariable processes, the loop with the controller under consideration is opened while all the other loops are kept closed. For a 2 = 2 process, the transfer function seen by the controller of the first loop, GC1Ž s ., is: G 1 Ž s . s G 11 Ž s . y
G 12 Ž s . G 21 Ž s . GC 2 Ž s . 1 q GC 2 Ž s . G 22 Ž s .
.
The transfer function seen by the second controller GC2 Ž s . is: G 2 Ž s . s G 22 Ž s . y
G 12 Ž s . G 21 Ž s . GC1 Ž s . 1 q GC1 Ž s . G 11 Ž s .
.
As can be seen, there exists an interdependence between GC1Ž s . and GC2 Ž s .. Indeed, their individual tuning depends on the choice of the other controller. Consequently, to evaluate the unknowns GC1Ž s . and GC2 Ž s ., further information is necessary. This supplementary information is given by the closed-loop specifications which can be translated into open-loop char-acteristics.
Ž 1 q t 11 s . Ž 1 q t 21 s .
Z1 Ž s . Y2 Ž s . Z2 Ž s .
1 y t 01 s
s
1 y t 02 s
s
Ž 1 q t 12 s . Ž 1 q t 22 s .
The open-loop characteristics deduced from the tracking closed-loop specifications are: G 0 L 1 Ž s . s GC1 Ž s . G 1 Ž s . 1 y t 01 s
s
s Ž t 11 q t 21 q t 01 q t 11t 21 s . G 0 L 2 Ž s . s GC 2 Ž s . G 2 Ž s . 1 y t 02 s
s
s Ž t 12 q t 22 q t 02 q t 12t 22 s . The transfer functions GC1Ž s . and GC2 Ž s . of the controllers are obtained by solving the previous equations. Due to their quadratic form, it leads to two Bode plots for both controllers. Knowing the required sign of the controller gain and the necessary presence of an integrator, the appropriate solution can be selected without any difficulty. The last step in tuning the controllers consists of defining the transfer functions GCP1Ž s . and GCP2 Ž s . that best approximate the frequency response GC1Ž j v . and GC2 Ž j v ., while ensuring a stable closed-loop system. The structure of each controller is imposed and an evaluation of the
Fig. 8. Distributed SISO controllers equivalent system.
A. Pomerleau et al.r Powder Technology 108 (2000) 103–115
parameters is performed using a constrained objective function. The following structure is used: GCP i Ž s . s
K C Ž 1 q Ti s . Ž 1 q Td1 s q Td 2 s 2 . Ti s Ž 1 q Tf 1 s q Tf 2 s 2 .
.
2
GC1 Ž s . s GC 2 Ž s . s
1.28 Ž 1 q 85.4 s . Ž 1 q 19.15s q 73.4 s . 85.4 s Ž 1 q 17.81 s q 90.6 s . 0.32 Ž 1 q 2.4 s . Ž 1 q 4.77s q 9.06 s 2 . 2.4 s Ž 1 q 6.4 s q 7.8 s 2 .
7. Control algorithm 2: algebraic multivariable control Multivariable algebraic control can be done either with a full multivariable controller or by using a decoupler in cascade with monovariable controllers. The latter one, as shown in Fig. 9, is used here. The methodology consists of: – evaluating the decoupler, – tuning the monovariable controllers on the model including the multivariable process model and the decoupler. For many processes Že.g., grinding circuit., it is impossible to realize a perfect decoupler. This happens with non-minimal phase processes or with systems with longer delays in the direct transfer functions. For the system under study, the decouplers are: D 12 Ž s . s y s
G 22 Ž s .
1.633 Ž 1 y 2.5s . Ž 1 q 2 s . Ž 1 q 21 s .
D 21 Ž s . s y
Ž 1 q 92 s .
y7 s
e
G 21 Ž s . G 11 Ž s .
0.475 s
0.6925 ey5 s
Ž 1 q 18 s . Ž 1 y 2.5s . Ž 1 q 2 s . Ž 1 q 21s . ey7 s q 1.974 Ž 1 q s . Ž 1 q 28 s . Ž 1 q 92 s . Ž 1 y 3.5s . = Ž 1 q 4.5s . Ž 1 q s . 0.393 Ž 1 q 92 s . eys G2 Ž s . s Ž 1 q 2 s . Ž 1 q 21s . Ž 1 y 3.5s . Ž 1 q 18 s . q 1.120 Ž 1 q 4.5s . Ž 1 q s . Ž 1 y 2.5s . ey8 s = Ž 1 q s . Ž 1 q 28 s . These transfer functions can be approximated by: G1Ž s . s G2 Ž s . s
2.666ey5 s
Ž 1 q 85s . 1.513 Ž 1 q 5s . eys Ž 1 q 3s .
2
The distributed internal model controllers are then tuned into these transfer functions. 8. Control algorithm 3: multivariable state-space internal model predictive control
G 12 Ž s .
0.393 Ž 1 q 28 s .
Under the assumption that the system is satisfactorily decoupled for neglecting the remaining coupling for tuning consideration, the transfer functions for the system under study are given by G1Ž s . s
The resulting controllers are:
109
Ž 1 y 3.5s . Ž 1 q 18 s . 0.69 Ž 1 q 4.5s . Ž 1 q s .
D 12 Ž s . realizes a perfect decoupler, while D 21Ž s . is a pseudo-decoupler.
Here the controller is designed to control the process model and not the process itself. The difference between the model and the process outputs, respectively, y and y P submitted to the same control action u, is then interpreted as an overall estimation d of the disturbance d. This estimated value embeds in fact both disturbances and modelling errors. The value d is used to modify the setpoint Žreference value r .. This biased setpoint is the desired value yd that must be tracked by the process model. The controller calculates u in such a way to cancel the difference e between yd and the model output ym .
Fig. 9. Internal model controller with explicit decoupler.
A. Pomerleau et al.r Powder Technology 108 (2000) 103–115
110
Fig. 10 summarizes the basic principles of this control architecture. It is generally named Internal Model Control ŽIMC.. At time i, the sequence of control actions to be applied in the future Ž u i , u iq1 , . . . u iqHy1 ., where u k is the vector of control actions at time k and H the control horizon, is selected such that the following quadratic performance criterion is minimized: iqH
Ji Ž u i . . . u iqHy1 . s
e kT Q e k
Ý ksiq1
iqH
q
Ý Ž u ky 1 y u ky2 . T ksiq1
=R Ž u ky 1 y u ky2 . In the equation, e k is the vector of differences between the desired output yd and the model output y at time k in the future. The matrix Q is used to scale the variables and to define the relative importance of the errors for the different output values. The matrix R has the same function, i.e., to weigh the incremental variations of the control actions. It is important to note that the control variable penalty term in the criterion Ji is expressed with respect to the increment of u and not to its variation around its nominal value. This formulation warrants the cancellation of the steady-state deviation of the process output to the setpoint. At each time i, the criterion Ji is minimized with respect to u i , . . . u iqHy1 , then at the next time i q 1, again the criterion Jiq1 is minimized. This is why this strategy is also called a moving horizon control. The control scheme of Fig. 10 can be improved by adding reference models which allow an easy tuning of the closed loop dynamics. Reference models can be added to specify the setpoint tracking dynamics and the disturbance rejection dynamics. The reference models Md and Mr can
also be viewed as filters which modify the dynamics of setpoint and disturbance changes. Although the control actions weighting factors can be used to constrain the variations of these control actions, strict inequality constraints can be added to the criterion. They can be formulated as: u min F u F u max < D u < F D max where D u is the incremental variation of u. Also, inequality constraints can be formulated on the process model outputs y or even on the process model states. The resulting algorithm is a controller which uses the state vector x of the process model and the reference trajectory yd . For the case under study, we have: 0.1 0 1 0 Rs Qs 0 0.1 0 1 1 0 1 q 25s Md s 1 0 1 q 5s 1 0 1 q 25s Mr s 1 0 1 q 5s 9. Adaptive predictive control algorithm Industrial processes are generally characterized by nonlinear time-varying dynamics w9x. Most of nonlinear systems, however, may be modelled by locally linear approximations. The input–output signals are assumed to be large enough to avoid small signal nonlinearities Žvalve stiction, etc.., but not too large to reach large signal nonlinearities Žhysteresis, saturation, etc...
Fig. 10. Full multivariable internal model controller.
A. Pomerleau et al.r Powder Technology 108 (2000) 103–115
111
Fig. 11. Distributed PID controller results.
The plant input–output behaviour which is to be adaptively modelled is described by the following disturbed linear discrete model: A Ž qy1 . y Ž t . s qyk B Ž qy1 . u Ž t . q Õ Ž t . ÕŽ t. s
C Ž qy1 . D Ž qy1 .
j Ž t.
where A, B and C are polynomials in the backward shift operator of order na, nb and nc. The parameter adaptation mechanism is the heart of the adaptive control system: the performance of the control algorithm depends directly on the quality of this mechanism. The presented method deals with parameter estimation based on the information collected on the process, i.e., the input–output measurements.
112
A. Pomerleau et al.r Powder Technology 108 (2000) 103–115
Fig. 12. Internal model with explicit decoupler results.
For identification purposes, the process output is expressed in terms of a parameter vector u and a regressor or data vector: y Ž t . s u Tf Ž t . q Õ Ž t .
u s w a1 ,a2 , . . . ,a n a ,b 0 ,b 1 , . . . ,bn b x
T
f Ž t . s yy Ž t y 1 . , . . . ,y y Ž t y na . , u Ž t y k . , . . . , u Ž t y k y nb .
T
A. Pomerleau et al.r Powder Technology 108 (2000) 103–115
It is always advisable to provide the parameter estimator with well-conditioned signals. The signal conditioning is provided by scaling, filtering and data normalization. Low-pass filtering of the measurements may be used in order to reduce the high frequency modes due to both
113
noise and unmodelled dynamics. This ensures smooth behaviour of parameter estimation and reduces their variations. The long-range predictive control objective consists of minimizing the sum of squares of the errors between
Fig. 13. Multivariable predictive controller results.
A. Pomerleau et al.r Powder Technology 108 (2000) 103–115
114
Fig. 14. Adaptive decentralized predictive controller results.
predicted and desired output trajectories with an additional term weighting projected control increment: JsE
½
Ny 2 Ý Ž yŽ tqj. ywŽ tqj. .
jsNm Nu
ql
Ý jsNm
2
D u Ž t q j y 1 . rt
5
subject to D uŽ t q j y 1. s 0 when j ) Nu, where E is the mathematical conditional expectation given information up to time t, w Ž t . is the desired output, Ny the maximum prediction horizon, Nm the minimum prediction horizon, Nu is the control horizon and l is a weighting upon future control increments. This criterion is roughly the same as the previous one except that it is monovariable and weighting factors are limited only to one l.
A. Pomerleau et al.r Powder Technology 108 (2000) 103–115
115
10. Results and discussion A benchmark has been used to test the various control algorithms. It consists of the following: Time Žmin. 50 225 300
Event CL setpoint SD setpoint CL setpoint
575 825 1075 1300 1500
Ore feed Number of hydrocyclones Diameter of apex CL setpoint SD setpoint
Based on the observations of the results ŽFigs. 11–14., the following conclusions can be drawn: – For setpoint changes, the algorithms that have a decoupler eliminate totally or partially the disturbance on the variable that was not imposed a setpoint change as was expected. – The algebraic internal model and the multivariable predictive controller perform equally good. The advantage of the former is to be easily implementable on DCS Ždistributed control system. that has bloc programmation available. It may require though better know-how from the designer. – The adaptive controller performs better than the fixed controllers in face of parameter variations occurring closed to the inputs of the process Žore feed change., but no better for parameter variations near the output Žnumber of hydrocyclones..
11. Conclusion From the observations of the results on the grinding process, which can be extrapolated to any process from a control point of view, one can conclude the following. – Fixed structure controllers Že.g., PID. perform as well as model-based controllers for processes where the delay is relatively small compared to the dominant time constant Ž u - 5T . since the models are rarely of order higher than two. – Distributed controllers perform as well as multivariable controllers in regulation if the coupling of the process is taken into account in the design and the proper pairing is used for the dynamics required. – Algebraic multivariable controllers perform as well as optimal multivariable controllers and they have exactly the same limitations Že.g., perfect decoupling might be impossible..
from 191.7 to 211.7 trh from 76.9% to 75.65% from 211.7 to 231.7 trh from 75.65% to 74.4% from 100% hard coarse to 50% hard coarser50% soft fine from 4 to 3 from 4.64 to 4.44 cm from 231.7 to 171.7 trh from 74.4% to 75.65% – Algebraic tuning methods may require more knowhow than optimal controllers, but are easier to implement on industrial DCS. – Predictive controller and algebraic controllers perform equally in a stochastic environment if there is no noise model available. – Adaptive controllers perform better than fixed controllers for parametric disturbances or soft nonlinearities. It must be noted though that identification is very difficult in regulation where external disturbances act on the process. They also have the advantage of facilitating the tuning principally for distributed controllers.
References w1x D. Hodouin, Y. Marcotte, A. Pomerleau, F. Flament, Predictive Control of Grinding Circuit: An Evaluation by Dynamic Simulation. w2x A. Desbiens, A. Pomerleau, K. Najim, Int. J. Miner. Process. 41 Ž1994. 17–31. w3x A. Desbiens, K. Najim, A. Pomerleau, D. Hodouin, Distributed Partial State Reference Model Adaptive Control — Practical Aspects and Application to a Grinding Circuit, Optimal Control Applications and Methods 18 Ž1997. 29–47. w4x F. Flament, R. del Villar, R. Lanthier, Computer aided design of a control strategy for an industrial grinding circuit, in: R. Poulin, R.C.T. Pakalnis, A.L. Mular ŽEds.., Proc. of the 2nd Can. Conf. on Computer Applications in the Mineral Industry, Vol. 1, 1991, pp. 337–348. w5x D. Hodouin, Y. Dube, ´ R. Lanthier, Stochastic simulation of filtering and control strategies for grinding circuit, Int. J. Miner. Process. 22 Ž1988. 261–274. w6x E.H. Bristol, On a new measure of interaction for multivariable process control, IEEE Trans. Autom. Control AC 11 Ž1966. 133. w7x H.P. Huang, M. Ohshima, I. Hashimoto, J. Proc. Cont. 4 Ž1994. 15–17. w8x A. Desbiens, A. Pomerleau, D. Hodouin, Frequency based tuning of SISO Controllers for two by two Processes, IEE Proc. on Control Theory and Applications 43 Ž1996. 49–56. w9x K. Najim, Process Modelling and Control in Chemical Engineering, Marcel Dekker, New York, 1989, 520 pp.
A survey of grinding circuit control methods: from decentralized PID controllers to multivariable predictive controllers Andre´ Pomerleau a b
a,)
´ Gagnon , Daniel Hodouin b, Andre´ Desbiens a , Eric
b
Department of Electrical Engineering, LaÕal UniÕersity, Quebec, Canada G1K 7P4 Department of Mining and Metallurgy, LaÕal UniÕersity, Quebec, Canada G1K 7P4 Accepted 20 September 1999
Abstract A conventional grinding circuit consisting of one open-loop rod mill and one closed-loop ball mill is essentially a two-input= twooutput system, assuming that the classifier pump box level is controlled by a local loop. The inputs are the ore and water feed rates and the outputs are the product fineness and the circulating load. The design problem is to find a control algorithm and a tuning procedure which satisfy specified servo and regulatory robust performances. A first approach is to use decentralized PID controllers and systematic tuning methods which take into account loop interactions. Another technique consists of adding decouplers or pseudo-decouplers to the decentralized controllers. Finally, the design of a fully multivariable controller is a possible option. To face the problem of performance robustness related to change of process dynamics, two options are studied. A design criterion involving the minimization of a penalized quadratic function on a future trajectory can be used. A second alternative is to track process dynamics changes using adaptive process modelling. The paper will present a comparison of these various strategies, for a simulated grinding circuit. A benchmark test, involving a sequence of disturbances Žgrindability, feed size distribution, change of cyclone number . . . . and setpoint changes, is used to compare the performances of the controllers. q 2000 Elsevier Science S.A. All rights reserved. Keywords: Grinding circuit; Distributed PID; Internal model control; Multivariable predictive control; Distributed adaptive predictive control
1. Introduction The objective of this paper is to study the performances of various control algorithms acting on a grinding circuit in face of different disturbances and setpoint changes. Even though numerous articles w1–3x have studied this control problem, it is hard to compare the results since the benchmarks are different. In this article, the characteristics of the process, from a control point of view are first founded. From there, various control algorithms will be applied and general conclusions on the performances will be drawn. The algorithms studied are: –
Decentralized PID controllers which are characterized by this fixed structure. Algebraic internal model with explicit decouplers. They are characterized by a model-based structure. Full multivariable predictive controllers. Distributed adaptive predictive controllers. These are
– – – )
Corresponding author.
characterized by variable parameters as function of the process model variation.
2. Process description The grinding process is a size reduction operation used in the mineral industry to liberate the valuable minerals from the discardable gangue. Fig. 1 shows the flowsheet of the plant: it is a standard two-stage grinding circuit which consists of a rod mill in open circuit and a ball mill in closed circuit with a hydrocyclone classifier. Two local controllers are assumed to be part of the process. Sump level is maintained at a constant value using a variable speed pump controlled by a PID. The percentage of solids in the feed is kept around a desired set-point using a feed water addition rate proportional to the fresh ore feed rate. Grinding typically represents around 50% of the total expenditure of the concentrator operation. This step prepares the ore for the concentration step such as flotation or magnetic separation. Therefore, it determines the feed characteristics to these subsequent separation stages and
0032-5910r00r$ - see front matter q 2000 Elsevier Science S.A. All rights reserved. PII: S 0 0 3 2 - 5 9 1 0 Ž 9 9 . 0 0 2 0 7 - 7
104
A. Pomerleau et al.r Powder Technology 108 (2000) 103–115
Fig. 1. Grinding circuit plant flowsheet.
influences the efficiency of the valuable minerals recovery process. The true objective of grinding is to obtain a proper liberation of the minerals coupled with a proper particle size distribution in order to maximize recovery and concentrate value without overgrinding, while keeping specific comminution energy and reagent consumption as low as possible, taking into account the prevailing economic conditions. These objectives being highly complex, simpler ones must be formulated. Generally, as the liberation degree cannot be measured on line, the production objective is to maintain at a set-point value the percentage of the product particles smaller than the reference size Žfor instance 70% smaller than 74 mm..
In parallel, the circulating load Žtons per hour of solids entering the ball mill. is to be maximized. High recirculation at a fixed product size means lower energy consumption. Excessive circulating load can result in classifier, mill or pump overload conditions. The selected controlled variables are therefore the size distribution and the circulating load. The fresh ore feed rate and the water addition to the sump box are the available manipulated variables. Since industrial design and testing are costly, the dynamic grinding circuit simulator DYNAFRAG is used w4x. It simulates the mineral processing plant with appropriate dynamic phenomenological models Žset of mass balance equations for the pulp, the water and the particles of different hardness.. The detailed mathematical formulation of the models involved in the simulator is given elsewhere w5x. From 100 to 150 simultaneous differential equations are to be solved. With DYNAFRAG, it is also possible to simulate disturbances such as sensor and actuator noises, natural variations in ore characteristics Žparticularly hardness changes., changes in circuit performance and mechanical upsets Žblockages in pumps and classifiers, etc... The simulator operates in a multi-computer environment. The ‘‘process’’ and the ‘‘control room’’ are physically separated. They are both simulated on their own micro-computer and are interconnected through digital-to-analog and analog-to-digital interfaces. This multi-computer architecture resembles industrial reality.
Fig. 2. Open-loop step response.
A. Pomerleau et al.r Powder Technology 108 (2000) 103–115
3. Process model Fig. 2 illustrates the open-loop responses of DYNAFRAG Žsampling period: 0.5 min.. Fig. 2a and b shows the size distribution ŽSD: percentage passing 74 mm. and the circulating load ŽCL. following an 8.8 m3rh step on the pump box water ŽPBW. at time t s 10 min. The responses to a 6 trh step on the rod mill feed ŽRMF. at time t s 10 min are given in Fig. 2c and d. Fig. 3 shows the same responses with white noise added to the measured variables. The same level of noise is kept during the closed-loop simulations. The linear process model identified at the operating point is: 0.6925ey5 s
CL s RMF CL s PBW
1 q 18 s 0.475 Ž 1 y 3.5s .
Ž 1 q 4.5s . Ž 1 q s .
SD s
y1.63 Ž 1 y 2.5s . ey8 s
RMF
Ž 1 q 28 s . Ž 1 q s . SD 0.393 Ž 1 q 92 s . eys s PBW Ž 1 q 21 s . Ž 1 q 2 s . The time constants are expressed in minutes. From these transfer functions, one concludes that a grinding
105
circuit is a highly coupled process where PBW addition has a very fast short-term effect on the size distribution, but where the RMF has the dominating long-term effect.
4. Selected control algorithms There are two main approaches to select the control loop architecture for multivariable processes: – decentralized Single-Input Single-Output ŽSISO. loops – multivariable controllers, where the latter can also be subdivided in two depending if the decoupler is explicitly expressed or not. Fig. 4 shows a multivariable process controlled by distributed SISO controllers. In this structure, the proper inputroutput pairing has to be chosen and the tuning of one controller depends on the tuning of the second controller. In a structure where the decoupler is explicitly expressed as shown in Fig. 5, a decoupler structure has to be chosen and the facility of decoupling depends on the choice of the coupling used. The controller tuning of the SISO controllers is done on the combined system including the process and the decoupler. Finally, in a fully multivariable controller as shown in Fig. 6, the designer does not need to choose the pairing, but the specifications have to be feasible.
Fig. 3. Open-loop step response Žnoise environment..
106
A. Pomerleau et al.r Powder Technology 108 (2000) 103–115
Fig. 4. Distributed SISO controllers.
The first studied algorithm is the standard decentralized PID controller which is usually the reference to compare different algorithms. PID controllers present some drawbacks. First, their performances are limited by the properties of the fixed structure. Second, the inputroutput pairing problem has to be solved for multivariable systems. Finally, an efficient tuning method has to be used to obtain the desired dynamics of the controlled loops. The second studied algorithm consists of distributed internal model controllers applied to the process preceded by a partial decoupler. The decoupler partially solves the pairing problem and facilitates the tuning. The internal model control structure which is model-based eliminates the problems of a fixed structure. The third studied algorithm is a full multivariable optimal predictive controller. It totally solves the pairing problem and takes into account the compromise between the performances and the manipulated variables energy through a quadratic criterion. Finally, a distributed adaptive predictive controller is used. Compared to the other algorithms, the adaptive part permits to take into account the nonlinearities of the
process and the time evolution of the parameters. It also solves the tuning problems of the distributed control algorithm. The time specification on the circulating load has been arbitrarily fixed to 25 min. This corresponds approximatively to a dynamic similar to the open-loop function between the CL and RMF. The specification on the size distribution has been fixed to 5 min. Here, since the open-loop transfer functions relating the size distribution to the RMF or to the PBW have zeroes, it is more difficult to relate the specification to the open-loop time constants.
5. Pairing of controlled and manipulated variables Distributed control consists of using n SISO to control an n = n process as shown in Fig. 4. The main difficulties of distributed control lie in the choice of the proper input–output pairing and in the fact that the tuning of each SISO controller depends on the tuning parameters of the second SISO controller. When manipulated variables are not properly selected, interactions between controlled and
Fig. 5. Distributed SISO controllers with decoupler.
A. Pomerleau et al.r Powder Technology 108 (2000) 103–115
107
Fig. 6. Full multivariable controller.
manipulated variables can result in undesirable control loop interactions, leading to poor control performances. The well-known Bristol’s relative gain array ŽRGA. w6x being based only on steady-state information does not consider the process dynamics and is inappropriate for many processes. This is especially the case when dynamics in transversal transfer functions are fast compared to direct transfer functions and present smaller static gain than
direct transfer functions. The presence of zeros in the transfer functions also degrades the performance of RGA. Consequently, a dynamic method must be employed. The basic and simple idea of the RDG ŽRelative Dynamic Gain. method is to extend the calculation of the Bristol’s relative gain for all frequencies. The method used here is the generalized dynamic relative gain w7x. For the process under study, the RDG frequency responses are
Fig. 7. RDG frequency response.
108
A. Pomerleau et al.r Powder Technology 108 (2000) 103–115
given in Fig. 7. For the dynamics required by the specification, the output Y1 , is paired with U1.
The system can be drawn as depicted in Fig. 8. For each closed loop, the control objectives are: Ž1. no permanent error and Ž2. a second-order set point tracking response with specified dynamics. Thus, the corresponding closedloop transfer functions can be selected as follows:
6. Algorithm 1: Decentralized extended PIDs Y1 Ž s . To tune the controllers, a frequency-based method, in which the open-loop characteristics are deduced from the closed-loop dynamics is used w8x. First, for tuning a controller, it is imperative to know the dynamics of the process to be controlled. With monovariable processes, the method consists of opening the loop and evaluating the transfer function seen by the controller. For multivariable processes, the loop with the controller under consideration is opened while all the other loops are kept closed. For a 2 = 2 process, the transfer function seen by the controller of the first loop, GC1Ž s ., is: G 1 Ž s . s G 11 Ž s . y
G 12 Ž s . G 21 Ž s . GC 2 Ž s . 1 q GC 2 Ž s . G 22 Ž s .
.
The transfer function seen by the second controller GC2 Ž s . is: G 2 Ž s . s G 22 Ž s . y
G 12 Ž s . G 21 Ž s . GC1 Ž s . 1 q GC1 Ž s . G 11 Ž s .
.
As can be seen, there exists an interdependence between GC1Ž s . and GC2 Ž s .. Indeed, their individual tuning depends on the choice of the other controller. Consequently, to evaluate the unknowns GC1Ž s . and GC2 Ž s ., further information is necessary. This supplementary information is given by the closed-loop specifications which can be translated into open-loop char-acteristics.
Ž 1 q t 11 s . Ž 1 q t 21 s .
Z1 Ž s . Y2 Ž s . Z2 Ž s .
1 y t 01 s
s
1 y t 02 s
s
Ž 1 q t 12 s . Ž 1 q t 22 s .
The open-loop characteristics deduced from the tracking closed-loop specifications are: G 0 L 1 Ž s . s GC1 Ž s . G 1 Ž s . 1 y t 01 s
s
s Ž t 11 q t 21 q t 01 q t 11t 21 s . G 0 L 2 Ž s . s GC 2 Ž s . G 2 Ž s . 1 y t 02 s
s
s Ž t 12 q t 22 q t 02 q t 12t 22 s . The transfer functions GC1Ž s . and GC2 Ž s . of the controllers are obtained by solving the previous equations. Due to their quadratic form, it leads to two Bode plots for both controllers. Knowing the required sign of the controller gain and the necessary presence of an integrator, the appropriate solution can be selected without any difficulty. The last step in tuning the controllers consists of defining the transfer functions GCP1Ž s . and GCP2 Ž s . that best approximate the frequency response GC1Ž j v . and GC2 Ž j v ., while ensuring a stable closed-loop system. The structure of each controller is imposed and an evaluation of the
Fig. 8. Distributed SISO controllers equivalent system.
A. Pomerleau et al.r Powder Technology 108 (2000) 103–115
parameters is performed using a constrained objective function. The following structure is used: GCP i Ž s . s
K C Ž 1 q Ti s . Ž 1 q Td1 s q Td 2 s 2 . Ti s Ž 1 q Tf 1 s q Tf 2 s 2 .
.
2
GC1 Ž s . s GC 2 Ž s . s
1.28 Ž 1 q 85.4 s . Ž 1 q 19.15s q 73.4 s . 85.4 s Ž 1 q 17.81 s q 90.6 s . 0.32 Ž 1 q 2.4 s . Ž 1 q 4.77s q 9.06 s 2 . 2.4 s Ž 1 q 6.4 s q 7.8 s 2 .
7. Control algorithm 2: algebraic multivariable control Multivariable algebraic control can be done either with a full multivariable controller or by using a decoupler in cascade with monovariable controllers. The latter one, as shown in Fig. 9, is used here. The methodology consists of: – evaluating the decoupler, – tuning the monovariable controllers on the model including the multivariable process model and the decoupler. For many processes Že.g., grinding circuit., it is impossible to realize a perfect decoupler. This happens with non-minimal phase processes or with systems with longer delays in the direct transfer functions. For the system under study, the decouplers are: D 12 Ž s . s y s
G 22 Ž s .
1.633 Ž 1 y 2.5s . Ž 1 q 2 s . Ž 1 q 21 s .
D 21 Ž s . s y
Ž 1 q 92 s .
y7 s
e
G 21 Ž s . G 11 Ž s .
0.475 s
0.6925 ey5 s
Ž 1 q 18 s . Ž 1 y 2.5s . Ž 1 q 2 s . Ž 1 q 21s . ey7 s q 1.974 Ž 1 q s . Ž 1 q 28 s . Ž 1 q 92 s . Ž 1 y 3.5s . = Ž 1 q 4.5s . Ž 1 q s . 0.393 Ž 1 q 92 s . eys G2 Ž s . s Ž 1 q 2 s . Ž 1 q 21s . Ž 1 y 3.5s . Ž 1 q 18 s . q 1.120 Ž 1 q 4.5s . Ž 1 q s . Ž 1 y 2.5s . ey8 s = Ž 1 q s . Ž 1 q 28 s . These transfer functions can be approximated by: G1Ž s . s G2 Ž s . s
2.666ey5 s
Ž 1 q 85s . 1.513 Ž 1 q 5s . eys Ž 1 q 3s .
2
The distributed internal model controllers are then tuned into these transfer functions. 8. Control algorithm 3: multivariable state-space internal model predictive control
G 12 Ž s .
0.393 Ž 1 q 28 s .
Under the assumption that the system is satisfactorily decoupled for neglecting the remaining coupling for tuning consideration, the transfer functions for the system under study are given by G1Ž s . s
The resulting controllers are:
109
Ž 1 y 3.5s . Ž 1 q 18 s . 0.69 Ž 1 q 4.5s . Ž 1 q s .
D 12 Ž s . realizes a perfect decoupler, while D 21Ž s . is a pseudo-decoupler.
Here the controller is designed to control the process model and not the process itself. The difference between the model and the process outputs, respectively, y and y P submitted to the same control action u, is then interpreted as an overall estimation d of the disturbance d. This estimated value embeds in fact both disturbances and modelling errors. The value d is used to modify the setpoint Žreference value r .. This biased setpoint is the desired value yd that must be tracked by the process model. The controller calculates u in such a way to cancel the difference e between yd and the model output ym .
Fig. 9. Internal model controller with explicit decoupler.
A. Pomerleau et al.r Powder Technology 108 (2000) 103–115
110
Fig. 10 summarizes the basic principles of this control architecture. It is generally named Internal Model Control ŽIMC.. At time i, the sequence of control actions to be applied in the future Ž u i , u iq1 , . . . u iqHy1 ., where u k is the vector of control actions at time k and H the control horizon, is selected such that the following quadratic performance criterion is minimized: iqH
Ji Ž u i . . . u iqHy1 . s
e kT Q e k
Ý ksiq1
iqH
q
Ý Ž u ky 1 y u ky2 . T ksiq1
=R Ž u ky 1 y u ky2 . In the equation, e k is the vector of differences between the desired output yd and the model output y at time k in the future. The matrix Q is used to scale the variables and to define the relative importance of the errors for the different output values. The matrix R has the same function, i.e., to weigh the incremental variations of the control actions. It is important to note that the control variable penalty term in the criterion Ji is expressed with respect to the increment of u and not to its variation around its nominal value. This formulation warrants the cancellation of the steady-state deviation of the process output to the setpoint. At each time i, the criterion Ji is minimized with respect to u i , . . . u iqHy1 , then at the next time i q 1, again the criterion Jiq1 is minimized. This is why this strategy is also called a moving horizon control. The control scheme of Fig. 10 can be improved by adding reference models which allow an easy tuning of the closed loop dynamics. Reference models can be added to specify the setpoint tracking dynamics and the disturbance rejection dynamics. The reference models Md and Mr can
also be viewed as filters which modify the dynamics of setpoint and disturbance changes. Although the control actions weighting factors can be used to constrain the variations of these control actions, strict inequality constraints can be added to the criterion. They can be formulated as: u min F u F u max < D u < F D max where D u is the incremental variation of u. Also, inequality constraints can be formulated on the process model outputs y or even on the process model states. The resulting algorithm is a controller which uses the state vector x of the process model and the reference trajectory yd . For the case under study, we have: 0.1 0 1 0 Rs Qs 0 0.1 0 1 1 0 1 q 25s Md s 1 0 1 q 5s 1 0 1 q 25s Mr s 1 0 1 q 5s 9. Adaptive predictive control algorithm Industrial processes are generally characterized by nonlinear time-varying dynamics w9x. Most of nonlinear systems, however, may be modelled by locally linear approximations. The input–output signals are assumed to be large enough to avoid small signal nonlinearities Žvalve stiction, etc.., but not too large to reach large signal nonlinearities Žhysteresis, saturation, etc...
Fig. 10. Full multivariable internal model controller.
A. Pomerleau et al.r Powder Technology 108 (2000) 103–115
111
Fig. 11. Distributed PID controller results.
The plant input–output behaviour which is to be adaptively modelled is described by the following disturbed linear discrete model: A Ž qy1 . y Ž t . s qyk B Ž qy1 . u Ž t . q Õ Ž t . ÕŽ t. s
C Ž qy1 . D Ž qy1 .
j Ž t.
where A, B and C are polynomials in the backward shift operator of order na, nb and nc. The parameter adaptation mechanism is the heart of the adaptive control system: the performance of the control algorithm depends directly on the quality of this mechanism. The presented method deals with parameter estimation based on the information collected on the process, i.e., the input–output measurements.
112
A. Pomerleau et al.r Powder Technology 108 (2000) 103–115
Fig. 12. Internal model with explicit decoupler results.
For identification purposes, the process output is expressed in terms of a parameter vector u and a regressor or data vector: y Ž t . s u Tf Ž t . q Õ Ž t .
u s w a1 ,a2 , . . . ,a n a ,b 0 ,b 1 , . . . ,bn b x
T
f Ž t . s yy Ž t y 1 . , . . . ,y y Ž t y na . , u Ž t y k . , . . . , u Ž t y k y nb .
T
A. Pomerleau et al.r Powder Technology 108 (2000) 103–115
It is always advisable to provide the parameter estimator with well-conditioned signals. The signal conditioning is provided by scaling, filtering and data normalization. Low-pass filtering of the measurements may be used in order to reduce the high frequency modes due to both
113
noise and unmodelled dynamics. This ensures smooth behaviour of parameter estimation and reduces their variations. The long-range predictive control objective consists of minimizing the sum of squares of the errors between
Fig. 13. Multivariable predictive controller results.
A. Pomerleau et al.r Powder Technology 108 (2000) 103–115
114
Fig. 14. Adaptive decentralized predictive controller results.
predicted and desired output trajectories with an additional term weighting projected control increment: JsE
½
Ny 2 Ý Ž yŽ tqj. ywŽ tqj. .
jsNm Nu
ql
Ý jsNm
2
D u Ž t q j y 1 . rt
5
subject to D uŽ t q j y 1. s 0 when j ) Nu, where E is the mathematical conditional expectation given information up to time t, w Ž t . is the desired output, Ny the maximum prediction horizon, Nm the minimum prediction horizon, Nu is the control horizon and l is a weighting upon future control increments. This criterion is roughly the same as the previous one except that it is monovariable and weighting factors are limited only to one l.
A. Pomerleau et al.r Powder Technology 108 (2000) 103–115
115
10. Results and discussion A benchmark has been used to test the various control algorithms. It consists of the following: Time Žmin. 50 225 300
Event CL setpoint SD setpoint CL setpoint
575 825 1075 1300 1500
Ore feed Number of hydrocyclones Diameter of apex CL setpoint SD setpoint
Based on the observations of the results ŽFigs. 11–14., the following conclusions can be drawn: – For setpoint changes, the algorithms that have a decoupler eliminate totally or partially the disturbance on the variable that was not imposed a setpoint change as was expected. – The algebraic internal model and the multivariable predictive controller perform equally good. The advantage of the former is to be easily implementable on DCS Ždistributed control system. that has bloc programmation available. It may require though better know-how from the designer. – The adaptive controller performs better than the fixed controllers in face of parameter variations occurring closed to the inputs of the process Žore feed change., but no better for parameter variations near the output Žnumber of hydrocyclones..
11. Conclusion From the observations of the results on the grinding process, which can be extrapolated to any process from a control point of view, one can conclude the following. – Fixed structure controllers Že.g., PID. perform as well as model-based controllers for processes where the delay is relatively small compared to the dominant time constant Ž u - 5T . since the models are rarely of order higher than two. – Distributed controllers perform as well as multivariable controllers in regulation if the coupling of the process is taken into account in the design and the proper pairing is used for the dynamics required. – Algebraic multivariable controllers perform as well as optimal multivariable controllers and they have exactly the same limitations Že.g., perfect decoupling might be impossible..
from 191.7 to 211.7 trh from 76.9% to 75.65% from 211.7 to 231.7 trh from 75.65% to 74.4% from 100% hard coarse to 50% hard coarser50% soft fine from 4 to 3 from 4.64 to 4.44 cm from 231.7 to 171.7 trh from 74.4% to 75.65% – Algebraic tuning methods may require more knowhow than optimal controllers, but are easier to implement on industrial DCS. – Predictive controller and algebraic controllers perform equally in a stochastic environment if there is no noise model available. – Adaptive controllers perform better than fixed controllers for parametric disturbances or soft nonlinearities. It must be noted though that identification is very difficult in regulation where external disturbances act on the process. They also have the advantage of facilitating the tuning principally for distributed controllers.
References w1x D. Hodouin, Y. Marcotte, A. Pomerleau, F. Flament, Predictive Control of Grinding Circuit: An Evaluation by Dynamic Simulation. w2x A. Desbiens, A. Pomerleau, K. Najim, Int. J. Miner. Process. 41 Ž1994. 17–31. w3x A. Desbiens, K. Najim, A. Pomerleau, D. Hodouin, Distributed Partial State Reference Model Adaptive Control — Practical Aspects and Application to a Grinding Circuit, Optimal Control Applications and Methods 18 Ž1997. 29–47. w4x F. Flament, R. del Villar, R. Lanthier, Computer aided design of a control strategy for an industrial grinding circuit, in: R. Poulin, R.C.T. Pakalnis, A.L. Mular ŽEds.., Proc. of the 2nd Can. Conf. on Computer Applications in the Mineral Industry, Vol. 1, 1991, pp. 337–348. w5x D. Hodouin, Y. Dube, ´ R. Lanthier, Stochastic simulation of filtering and control strategies for grinding circuit, Int. J. Miner. Process. 22 Ž1988. 261–274. w6x E.H. Bristol, On a new measure of interaction for multivariable process control, IEEE Trans. Autom. Control AC 11 Ž1966. 133. w7x H.P. Huang, M. Ohshima, I. Hashimoto, J. Proc. Cont. 4 Ž1994. 15–17. w8x A. Desbiens, A. Pomerleau, D. Hodouin, Frequency based tuning of SISO Controllers for two by two Processes, IEE Proc. on Control Theory and Applications 43 Ž1996. 49–56. w9x K. Najim, Process Modelling and Control in Chemical Engineering, Marcel Dekker, New York, 1989, 520 pp.