Optimization of control parameters of a hot cold controller by means of Simplex type methods

This paper describes a hot/cold controller for regulating crystallization operations. The system was identified with a common method (the Broida method) and the parameters were obtained by the Ziegler-Nichols method. The paper shows that this empirical method will only allow a qualitative approach to regulation and that, in some instances, the parameters obtained are unreliable and therefore cannot be used to cancel variations between the set point and the actual values. Optimization methods were used to determine the regulation parameters and solve this identcation problem. It was found that the weighted centroid method was the best one.

This paper describes a hot/cold controller for regulating crystallization operations. The system was identified with a common method (the Broida method) and the parameters were obtained by the Ziegler-Nichols method. The paper shows that this empirical method will only allow a qualitative approach to regulation and that, in some instances, the parameters obtained are unreliable and therefore cannot be used to cancel variations between the set point and the actual values. Optimization methods were used to determine the regulation parameters and solve this identcation problem. It was found that the weighted centroid method was the best one.
Until the 1970s, regulation, which is important in bulk chemistry, was based on traditional controllers: analogue, pneumatic or electrical. When microprocessors began to appear, numerical PID controllers were available and proved much easier to use and preset. In addition, controllers quickly became less expensive. To maintain a temperature in chemistry, often involves complex regulation, either in cascade or of the hot/cold type. In both cases, the transfer functions of the systems have to be known so that the parameters involved in the regulation can be identified. This requires the use of empirical methods, such as the Broida method, which is associated with the Ziegler-Nichols method. These methods are easy to implement in terms of obtaining the parameters of the controller, but they require numerous 'tries' in order to reach the objective. This paper reports on the use of optimization methods to determine the optimal parameters of a hot/cold controller, which was being used for the regulation of the temperature in a reactor-crystallizer. This paper describes the experimental set-up for a laboratory-scale crystallization apparatus and shows how the parameters for the PID controller were determined using the Broida and Ziegler-Nichols methods. The optimization of the settings are described, and the results obtained by different optimization methods are compared: simplex, modified simplex (Nelder and Mead), multi-move and weighted centroid method (WCM).
The hot/cold controller The aim was to regulate temperature control during a crystallization operation [1]. Temperature is controlled over time, so nucleation and growth rate can also be regulated. The control over temperature must be as accurate as possible and requires a good regulation around the desired value. This value varies over time and follows a polynomial law of the type T a / bt3, where T is the temperature, the duration of the operation, and a and b are the parameters. The Pyromat 300 and Programat controllers, manufactured by Schlumberger, were used. The program emitter allows the set value to be changed over time.
The Pyromat 300 series are programmable PID numerical controllers. There is a channel dedicated to the cooling system and another dedicated to the heating system. The two channels are completely independent and can be set to analogue or numerical PID modes, They offer three possibilities for the heating and cooling channels: modulated/modulated, modulated/analogue, and analogue/modulated. The 'cold' and 'hot' channels can overlap for up to -10 to 20% of the command signal. Therefore, both the heating and cooling rate can be controlled.

