A Fast Fourier Transform Approach for On-Line ... - ACS Publications

Jan 23, 1998 - A monitoring procedure incorporating the fast Fourier transform (FFT) technique is presented to identify the maximum closed-loop log mo...
0 downloads 0 Views 125KB Size
Ind. Eng. Chem. Res. 1998, 37, 1045-1050

1045

A Fast Fourier Transform Approach for On-Line Monitoring of the Maximum Closed-Loop Log Modulus Jun Ju and Min-Sen Chiu* Department of Chemical Engineering, National University of Singapore, Singapore 119260, Singapore

A monitoring procedure incorporating the fast Fourier transform (FFT) technique is presented to identify the maximum closed-loop log modulus (Lc,max) on-line. The proposed method addresses several potential problems present in the previous methods, which are (1) too many relay tests are required, (2) the frequency search range is confined to the third quadrant, and (3) the identified value of Lc,max cannot be used on-line to redesign the controller, if there is sufficient incentive to do so. Simulation results also show that the proposed method can identify Lc,max rather accurately as compared to the previous methods. Introduction Despite the recent advances in multivariable control design methods, the control scheme used in chemical process industries is still dominated by simple multiloop PID/PI controllers. Every plant has at least a few hundred control loops operating in multiloop control structure. Due to this high degree of automation, the proper design of controller parameters is crucial to achieve control objectives from product specifications to environmental regulations. However, it is well-documented that many existing multiloop PID/PI control loops perform poorly. This is because the design of the controller is normally based on a process model which is rarely a perfect representation of the actual process. Even if a perfect model is available, its parameters are uncertain, and often they depend on the operating conditions. Therefore, the controller’s performance is expected to deteriorate as the process dynamics change with time. Thus, it is important to have an efficient and reliable tool to determine on-line if the controller is performing according to the design specification. The goal of monitoring is to provide the information that can be used to assess the current status of the existing controller and to assist the plant operators and process engineers in deciding whether a redesign of the controller is necessary. However, most of the monitoring procedures reported in the literature are based on time domain techniques (A° stro¨m et al., 1986; Harris, 1989; Stanfelj et al., 1993) and the deterioration of the control system cannot be detected until a significant deviation in the controlled variable is observed. Recently, Chiang and Yu (1993) proposed a frequency domain monitoring procedure to identify on-line the maximum closed-loop log modulus (Luyben, 1990), Lc,max, for single-input single-output (SISO) systems. Lc,max is considered because its value can be used to gauge the compromise made between performance and robustness. As mentioned above, the dynamics of most chemical processes vary with time and therefore a certain tradeoff must be made between these two conflicting objectives in the controller design. Based on this, a +2n dB Lc,max is recommended as the design criterion for n × n multiloop control systems (Luyben, * To whom all correspondence should be addressed. Phone: (65)874 2223. Fax: (65)779 1936. E-mail: [email protected].

1986). A large number of tests show that this criterion yields a reasonable controller design. In order to evaluate Lc,max on line, two to three relay experiments are required in Chiang and Yu’s method, which may not be effective for a simple SISO system. In addition, their method has another two limitations, i.e., the frequency search range for Lc,max is restricted to the third quadrant and the identified value of Lc,max cannot be used in the redesign of the controller on-line. Ju and Chiu (1997) further extend their work to multiloop control systems. However, the number of relay experiments required in their method increases exponentially with the dimension of the multiloop control systems, e.g., 9 relay tests for a 2 × 2 system and 21 relay tests for a 3 × 3 system. Therefore, it can be time-consuming to apply their methods in practice due to the tedious procedure involved. The purpose of this work is to find a different monitoring procedure to address several potential problems present in the previous methods as described above. The recently developed fast Fourier transform (FFT) based frequency response identification technique (Hang et al., 1995) is capable of achieving this objective. The theory behind the FFT technique reveals that as many accurate frequency response points as desired can be obtained from one relay experiment. Because of the simplicity and accuracy of the proposed FFT-based monitoring procedure, it has a great potential in industrial application. Preliminaries This section presents the concepts of relay feedback and Lc,max and mathematical tools necessary for further developments of this paper. Relay Feedback Systems and Lc,max. A° stro¨m and Ha¨gglund’s relay tuning (1988) is a simple method for exciting a dynamic system and estimating useful process characteristics which can be used in the controller design. Figure 1 illustrates the relay feedback system in which a relay (N) is inserted in the feedback loop. The key idea is that many process systems will exhibit stable limit cycles when subjected to relay feedback. According to Luyben (1990), the maximum closed-loop log modulus, Lc,max, is defined as follows for multivariable control systems:

