Particle size distribution soft-sensor for a grinding circuit

Particle size distribution soft-sensor for a grinding circuit

POWDER TECI'WOLGGY ELSEVIER Pow,ler Technology99 (1998) 15-21 Particle size distribution soft-sensor for a grinding circuit A. Casali ~'*, G. Gonzal...

792KB Sizes 21 Downloads 310 Views

POWDER TECI'WOLGGY ELSEVIER

Pow,ler Technology99 (1998) 15-21

Particle size distribution soft-sensor for a grinding circuit A. Casali ~'*, G. Gonzalez b, F. Tones c, G. Vallebuona a, L. Castelli c, P. Gimenez ~ aMining Engineering Department, Universit)"of Chile, P.O. Box 2777, Santtago, Chile hElectrical Engineering Department, Universityof Chile, Santiago, Chile CODELCO-Chile, Andina Dwislon, Santiago, Chile

Received 14 November 1997; received in revised form I May 1998

Abstract Particle size measurement, specially in product streams, is fundamental to control of a grinding circuit. In this paper, the design of a softsen,;or in an actual grinding plant is addressed, with the purpose of increasing the availability of such measurement. The system is designed to substitute the momentary unavailability of the real sensor, either because it has failed o!"it is in calibration or because it is time shared by several parallel grinding lines. An ARMAX model structure is determined using the stepwise regression method. Starting with a list of the presumably correlated candidate variables, the method selects, one by one and in a systematic manner, the variables that best model the measurement to be replaced. Not only the single measurements are included in the list of candidates, but also combinations of such measurements having physical significance. A phenomenological model has been built with the purpose of determining these combinations. It turns out that, in most of the cases, the components of combined measurements were selected by the method, instead of the single measurements. Actual grinding plant data is used to determine the model structure, to estimate the model parameters and to test the predictive capability of the soft-sensors developed. Good results are obtained during periods of up to 24 h of the particle size soft-sensor operation. © ! 998 Elsevier Science S.A. All rights reserved. Keyword~: Soft-sensor; Panicle size distributton; Grinding

1. Introduction Adequate control of a grinding circuit usually relies on the availability of particle size distribution measurements on-line in the product streams. As a consequence of the usually harsh environment, the momentary unavailability of a sensor is possible either because it has failed or because it has been removed for maintenance or calibration. This is also the case when the sensor is time shared by several parallel grinding lines. When the sensor becomes unavailable the control of the involved variable, either automatic or manual, is lost, with the ensuing negative operational and economic consequences. Soft-sensors, or virtual sensors, have been used to provide and estimate the missing sensor signals, based on the relation between these signals and other measurements [ 1-5 ]. This relation is incorporated in a model having as input the corresponding related measurements. While the real sensor is in operation, the parameters of the soft-sensor model may be estimated on line (e.g., using a recursive least squares method). When the sensor is unavailable, the model output * Corresponding author.

becomes the soft-sensor prediction for the missing measurement. Soft-sensors are implemented by software and in many cases are based on ARMAX-type models [6]. In this paper, the design of a soft-sensor for the percentage of material with a dimension of + 65 mesh in the product of a section of the CODELCO-Andina grinding plant is addressed. In particular, the procedure for the determination of a model structure using the stepwise regression method is improved [7] through the inclusion of components having physical significance in the list of candidate components to be used as inputs to the soft-sensor model. This inclusion of the plant model knowledge in the components of the ARMAX model turns the original black box model into a kind of "grey box' model.

2. Derivation of soft-sensor models A diagram including instrumentation of one of the grinding circuits for which the particle size distribution soft-sensor has been designed and tested is given in Fig. I. The available online measurements, considered in this work, are fresh ore feed rate, rod mill water feed rate, rod mill power, sump water

0032-5910/98/$ - see front matter © 1998ElsevierScience S.A. All rights reserved. Pll S003 2-5 9 ! 0 ( 98 ) 00084-9

16

A, C¢zsaliet al. / Po~vder Technology 99 (lO98) 15-21

i

al t-:

T

i

Fig. I. Grinding circuit flowsheet. MF, fresh ore feed rate; Wl, rod mill water t'~d rate; JR, rod mill power; W2, s~mp water addition rate; PS, sump pump speed; SC, solid concentration by weight; JB, hall mill power: 165, percentage over 65 m~sh.