The heating loop
The heating loop consists of a Ptl00 probe (TP) immersed in a 10 thermostatted bath filled with water, with a circulation pump (flow rate 1000 l/h), a heating resistance (HR) of 1000 W connected to a power unit (PU) through a controller (C) (see figure 1). The power unit is controlled by the analogue output (4-20 mA) of the 'hot' channel of the controller. This signal is converted by the power unit's interface to a percentage (0-100%) of the total heating power. The cooling loop consists of the same controller (C) and Ptl00 probe (TP) immersed in the thermostatted bath. Tap water flows through a coil (average temperature 16 C) when an electrovalve (EV) is opened. For crystallization processes requiring operations below 20C, a water/glycol mix cooled to -10C is used instead of tap water. A static relay is activated to allow the passage of this mix. This regulation loop is shown in figure 2. The cooling system is controlled through the 'cold' channel with a discontinuous modulated signal (0-10 V). Before using controllers without autosetting capabilities, the regulation parameters must be determined as follows: The proportional band, expressed as the variation in a percentage of the total measurement range.
The integral time (in seconds) which represents the necessary time lapse before the signal output variation equals the variation caused by the proportional action. The derived time which allows the effects of the proportional action to be amplified.
The modulation time.
These parameters are set during the process identification phase. the evolution of the controlled value over time allows the transfer function to be calculated and also, therefore, the most favourable conditions for the stability of the system, i.e. the parameters to be set on the controller. The identification of a process is the key to finding its transfer function. Identification is often a delicate task, but it is necessary to calculate the controller's actions. It is often accomplished by creating perturbations. The investigation of the system's response was conducted in this study in an open loop, using the Broida method [2]. Implementation of the Broida method For system identification purposes, the controller must be set to manual mode, i.e. disconnected. The system must also be in the permanent regime, with the output and command signals, O(s)and C(s), constant. This is an open loop process. A perturbation is then applied to the system; the most common perturbation used is an incremental step. The command signal is incremented, and the output signal is recorded.
The transfer function for a first order system with pure delay can be expressed as follows: where Gs-static gain; "r-time constant; 0-pure delay.
The theoretical Broida graph (see figure 3) (first order system with pure delay) shows the output value to be equal to 0"28 times the final variation at the time (t 0); and, at (t2-0), 0"40 of the final variation: These equations lead to the value of the time constant and 0" 7"-5"5(t --tl) (4) 0-2"8. tl 1.8.t2 (5) Parameter setling The majority of regulation loops can be accomplished simply with an empirical adjustment of the parameters of universal controllers. Other loops are more complicated and it is necessary to know the dynamics of the process concerned. Most systems are naturally stable and, when displaced out of their equilibrium state by a perturbation or a corrective action, they tend to attain a new state of equilibrium after a certain time. This dynamic phase of the system is a characteristic function of the process.
The optimal parameter values for a PID controller depend only on the characteristics of the system to be regulated and on the controller which is being used.   characteristics of a process are determined with an empirical method. The influent variable is changed in steps covering the full range of possible values, and a correlation is derived from observing the resulting effects on the controlled value.

Indicial response
Indicial response of the system to incremental steps was studied in order to determine the types of action that a controller needs to perform. This analysis is conducted in an open loop (the controller does not react to the temperature read from the probe). The programmer is disconnected from the controller (see figure 4).
An incremental step AY, corresponding to a variation of the heating or cooling power of Y% at the controller's output, is applied. This is equivalent to sending a constant signal to the regulating device. The resulting increase (AX) of the response signal (controlled value) is measured until the system stabilizes. Static gain is:  gain can be derived from tl and t2 (see table 3). Table 4 shows the formulae established by Ziegler and Nichols to determine the values of the PID parameters of a controller. The proportional band Xp is expressed as a percentage of the setting range. The parameters Xp and Ti were determined for both channels. The results are shown in table 5.
where E =total temperature range (-200"0 C to +200"0 C); YM total variation of the output variable; AX variation of the input variable; A Y variation of the output variable.
The operational parameters are shown in table 1.

Response of the two channels
The temperature in the thermostatted bath is recorded, the characteristic durations tl and t2 are determined from the recording of the indicial responses (see table 2). The values of the latency time, the pure delay and the static (2) If 7-/0 > 20, then the regulation must be Boolean in nature.
As the controller requires common setting parameters for both channels (for the integral and derived time constants), the best adapted regulation was considered for this experiment to be of the type PI: XPI-I 15; XPc 6; Ti (720 + 414)/2 567; Td 0 The overlapping segment, R, was arbitrarily set to zero, because no examples of application to hot/cold regulation were found in the literature. The introduction of these five parameters: XP4, XPc, Ti, Td and R in the controller, allowed the authors' first attempts at controlling the temperature in the reactor/crystallizer.
The indicial response of the system to a programmed set value was recorded. In a previous study [3], and in experiments with the calculated parameters, it was observed that the temperature in the reactor does not follow the set value and that considerable discrepancies appear when the process is run automatically; these discrepancies may reach 10C. To reduce these unwanted effects, it was decided to optimize the regulation parameters in this experiment.
Optimization of the PID parameters on the simplex method The aim here was to optimize the controller's parameters to obtain a value closer to the programmed value using the simplex method [4,5]. Starting with an experiment already conducted in the experimental domain, k other experiments are conducted to obtain a regular geometric shape called a 'simplex'. The responses for these k + other experiments are measured. The experiment with the least desirable result is discarded, and a new one is substituted, with co-ordinates obtained by calculating the symmetric point of the discarded point with respect to the centre of gravity of those remaining. The coordinates for the starting simplex, expressed with reduced variables, are given in table 6. The reduced and real coordinates are linked as follows: x i,j X ,j --Xi,j / xj (7) new point cannot be added because it would be situated beyond the allowed boundaries (out of the domain), the symmetric of the second worst point is used instead. The following notation is used for the points forming the simplex: FR Si -a Zos -t-b Z Ei (10) where ranged from to k + and ranged from to k; xi,j the real coordinate of the jth variable at the point i; Xl,j =the real coordinate of the jth variable at the starting point; Xi,j =the reduced variable of the jth variable at the point i; Axj the chosen incremental step for thejth variable.
Simplex evolution follows several rules: 1) Rule 1: The worst point is discarded and is replaced by its symmetric point with respect to the centre of gravity of the remaining points.
(2) Rule 2: When a better point appears in k+ consecutive simplexes, the value must not be overestimated, and it must be verified. This prevents endlessly producing a false optimal point. If a new experiment confirms the original value, the process is resumed; otherwise the faulty point is eliminated.
(3) Rule 3: When the evolution takes place on a ridge, i.e. the newest point is always the worst one, or when a  Arbitrary values, and 3, were chosen for the constant a and b; this was to ensure that the weight of each term in the expression would be equivalent. The function was chosen to favour regulation rather than tracking. The surface is determined as follows: the surface contained between the temperature curve and the line made by the set value drawn on a fixed-density paper is cut and weighted; the weight is converted in a surface expressed in square millimeters. The object of this optimization is to minimize the FR.