S0888-5885(97)00412-0 CCC: $15.00 © 1998 American Chemical Society Published on Web 01/23/1998

1046 Ind. Eng. Chem. Res., Vol. 37, No. 3, 1998

Figure 1. A° stro¨m-Ha¨gglund relay feedback system.

(

Lc,max ) 20 log max ω

|

|)

W 1+W

(1)

and

W ) -1 + det(I + GK)

(2)

where G is the process transfer function matrix, K is the controller transfer function matrix, and I is the identity matrix with the same dimension as the matrix GK. It can be shown from the stability analysis that the closer W approaches (-1, 0), the closer the feedback control system is to closed-loop instability. Therefore, a large Lc,max would indicate that the control system can tolerate a small amount of uncertainties. The identified Lc,max from the proposed method can be compared with +2n dB, which is the design criterion used in the n × n multiloop control system by Luyben (1986). It is then clear that Lc,max provides a good measure to assess the current status of the existing control system and to help in deciding whether or not there is a need to redesign the controller. Discrete Fourier Transforms and Fast Fourier Transforms. Suppose that a continuous function f(t) is sampled at the time interval of T and the magnitude of the sampled values is denoted by fm ≡ f(mT) (m ) 0, 1, 2, ...). If the function f(t) is nonzero only in a finite interval of time, that is, f(mT) ) 0 for m g L, the definition of the discrete Fourier transform (DFT) of the L points fm is given by (William et al., 1989) L-1

Fi ) F(ωi) ≡ T

∑ m)0

fme-2πjmi/L ) T

L-1

2πi LT

relay experiments which can be tedious in industrial application. In addition, their method has another two limitations: (1) the frequency search range is restricted to the third quadrant; (2) the identified value of Lc,max cannot be used in the redesign of the controller on-line. In this paper, the FFT technique is employed such that Lc,max can be obtained directly from only one relay experiment. As discussed in the previous section, the signal to be transformed by using FFT should be nonzero only in a finite interval of time. Hence, both the signals v(t) and u(t) shown in Figure 2 cannot be transformed to the frequency domain by directly using FFT. However, this problem can be overcome by decomposing v(t) or u(t) into the transient parts ∆v(t) or ∆u(t) and the periodic stationary cycle parts vs(t) or us(t) as (Hang et al., 1995)

-

i ) 0, 1, ..., L - 1 (3)

It is obvious that the more sampling points there are, the longer the computing time will be required. The efficient tool to calculate the DFT is the fast Fourier transform (FFT). There are different versions of FFT implementations available. The FFT routine used in this work is based on the algorithm written by N. Brenner (William et al., 1989). FFT-Based On-Line Monitoring of Lc,max SISO Control Systems. Figure 2 is the configuration of the monitoring system proposed by Chiang and Yu (1993). The objective of this monitoring scheme is on-line identification of Lc,max. For the purpose of monitoring, a series of relays are placed between the process and controller. In the monitoring mode, the controller output goes directly into the relay and the input of the process is the output of the relay. In the control mode, there is no relay and the controller output goes into process directly. Chiang and Yu devised a monitoring procedure to obtain Lc,max on-line for SISO systems. However, their method involves two to three

(4)

u(t) ) ∆u(t) + us(t)

(5)

where ∆v(t) and ∆u(t) will decay to zero when the system reaches the stationary oscillation. Now, the FFT technique can be incorporated into the identification of Lc,max, as will be shown in the subsequent discussion. Based on the monitoring scheme in Figure 2, the following equation holds.

i

i ) 0, 1, ..., L - 1

v(t) ) ∆v(t) + vs(t) and

fme-jω mT ∑ m)0

where the frequency ωi is defined as

ωi )

Figure 2. SISO monitoring scheme: m stands for monitoring mode and c stands for control mode.

∆V(s) + Vs(s) V(s) 1 )) G(s) K(s) (6) )N U(s) ∆U(s) + Us(s)