addition rate, sump pump speed, solid,; concentration, ball mill power and percent over 65 mesh. Apart from these variables, there are others that were not considered for this work. Such is the case of the sump level which was not considered because it was under local control. Other measurements such as the hydrocyclones feed pressure and volumetric flow rate were discarded due to their high correlation with each other and with the pump speed. Since the variable f65 model is used for the whole combined grinding section (see Fig. ! ), it was considered convenient to represent the circuit by the equivalent grinding circuit shown in Fig. 2. For this global circuit, measurements of the same variable for each one of the bull mill grinding lines have been grouped together. This is done by adding up the mill powers (JB), the water addition rates (W2) and the pump speeds (PS), and by averaging the solids concentrations (SC). The problem of finding a soft-sensor model for a plant may be separated into two stages. The first stage, structure determination, is performed off-line, using plant data, and its objective is to determine which are the components in the

ARMAX model that contribute most significantly to the modeling of the desired variable, f65 in this case. The other stage is the parameter estimation stage which is done on-line while the soft-sensor is installed in the plant. The objective here is to adjust the parameters of the components so that the mean square estimation error or residue between the actual measurement and the value estimated by the soft-sensor is minimized. 2. !. Candidate components

The candidate components for an ARMAX soft-sensor model fall into two groups, (i) single measurement and (ii) combined measurement. The on-line measurements available in the grinding plant likely correlated with the percentage over 65 mesh must he considered, including the auto-regressive terms i.e., measurements delayed in time. In some cases, as was already explained, these on-line measurements have been grouped together before being included in the list of possible components. In order to include combinations of measurement variables having physical meaning, process models, phenomenological considerations and knowledge of the grinding circuit have been taken into account. The phenomenological mathematical model developed in order to lind the combinations of measured variables having physical meaning, is based on (!) three well-known hydrocyclone models: (i) the slurry volumetric flow split equation [8l: (ii) the water split equation and the cut size equation [91" and (iii) the corrected efficiency equation [10]. (It) The population balance approach for grinding [ I 1,12]. (Ill) All possible solid (total and per size) and liquid mass balance equations. ( IV ) The linearization in a truncated Taylor series of all components that were lbund to be nonlinear in the parameters. The details regarding the development of the model are planned to be presented in a future paper. From this model, more than a hundred combined components relating the variable to model, i.e., f65, with some of the available inputs have been found. All these have been

®

Fig. 2. Equivalent grinding circuit. W2s, sum of the sump water addition rates; PSs, sum of the sump pump speeds; SCa, average solid concentration by weight; JBs, sum of the ball mill powers.

A. Casali et al. / Powder Technology 99 (1998) 15-21

included in the components candidate list, together with the single variable components. Considering the process phenomena and knowing the layout of the actual sensors, possible delays for each input variable of zero, one, two and three sampling times were assigned. With all possible combinations the list of candidate components reached a total number of over three thousand. 2.2. Experimental

Striving to cover the normal fluctuations in all the measured variables related with the output, an experimental test was performed. The plant was excited in order to be able to determine dynamic ARMAX models for the soft-sensor. The set points of MF, W! and W2s controllers were varied in a pseudo-random binary manner, as shown in Fig. 3. In Fig. 4 the fluctuations of the other measured variables during the test are presented. The data corresponding to all variables have been statistically normalized. 2.3. Stepwise regression

The determination of model structure is based on a stepwise regression method [ 13 ] that finds the set of components 2

--

o

!7

from a list of candidate components that give the best ARMAX-type model for the desired variable, in this case the percentage + 65 mesh in the product section of the grinding plant. Starting with the list of candidate components, the procedure begins by selecting the one that has the highest correlation with the output. In this way a first partial model is determined. The residue of this model (using a simple linear regression coefficient) is correlated with the remaining candidate variables. In each of the next steps, the component to be included is the one giving the highest partial correlation with the previous step residue, calculated after a multilinear regression is performed with the variables previously included. With this procedure, the model is increased by adding one component at a time, such that the added component is the one that contributes the greatest improvement in the goodness of fit to the model. Before each new component is included in the model, all the included components are tested for their statistical significance. Since a prediction error estimator is used to estimate the coefficients, 0, the estimates have an estimated prediction error variance, o'(0). The ratio between this variance and the coefficient, for each component included, is used to decide the inclusion of the new term. If for any component the ratio (see Eq. ( 1 ) ) exceeds a given value (0.3 in this case), then it is nt~, included.

%%2

O,

-2,

<0.3

(I)

.q

~

2 o~

t-

This procedure is repeated until a component is either deleted or included in the determined model. In each step an F (Fisher) test is performed to decide if the procedure has been completed. If the F value, calculated as shown in Eq. (2), is lower than I !, then the previous step was the last one and the model structure determination is complete.

Wl

-2'

2 0 "2 ' 0

t

MF

b0

100

lbO

700

250

r:

~"C mr. eS

R -RL,

*

(2)

Fig. 3. Mampulatedvariables during a persistent exc,ation test.

t"65

,

0 -3 c



o4 I

JBs -

13 N

[ 7

3

SCp

0"i. -3 3~ !',, 0 4

-31 o

.

PSs

,'

T

50

......

100

1£0

200

Samples

Fig. 4. Other measured variablesduring a persistent excitation test.

250

18

A. Casali et al. / Powder Technology 99 (1998) 15-21

where Riz is the multiple linear regression coefficient in the step i and dfindicates the degree of freedom. 2.4. Model structure

With the procedure applied to a first set of 240 data (previously normalized) obtained from the experimental test, the following structure was determined (Model A): f65(t)=kl f 6 5 ( t - l)+k2 J B s ( t - l)[SCa(t)] 3

3.2. Model C (JBs measurement missing)

Models A and B include the sum of the ball mill powers (JBs) as a common input variable. To overcome the unavailability of this measurement (when one or more of the JB measurements are unavailable) the stepwise regression method was applied, excluding components depending on JBs from the candidate component list. The result is model C, given by f65(t)=kt f 6 5 ( t - 1)+k 2 SCa(t)

[Wl(t)l 2 +k3 J B s ( t - 3 ) S C a ( t - 3 ) + k 4 M F ( t _ 2 ) P S s ( t - 3 )

(5)

[Wl(t)l 2 +k3 M F ( t - 2 ) PSs(t-3) (3)

where t is the sampling time and k, the parameters to be estimated as the operating point changes. It should be noticed that only the single component f 6 5 ( t - I), i.e., of the autoregressive type, was selected by the stepwise regression method. Two different data sets were used to test the robustness of the structure, both having 240 samples. In both cases the same structure was determined.

In order to ensure the availability of the soft-sensor in case of failure or unavailability of one of the measurements affecting the components of model A, alternative soft-sensor models have been determined. Ti~ey have been found assuming that one or more of the measurements affecting the compo. nents in model A is not present. A family of models then results so that a soft-sensor system is able to continue predicting the missing particle size sensor measurement, even if one of the input measurements is not present. The structure determination for the alternative softsensor models was accomplished following two different approaches: (i) truncation of a model and ( ii ) the use of the stepwise regression procedure with a reduced candidate components list. 3. I. Model B (WI or MF or PSs measurements missbzg)

The first alternative model was obtained by truncation of the main model presented in Eq. (3), removing the last term of this equation which does not contribute too much to the goodness of fit. The resulting new model shown in Eq. (4) does not include the rod mill water feed rate ( W l ) , the sum of the sump pump speeds (PSs) orthe fresh ore feed rate (MF) as inputs, and is given by

+k3 J B s ( t - 3 ) S C a ( t - 3 )

This alternative model was obtained by truncation of model C given hy Eq. (5), removing the last term of this equation which does not contribute too much to the goodness of fit. The resulting new model, shown in Eq. (6), includes only the average solid concentration (sr'a ~ and the auto-regressive term, and is given by f65(t)=k~ f 6 5 ( t - I )+k., SCa(t)

3. Alternative models

f65(t)=k~ f 6 5 ( t - 1)+k2 J B s ( t - 1)[SCa(t)] 3

3.3. Model D (WI or MF or PSs or JBs measurements missing)

(4)

(6)

3.4. Model E (WI or SCa or PSs or JBs measurements missing)

All the previous models include the average solids concentration (SCa) as a common input variable. To overcome the unavailability of this measurement the stepwise regression method was applied, excluding components depending on SCa from the candidate component list. In this case, after selecting one phenomenological component, the method rejects all other components, except the auto-regressive term, f 6 5 ( t - ! ), so the resulting model is f65(t)=k~ f 6 5 ( t - l ) + k ,

MF(t-I) "MF(t- I )+W2s(t)

(7)

In summary, Table ! gives the soft-sensor availability when an input signal, of those integrating the models, is missing. Table I Soft-sensor availability Missing signal

Available models

MF Wl W2s PSs SCa JBs

B and D B, D and E A. B, C and D B, D and E E C, D and E

A. Casali et al. / Powder Technology 99 (1998) 15-21

4. Models performance

!9

4O Exper'i..'rent al

35

Mode, A

I

To examine the performance of all the models developed, the data from the whole experimental campaign were used. The campaign consisted of a first period of 8 h with a persistent excitation test, immediately followed by a 30 h ( 1800 samples) period under normal operating conditions.

3oi

~5 lo4

I

4.1. Parameter estimation

5 ~

To estimate the parameters of all the models, three groups of data from the persistent excitation part of the experimental campaign were used: one with 60 samples, another with ! 20 samples and one with 240 samples. In this estimation procedure the unnormalized form of the models (adding the constant ko ) and unnormalized data (as they are measured) were used. Table 2 ~hows the results obtained in the parameter estimation stage. It may be seen that the correlation coefficient R 2 for all models are similar, when the set of 12e md 240 samples are used. Besides, the parameters did not change too much fer these two sets, with the exception of ko. In contrast, the results for the set with 60 samples presented a lower correlation and the parameters are quite different from the ones obtained in the 120 and 240 samples sets. Therefore, it was decided that a period of 120 samples was good enough for the parameter estimation of all the models considered. 4.2. Predictive capabili~, test

To test the quality of all five soft.-sensor models, data from the experimental campaign were used. The objective was to check the predictive capability of each model with additional data and on a long run test, using the model structures previously determined.

0 .........

0

.

"

500

1000

1500

2000

2500

Samples

Fig. 5. Predictive p e r f o r m a n c e o f model A.

The parameters of each one of the models were estimated with the first set of 120 samples (unnormalized). With each model, and using the parameters thus determined, the corresponding output (f65) was then predicted for the next peri,,d of over 2000 samples and compared with the actual f65 measurements (since the real sensor had not actually failed). The f65 measurements and the prediction obtained with each model, during the complete 37 h sampling period (2220 samples), are shown in Figs. 5-9, corresponding to models A, B, C, D and E, respectively. As should have been foreseen, the performance of model A is better than the performance of all other models. However, the estimation quality of all of them is, at least during some periods, quite acceptable. In any case it is better to have a prediction with 6% error than not having any measurement at all once the real sensor has failed or been removed. The average error during the first 550 samples (at this point an actual failure of the real sensor occurred) was less than 2% for all the models. With respect to the validity of using constant values for the previously determined model parameters, it can be seen from Figs. 5-9 that there was no obsolescence of the parameters.

Table 2 P a r a m e t e r estimation and g o o d n e s s o f fit for different data sets ( u n n o r m a l i z e d data)

Model

Samples

R2

Parameters

I

6()

0.86 !

27.8

0.78

8.2e - 9

- 1.0e - 4

!

120

0.905

7.3

0.72

9.9e - 9

- 6.2e - 5

0.013

I

240

0.916

7.4

0.69

9.4e - 9

- 5.6e - 5

0.009

2

60

0.847

19.5

0.83

7.3e - 9

- 9.8e - 5

-

2

120

0.897

7.9

0.76

9.4e - 9

- 6.0e - 5

-

2

240

0.91 !

4.9

0.71

9.1e-9

-4.4e-

-

3

60

0.826

- 25.3

0.83

0.40

-

-

3

120

0.889

- 27.6

0.70

0.47

-

-

3 4

240 60

0.903 0.832

- 2 ! .7 - 27. !

0.69 0.80

0.39 0.42

0.011

-

4

! 20

0.896

- 29.7

0.66

0.49

0.0 ! 3

-

4

240

0.907

- 24. I

0.67

0.42

0.008

-

5

60

0.808

- 2.7

0.90

89.04

-

-

5

120

0.868

0.3

0.82

65.22

-

-

5

240

0.890

0.4

0.82

62.03

-

-

5

0.017

20

A. Casali et ai. / Powder Technology 99 (1998) 15-21

30 4.,-'i

25

r

,!

25

-.',

/ .,~.. . .

"t

"

t

~2o

E ~15

I

glO

°; 0

sJ 0

IF ~

t

.

0

500

.

.

Experimental

.

.

.

.

.

1000

.

I

- Model B

.

.

.

.

1500

2000

~

0

2500

-

0

500

.

.

.

Exper,mental - -- Voael 0 .

I000

Samples

t500

2000

2500

Samples

Fig, 8. Predictive perlbrmance of model D.

Fig, 6. Predictive performance of model B. 30 30

,

~

- .

25

/

t .,,

,

.

.

.

.

.)

,

/

2

5

~

~

~:,o E

01o,

l

b> °13-

I

5

_

5 .

I

i

I

_

Ex#er i'*'en" aE

[xper~mental ~ McOel C O+ C

~

o

.

5oo

.

.

.

.

.

.

.

"o0o

.

.

15oo

.

.

.

.

2000

2503

Samples

Fig. 7. Predictive peflbrmance of model C.

Mo<:el E

-

500

o

1000

1500

2030

2500

Samples

Fig. 9. Predictive performance of model E.

4.3. Reliability This is a good condition, since the parameters cannot be updated because the actual measurement is supposedly not available. It can be noticed, in Fig. 6 for example, that the prediction in the interval 1500-2000 is better than in the interval 1000- !500. However, some of the model parameters are presumably no longer adequate during some periods, because the plant characteristics affecting the model may have changed, mainly due to changes in the operating point or changes in the disturbances. When the plant returns to the origi~'d operating point, the fixed parameters turn out to be adequate again and the soft-sensor prediction error is once again small. Possible improvements of the soft-sensor prediction by using optimal estimates of the soft-sensor model parameters (once the actual sensor is no longer available) has been considered [ 3 ]. Shortly after the sensor becomes unavailable it may be adequate to freeze the model parameters estimated up to that moment and use them in the prediction stage. However, as time changes from when the soft-sensor began its operation, the optimal choice would be constant parameters equal to their respective expected values. Optimal values for the soft-sensor parameters could begin at the values they had at the moment of failure and change exponentially to their expected values. Since the results in the present case are quite acceptable, it was deemed unnecessary to use a more elaborate procedure.

Two cases must be considered concerning the reliability of the soft-sensor. In order to find out if the model with its detern)ined parameters is suitable at the moment when the soft-sensor is to be placed in operation, its performance during some period betbre the failure time must be assessed, h, the second case, when the soft-sensor is in operation, the reliability may be expressed in terms of the probability that the mean square error between the soft-sensor output al)d the now unmeasured actual variable does not exceed a given level. A simple way to obtain an estimate of the soft-sensor reliability is to determine the period during which the prediction error is less than a given threshold, 6% in this case. For doing this, data from the experimental campaign has been used. The number of sampling times when that error did not exceed a 6% (relative) divided by the total number of samples gives an estimate of the reliability of each model. As the soft-sensor should be based on the family of five models, the reliability of the whole system can also be estimated, considering that the soft-sensor system is available when at least one of the models is available. Table 3 contains the results obtained during the first 1500 samples period. 5. Conclusions

Good results have been obtained in the case of a soft-sensor for predicting the percentage over 65 mesh, designed for the

A. Casali et al / Powder Technology 99 (19981 15-21

Table 3 Soft-sensor reliability

Soft-sensor

Availability [~ l

Model A Model B Model C Model D Model E Soft-sensor system

63 29 61 49 45 83

Andina grinding circuit and developed on the basis of ARMAX type models which comprise on-line input measurements and physically significant combinations of such measurements. Thus, in order to determine the model structures, a stepwise regression method has been used but with a list of candidate components which includes combinations ot" measurements, derived from a phenomenological analysis. The models obtained in this way are much better than if only single measurements are used. A family of five soft-sensor models, all of them able to predict similarly the real measurement, was developed so that for the case of a missing input measurement a workable soft-sensor was available. During the first 550 rain, the average prediction error did not exceed 2% for all the models. The structures of the models are independent of the data set used in their construction if these sets correspond to a persistent excitation test. From the results, it can be observed that for the normal operating point there exists a set of parameters for each soft-sensor model that remains valid for almost all the testing period. In other words, there was no obsolescence of the parameters. An important factor from this result is the use of phenomenological components which causes an increased range of validity for a soft-sensor model for a given set of parameters [ 7 l. The plant seems to be for most of the time within a range of operating points that can be considered as normal operating conditions, where most of the models with fixed parameters work well. It is only during those periods that deviate significantly from the normal operating region where different parameter values are needed. It was found that a set of 180 samples (obtained in a persistent excitation test) is enough for the parameter identification of all the models. It is remarkable that one of the models (model B) is even capable of predicting correctly the starting points of the real sensor, coming from 0 to the normal operating point as shown in Fig. 6. By having a soft-sensor with a set of models to replace the percentage over 65 mesh, the availability of the measurement, including the soft-sensor substitute, is increased compared with the case when only one soft-sensor model is present. In fact, when the actual sensor i~ not available it may be replaced by one of the soft-sensor models (five in this case) depending on the availability of their input measurements. Even if one

21

of the input measurements of the best soft-sensor (i.e., using model ! ) fails during its operation, it may still be replaced by one of the alternative soft-sensors (using models 2, 3, 4 or 5 ). Thus the availability of the measurement over a given period, may be computed as the ratio between the sum of the period during which the sensor or some of the soft-sensor models is available, and the complete period.

Acknowledgements The research and development leading to this paper has been financed by CONICYT (National Science and Technology Council of Chile) through Project FONDEF MI-17. The authors express the thanks to CODELCO CHILEAndina Division for their valuable contribution during the experimeLttai campaigns in their grinding plant.

References [ ! l s. Kawatra, A. Seitz, An analysis of wet grinding circuit control using inferential particle size measurement. Preprints 82-18, SME-AIME Annual Meeting, Dallas, TX, USA, 1982. [ 2 ] G.D. Gonz,41ez,W. Meyer, Replacement of a particle size distribution sensor in grinding plant control by an estimator, Fourth Conference on Control Engineering Proceedings. The Inst. of Engineers. Gold Coast, Australia, 1990, pp. 65-69. [ 31 G.D. Gonz,41ez,M.A. Aguilera. L. Casteili. Developmentof a density soft-sensor for a mineral grinding plant, 12th IFAC World Congress, Vol. 5, Sydney, Australia. 1993, pp. 355-358. {4 ] A. Casali, G.D. Gonzfilez, F. Tortes, l. Cerda, L. Castelh, P. Gtm~nez, Pulp density soft-sensor for a grinding circuit, APCOM XXV, Brisbane, Australia, 1995, pp. 371-376. [ 51 M.T. Tham, Soft-sensors for process estimation and inferential control, J. Process Control I ( !99 ! ) 3- !4. 161 L. Ljung, System Identflication: Theory for the User, Prentice Hall, 1987. [ 71 G.D. Gonzfilez, A. Casali, R. Odgers, R. Barrera. F. Tortes, L. Castelh, P. Gim~nez. Soft-sensor design considering composite measurements and the effects of sampling periods, COPPER'95, Vol. It, Santiago, Chile, 1995, pp. 2 i 3-224. 181 D.A. Dahistrom, Fundamentals and applications of the liquid cyclone. Chem. Eng. Prog., Syrup, Series No. 15, Mineral Engineering Techniques, Vol. 50, 1954, p. 41. 191 A. Lynch, T. Rao, Modelling and scale-up of hydrocyclone classifiers, I I th Int. Mineral Processing Congress Proceedings, Cagliari. 1975. IlOl L.R. Plitt. J.A. Finch, B.C. Flintoff, Modeling the hydrocyclone classilier, European Symposium Particle Technology. European Fed. of Chemical Engineering, Dechema, Amsterdam, 1980. pp. 790--804. [ 11 ] J.A. Herbst, G.A. Grandy, T.S. Mika, On the developmentand use of lumped parameter models for continuous open and closed circuit grinding systems. Transactions, Inst. of Mining and Metallurgy, Section C. Vol. 80, 1971, pp. C!93-198. [12] J.A. Herbst. D.W. Fuerstenau. Scale-up procedure for continuou~ grinding mill design using population balance models. Int. J. Miner. Processing ( i 980) 1-3 I. [13] S.A. Billings. W.S.F. Voon, A prediction-error and stepwise-regression estimation algorithm for non-linear systems, Int. J. Control 44 ( 3 ) (1986) 803-822.