Variable choice
Five variables are likely to affect the quality of the regulation in some way: (1) The proportional width of the 'hot' channel, expressed in %.
(2) The proportional width of the 'cold' channel, expressed in %.
(4) The derived times (expressed in seconds). (5) The overlapping segment of the channels (in %). So, the system has five variables; the geometric representation of the simplex is therefore a six-vertex figure in five-dimensional space.  AXPH=5; AXPc=5; ATi=50; ATd=25; AR=20 Table 7 shows the reduced and real co-ordinates for the first simplex.

Boundaries
The applicable domain is limited by a number of constraints. First, the proportional bands of the 'hot' and 'cold' channels must be an integer between 0% and 100%; they have to be integer numbers above 10%. Second, the integral time constant must be an integer between 0 and 1999 s. Third, the derived time constant must be between 0 and 999 s; it must be an integer above 10. Finally, the overlapping of the bands must vary between -10% and +20% of the total measurement range, in increments of 0" 1%.

Results and interpretation
The initial simplex evolution can be traced in table 8, which shows the numerical values for all five coordinates and the value of the response function. Three values that proved useful for obtaining indications about the best moment to stop the evolution are also given: (1) The deviation between the responses at the best and worst points: Yw-YB.
(2) The value of the standard deviation (SD): In simplex 3, the new point, 8, remains the worst one, but because it presents a better response than point 4 whose symmetric it is, it is kept in the simplex, and the symmetric from the next to worst point, 3, is taken instead. The iterations continue normally until simplex 10. At this point, when the point 12 is discarded (Y12 330), point 16 is obtained but eliminated because its response is poorer (Y16 370). The next to worst point, 15, would not be discarded because this would lead to point 6, which has already been tested. The third worst point, 7 (Y7 302), was thus removed, and point 17 was obtained but also discarded because of its poorer response (f17 336). The fourth worst point, 14, was then discarded.
In simplex 12, as the new point is the worst, it is discarded; the evolution of the simplex is then stopped. The responses have been improved from 814 to 240. The value of the three criteria has also been significantly reduced: 78 instead of 482 for the deviation between the extreme responses, 32 instead of 165 for the standard deviation and 12 instead of 62 for the root-mean square deviation.

Optimization by the Nelder and Mead method
This derived method uses a variable-increment simplex [8][9][10]. If the symmetric becomes the best of all the points of the new simplex, then the direction of the last move is considered to be the most profitable. An extension is then applied. The following cases must be distinguished: If the response is worse than YB but better than Y2v, the process resumes without changing the size of the simplex. If the response is worse than Y2v, two cases are distinguished: (a) Ye is better than Yw: a contraction is applied at the side containing R and point Ca is added; (b) Ya is still worse than Yw: a contraction is applied at the side containing W and point Cw is added.
The coordinates of the points making up expanded (or contracted) simplex are given by the formula: (13) where "y 2 for E; -y 0"5 for Ca; "/= -0"5 for Cw.   For each simplex, the best vertex is expressed in bold Table 9 shows the evolution of the original simplex, with numerical values for the five sets of coordinates, and the value of the response function 'FR', including data about the three ending criteria. In order to build the second simplex, as the response of the symmetric point is the best of all, an expansion is applied, and retained. No further size alteration occurs until the fifth simplex. The sixth simplex was built from a contraction, as the symmetric point for point 5 (Y5 419) led to a better response (Y12 370), but it was the worst of this new simplex. Thus, a contraction was applied on point 12's side.
In the seventh simplex, the symmetric point for point 9 (Y9 308) was found not to be an improvement over the previous ones (Y5 338). The evolution of the simplex was discontinued at this point because: as with the initial method, the response was improved from 814 to 240 for the best point of the last simplex; the variation between the points of the last simplex is of the same order of magnitude: 68 (instead of 78); the SD has the same value: 31; the RMS was lower: 9 (instead of 12). So, 15 experiments were enough to obtain comparable results. The best point was obtained in the first expansion in the second simplex.
Optimization by the multi-move method This method [11] divides the experiments into groups of good and bad points rather than a single bad point and the rest. This could allow for a planning and lead to a reduced number of experiments. Table 10 gives the data used in the evolution of the original simplex. In the first basic simplex, the Hendrix method requires no alteration because the group of bad points only comprises point 1. In the second simplex, however, the group of bad points includes points 2, 3, 4 and 5 (see table 11).