where V(s) and U(s) denote the Laplace transforms of v(t) and u(t), respectively. Similarly, ∆V(s), ∆U(s), Vs(s), and Us(s) denote the Laplace transforms of the transient parts and periodic parts, respectively. Again, let L denote the smallest integer such that both v((L - 1)T) and u((L - 1)T) have reached the stationary state or both ∆v((L - 1)T) and ∆u((L - 1)T) are approximately zero. Hence, the DFT of ∆v(mT) and ∆u(mT) can be computed at different frequencies with the standard FFT technique, i.e. L-1

∆V(jωi) ) T

∆v(mT) e-jω mT ) FFT(∆v(mT)) ∑ m)0 i

i ) 0 to L - 1 (7)

and L-1

∆U(jωi) ) T

∆u(mT) e-jω mT ) FFT(∆u(mT)) ∑ m)0 i

i ) 0 to L - 1 (8)

where the frequency ωi is defined as before. For the periodic parts, the following Lemma (Kuhfittig, 1978) can be used to calculate their frequency responses.

Ind. Eng. Chem. Res., Vol. 37, No. 3, 1998 1047

Lemma. Suppose fp(t) is a periodic function with period Tc, i.e., for all t, fp(t) ) fp(t + T). Let Fp(s) denote the Laplace transform of fp(t); then

Fp(s) )

1 1 - e-sTc

∫0T fp(t) e-st dt c

(9)

Hence, Vs(jω) and Us(jω) in eq 6 can be computed at each frequency ωi using digital integral as follows:

Vs(jωi) )

Nc

T

vs(mT) e-jω mT ∑ m)0

(10)

i

-jωiTc

1-e

and

Us(jωi) )

Nc

T

us(mT) e-jω mT ∑ m)0

(11)

i

-jωiTc

1-e

where Nc ) (Tc - T)/T. Hence, at frequency ωi, eq 6 becomes

-

1 N

∆V(jωi) + )∆U(jωi) +

T -jωiTc

1-e

T -jωiTc

1-e

Figure 3. Monitoring scheme for a n × n multiloop control system: m stands for monitoring mode and c stands for control mode.

the process. Repeat until all possible pairs of control loops are tested.

Nc

∑ vs(mT) e-jω mT

l

i

m)0

)

Nc

∑ us(mT) e-jω mT i

m)0

G(jωi) K(jωi) (12) It is now clear that eq 12 can compute the frequency responses of G(s) K(s) at L frequencies in the frequency interval [0 2π(L - 1)/LT]. Subsequently, these frequency responses can be used in eq 1 to calculate Lc,max. Thus, one can use only one relay experiment to identify Lc,max on-line instead of several relay experiments required in Chiang and Yu’s method (1993). Moreover, it should be noted that with the FFT based approach, the frequency search range for Lc,max can extend outside the third quadrant. Consequently, the proposed method relaxes one of the limitations inherent in Chiang and Yu’s method (1993). Equation 12 also suggests that the process G can be identified in closed-loop fashion since the information of the controller is known. This implies that the monitoring information can be used to redesign the controller on-line if the identified Lc,max deviates from the design criterion, e.g., +2 dB according to Luyben (1986). These features will be illustrated further in the Examples section. Multiloop Control Systems. The framework for online identification of Lc,max for a n × n multiloop control system is illustrated in Figure 3, where G ) [gij]i,j)1-n and K ) diag[ki]i)1-n. The theoretical analysis of various relay tests required in this monitoring system is given by Ju and Chiu (1997) and will not be retraced here. With a slight modification to their result, the following gives the monitoring procedures of identifying Lc,max for n × n multiloop control systems via the FFT method. Step 1. Place a relay between the controller and the process in one control loop while keeping all other control loops in control mode. Repeat until n control loops take turns in the monitoring mode. Step 2. Select a pair of control loops and place a relay between each controller in the chosen control loops and

Step n - 1. Select n - 1 control loops and place a relay between each of the n - 1 controllers in the chosen control loops and the process. Repeat until all possible combinations of n - 1 control loops are tested. Step n. For all n control loops, place a relay between each controller and the process. It is noted that the proposed FFT based monitoring scheme only requires one test for each relay experiment, whereas three tests are involved in the previous method. In other words, the number of relay tests can be reduced from 9 to 3 for 2 × 2 systems and from 21 to 7 for 3 × 3 systems. This means that the experiment time is shorter if the FFT-based monitoring method is employed to identify Lc,max. In the subsequent sections, the identification procedures for 2 × 2 and 3 × 3 multiloop control systems via the FFT method are discussed in detail. 2 × 2 Multiloop Control Systems. There are two steps required in monitoring 2 × 2 multiloop control systems. Step 1. A relay is placed between the controller and the process in loop 1. After one relay is performed at t ) 0, the controller output v1(t) and process input u1(t) of loop 1 are recorded until v1(t) reaches the stationary state. Subsequently, the controller output v2(t) and process input u2(t) of loop 2 are recorded when loop 2 is switched to the monitoring mode. N1 and N2 are used to denote the relays applied in loop 1 and loop 2, respectively. Step 2. Two relays (N3 and N4) are placed between the controller and the process in both loops. The above analysis discussed for SISO systems can now readily be applied to evaluate the relays N1 to N4 in steps 1 and 2. For example, the following two relations result from the application of eq 12 to the signals v1(t) and u1(t), v2(t) and u2(t), respectively.

1048 Ind. Eng. Chem. Res., Vol. 37, No. 3, 1998

1

-

)-

N1

V1(jωi) U1(jωi)

∆V1(jωi) + ∆U1(jωi) + 1

-

(

Lc,max ) 20 log max

)

)-

N2

V2(jωi)

Nc1

T1 1 - e-jωiTc1

1

(13)

Nc1

T1

∑ u1s(mT1) e m)0

-jωimT1

-jωiTc1

1-e

F2(ωi) )

∆U2(jωi) +

(

Nc2

T2

v2s(mT2) e-jω mT ∑ m)0 i

1 - e-jωiTc2

2

u2s(mT2) e-jω mT ∑ m)0 i

-jωiTc2

1-e

A1 )

(

A2 )

(

A3 )

(

(14)

Nc2

T2

2

Nc2 ) (Tc2 - T2)/T2

Likewise, relays N3 and N4 can be obtained in a similar fashion. Therefore, the describing functions of all relays required in the monitoring scheme can be used to compute Lc,max by the next equation (Ju and Chiu, 1997).

( |

Lc,max ) 20 log max ωi

F1(ωi)

1 + F1(ωi)

|)

(15)

where

[(

1-

)(

1 N1

|)

(17)

)( )( )( ) )]/[( )( )( ) ] )( )( ) { ( )( ) } )( )( ) { ( )( ) } )( )( ) { ( )( ) }

A1 - A2 - A3 (18)

Nc1 ) (Tc1 - T1)/T1,

[(

[( )(

)

where T1 and T2, Tc1 and Tc2 are sampling times and stationary oscillation periods of stable limit cycles in loop 1 and loop 2, respectively. Also define

)(

)(

)( )]/ ) ( ) ( )]

1 1 1 1 111N1 N2 N3 N4 1 1 1 1 1- 11- 1× N2 N1 N4 N2 1 1(16) N3

F1(ωi) ) -1 -

1 + F2(ωi)

1 1 1 1 -1 -1 -1 -1 × N1 N2 N3 N10 1 1 1 1 1 -1 -1 -1 -1 -1 N11 N12 N1 N2 N3 -1 -

U2(jωi) ∆V2(jωi) +

F2(ωi)

where

v1s(mT1) e-jω mT ∑ m)0 i

ωi

|

1-

) (

)(

3 × 3 Multiloop Control Systems. There are three steps required in monitoring 3 × 3 multiloop control systems. Step 1. A relay (N1) is placed between the controller and the process in loop 1, i.e. loop 1 is in the monitoring mode. Subsequently, two relay experiments (N2 and N3) are conducted when loop 2 and loop 3 take turns in the monitoring mode. Step 2. Two relays (N4 and N5) are placed between the controller and the process in both loop 1 and loop 2, i.e., only loop 3 is in the control mode. Next, relay experiments (N6 and N7) are applied to both loop 1 and loop 3, while loop 2 is put on the control mode. Finally, relay experiments (N8 and N9) are conducted in both loop 2 and loop 3. Step 3. Three relays (N10, N11, and N12) are placed between the controller and the process in all three loops. Again, the frequency responses of these relays can be computed via the FFT method. Thus, the next formulation can be employed to compute Lc,max (Ju and Chiu, 1997).

1 1 1 -1 -1 -1 × N2 N3 N10 1 1 1 1+ -1 -1 N1 N11 β1 1 1 1 -1 -1 -1 × N1 N3 N11 1 1 1 1+ -1 -1 N2 N12 β2 1 1 1 -1 -1 -1 × N1 N2 N12 1 1 1 1+ -1 -1 N3 N10 β3

β1 )