Results and interpretation
The most important variation between the mean values is 140"8. Points 2, 3, 4 and 5 were therefore discarded. This process was repeated until point 16 was reached, with a response of Y16 362. If two groups of points were defined, then point 16 was the group of bad points. Consequently, point 16 was discarded; but then point 14 resulted. The simplex evolution was ended at that point because the same best response Y5--240 was obtained, although the variation between the points of the last simplex was higher (122 instead of 78). The SD and RMS were very close (39 for 32 and 15 for 12).
The results of this method, in the third simplex, were better than the original method because its 10th simplex had responses varying between 370 (Y16) and 260 (Y13).
Weighted centroid method ( WCM) This method moves away from the worst point towards the best one [12]. In order to achieve this the determination of the symmetric point takes into account a weighted centre of gravity (the weight of each point being a factor determined from the response at the point and the worst response). The co-ordinates for the new, weighted centre are: To avoid excessive degeneration of the simplex, the authors recommend limiting the value of 7, the ratio between the distance from the centre of gravity G to the point Gw and the distance from G to the best point B.
IIGw GII/IIB GII 3' ranges from 0 to 1. The highest advisable value 5 for this ratio is 0"3. If is lower than f or equal to it, the point GN is retained; otherwise, it is replaced by Gc, the revised centre of gravity, whose coordinates are given by XGc, Xa,j + (5(XB,j X,j) Gc, located between G and B, is not necessarily between G and Gw (unless only two variables are used for the optimization). For example, for a three-variable optimization, Gc is situated on an arc of a circle whose centre is G. In order to keep both the line of highest increase and the ponderation, the coordinates of point G are deter- The distance IJG-GII is equal to IIG-GcII, and the point G is situated between G and Gw. Results and interpretation Table 12 shows the evolution of the initial simplex. With the second simplex, because the ratio was higher than the limit 0"3, G was used as the centre of gravity. The evolution continued to the seventh simplex, where the symmetric point could no longer be considered because its response was much higher. The process ended there, with responses ranging from Y3 170 and Y9 275,  7 show the evolution of the three response criteria for each method with regard to each simplex. Figure 8 shows the evolution of the best responses; the weighted centroid method is clearly the best one.

Conclusion
A hot/cold controller order was used to regulate crystallization operations. The system was identified with a common method, the Broida method, and the parameters were obtained with the Ziegler-Nichols method. The paper shows that this empirical method will only allow a qualitative approach to regulation and that, in some instances, the parameters obtained are unreliable and therefore cannot be used to cancel variations between the set point and the actual values.
Optimization methods were used to determine the regulation parameters and solve this identification problem. These methods are the simplex approach and its derivatives (the Nelder and Mead method, multi-move method and the weighted centroid method). The best adapted was the weighted centroid method, which produced a very good response at the 13th point. The other methods yielded less satisfactory responses (the best response obtained in the others being 240). However, the coordinates for the points with this response were not identical; and 19 have the same response at 240, point 19 is better with the simplex method. It is also better than point 15 with multi-move and point 8 with Nelder-Mead, because the temperature with the parameters at point 19 never rises above 30"7C, while the other points lead to temperatures above 31 C. The most desirable characteristic for this regulation is a low deviation from the set value even if this implies longer stabilization delays. Even when the temperature is stable, deviations in the order of 0"2 C still occur.
The weighted centroid method produced good temperature regulation in a very short time. The PID parameters obtained were quite different from the parameters calculated. The Broida method is relatively unrefined and should be used in conjunction with an optimization method. The optimized PID parameters were used in the experiments and have given satisfactory results. If the set value is 30 C, the deviation is about 0.5 C; it rises to about C for a set value of 60 C.