[( (

)( )( )( )]/[( ) ) ( )( ) ( )( )]

1 1 1 1 1 -1 -1 -1 -1 -1 × N1 N2 N4 N5 N1 1 1 1 1 1 -1 -1 -1 -1 -1 N2 N2 N4 N1 N5

β2 )

[( (

)( )( )( )]/[( ) ) ( )( ) ( )( )]

1 1 1 1 1 -1 -1 -1 -1 -1 × N2 N3 N6 N7 N2 1 1 1 1 1 -1 -1 -1 -1 -1 N3 N3 N6 N2 N7

β3 )

[( (

)( )( )( )]/[( ) ) ( )( ) ( )( )]

1 1 1 1 1 -1 -1 -1 -1 -1 × N1 N3 N8 N9 N1 1 1 1 1 1 -1 -1 -1 -1 -1 N3 N3 N8 N1 N9

Examples SISO Control Systems. In order to illustrate the proposed FFT monitoring procedure and the improvements over the previous method (Chiang and Yu, 1993), the following control system is considered.

G(s) )

37.3e-10s , (7200s + 1)(2s + 1)

(

K(s) ) 10.54 1 +

1 (19) 39.96s

)

To apply the FFT method to the monitoring scheme, an ideal relay experiment is performed and the data of v(t) and u(t) are logged. Figure 4 shows the simulation results and the decomposition of two signals u(t) and

Ind. Eng. Chem. Res., Vol. 37, No. 3, 1998 1049

Figure 6. Closed-loop log modulus versus frequency.

control system, the identified Lc,max is 6.814 dB. Therefore, it is possible to redesign the controller to satisfy the +2 dB design criterion. To do so, the process dynamics (GFFT) are estimated first from eq 12 and a detuning factor f (Luyben, 1986) is introduced such that

K′ )

Figure 4. Simulation results and the decomposition of (a) v(t) and (b) u(t) from the relay experiment.

10.54 1 1 1+ , f f 39.96s

(

)

e-2s , (s + 0.5s + 1)(s + 1) 2

(

K(s) ) 0.38 1 +

v(t). The resulting Lc,max (shown in Figure 5) is 6.814 dB, which is very close to the actual value, 6.754 dB. By comparison, Chiang and Yu’s method can only obtain a very rough estimate of Lc,max, i.e., 10.438 dB. Therefore, the proposed method is much more accurate than their method for this example. The sensitivity of the proposed method to the measurement noise is also studied. If white noise is inserted between the controller K and relay element N and the ratio of the noise to signal is set to be 10%, the identified Lc,max is equal to 7.535 dB, which differs from the result of the noise free situation by 10%. As discussed previously, the proposed method has online tuning capability. For example, in this SISO

(20)

Thus, with the identified GFFT, Lc,max can be adjusted on-line to be equal to +2 dB through designing f. The calculated f from the proposed FFT method is 7.3, as compared to 7.8 from the true process. Last, it should be noticed that the on-line search method for Lc,max proposed by Chiang and Yu (1993) is confined to the third quadrant. Therefore, if the Lc,max is located outside the third quadrant, e.g., in the second quadrant, their method will fail. By contrast, the search space for Lc,max in the proposed FFT-based method is not limited in the third quadrant. To illustrate this feature, we choose a control system as follows:

G(s) )

Figure 5. Closed-loop log modulus versus frequency.

f>1

1 + 1.34s (21) 1.30s

)

It can be computed that the actual value of Lc,max for this control system is 3.319 dB, which occurs in the second quadrant. Figure 6 demonstrates that the resulting Lc,max from the proposed method is 3.354 dB. Therefore, from a practical point of view, the proposed method is more versatile than Chiang and Yu’s method. Multiloop Control Systems. The distillation column studied by Ogunnaike et al. (1983) is used to test the proposed monitoring procedure, where

[

G(s) ) 0.66e-2.6s 6.7s + 1 1.11e-6.5s 3.25s + 1 -34.68e-9.2s 8.15s + 1

-0.61e-3.5s 8.64s + 1 -2.36e-3s 5s + 1 46.2e-9.4s 10.9s + 1

-0.0049e-s 9.06s + 1 -0.01e-1.2s 7.09s + 1 0.87(11.61s + 1)e-s (3.89s + 1)(18.8s + 1)

]

1050 Ind. Eng. Chem. Res., Vol. 37, No. 3, 1998

Nomenclature DFT ) discrete Fourier transform FFT ) fast Fourier transform Fi ) discrete Fourier transform of the L points fm fm ) sampled value of f(t) at t ) mT G ) process transfer function K ) controller transfer function L ) number of sampled points Lc,max ) maximum closed-loop log modulus N ) describing function of relay SISO ) single-input single-output T ) sampling interval MIMO ) multi-input multi-output u(t) ) process input us(t) ) periodic stationary cycle part of u(t) v(t) ) controller output vs(t) ) periodic stationary cycle part of v(t) Figure 7. Closed-loop log modulus versus frequency.

[

K(s) ) 1 1.509 1 + 35.26s

0

0

1 -0.295 1 + 38.7s

(

0

)

0

(

0

)

0

]

1 2.629 1 + 14.211s (22)

(

)

The resulting Lc,max obtained via the proposed method is 4.435 dB (Figure 7), as compared to the actual value of 4.346 dB. Conclusions A FFT-based monitoring procedure to on-line identify Lc,max is proposed. Compared with the previous methods (Chiang and Yu, 1993; Ju and Chiu, 1997), it has several advantages. First, the method is very effective and accurate, and the testing time is shortened greatly. This is because only one relay experiment, instead of two to three relay experiments in Chiang and Yu’s method, is involved for SISO systems, and for MIMO systems, only one test is required for each relay experiment in the n-step monitoring procedure. Next, the frequency search range for identification of Lc,max is extended for SISO systems, whereas the previous method is restricted to the third quadrant. In addition, the identified information for Lc,max can be used on-line to redesign a new controller. All these features will make this FFT-based monitoring procedure very attractive in industrial application. Acknowledgment The authors thank Bing Wang for his assistance with some of the calculations.

Greek Symbols ∆u(t) ) transient part of u(t) ∆v(t) ) transient part of v(t) ωi ) frequency (2πi/LT), i ) 0, 1, ..., L - 1

Literature Cited A° stro¨m, K. J.; Ha¨gglund, T. Automatic Tuning of PID Controllers; Instrument Society of America: Research Triangle Park, NC, 1988. A° stro¨m, K. J.; Anton, J. J.; A° rze´n, K. E. Expert Control. Automatica 1986, 22, 277. Chiang, R. C.; Yu, C. C. Monitoring Procedure for Intelligent Control: On-Line Identification of Maximum Closed-Loop Log Modulus. Ind. Eng. Chem. Res. 1993, 32, 90. Hang, C. C.; Wang, Q. G.; Cao, L. S. Self-Tuning Smith Predictors for Processes with Long Dead Time. Int. J. Adapt. Control Signal Process. 1995, 9, 255. Harris, T. J. Assessment of Control Loop Performance. Can. J. Chem. Eng. 1989, 67, 856. Ju, J.; Chiu, M. S. Relay-Based On-Line Monitoring Procedures for 2 × 2 and 3 × 3 Multiloop Control Systems. Ind. Eng. Chem. Res. 1997, 36, 2225-2230. Kuhfittig, P. K. F. Introduction to the Laplace Transform; Plenum Press: New York, 1978. Luyben, W. L. Simple Method for SISO Controllers in Multivariable Systems. Ind. Eng. Chem. Process Des. Dev. 1986, 25, 654. Luyben, W. L. Process Modelling, Simulation, and Control for Chemical Engineers; McGraw-Hill: New York, 1990. Ogunnaike, B. A.; Lemaire, J. P.; Morari, M.; Ray, W. H. Advanced Multivariable Control of a Pilot Plant Distillation Column. AIChE J. 1983, 29, 632. Stanfelj, N.; Marlin, T. E.; MacGregor, J. F. Monitoring and Diagnosing Process Control Performance: The Single-Loop Case. Ind. Eng. Chem. Res. 1993, 32, 301. William, H. P.; Brian, P. F.; Saul, A. T.; William, T. V. Numerical Recipes: The Art of Scientific Computing (Fortran Version); Cambridge University Press: Cambridge, U.K., 1989.

Received for review June 9, 1997 Revised manuscript received November 10, 1997 Accepted November 11, 1997 IE970412P