# Extremum-seeking control of retention for a microparticulate system.

INTRODUCTIONPulp and paper processes have undergone a significant evolution over the last 10 years. Numerous changes in paper processes such as the increase of paper machine speeds and the use of new raw materials, amongst others, have introduced large disturbances in the processes that can have a significant impact on paper quality and process operability. In this new context, the use of new control and optimization techniques have been proposed as possible mechanisms to improve the regulatory performance of paper processes.

In this study, the primary focus is the operation of the wet-end of the paper machine. The wet-end plays a vital role in the efficiency of the process and on the quality of the paper produced. The main challenge with the regulatory control of the wet-end of the paper machine is the relatively limited understanding of the physical and chemical phenomena that govern the sheet formation process and, hence, the absence of credible mathematical model descriptions. In addition, the process remains highly susceptible to the normal variations in the raw materials. The objective is therefore to develop tools that facilitate the design of control systems and optimization applications that can regulate the system despite the absence of reliable models and highly uncertain processing conditions. In this paper, we propose a real-time optimization (RTO) scheme to optimize the retention of fibres and fines in the head-box configuration of a paper machine. The RTO scheme reduces the need for complex and costly experiments that would be necessary for model development and provides an adaptation mechanism to changes in the operating conditions.

In this study a model of the wet-end section of the paper making process is proposed. The model is used to validate the design of a control scheme that can optimize the process with respect to key process parameters in real time. An empirical model of the retention as a function of fines, fibres, and filler concentration is considered. The control scheme is based on an adaptive extremum-seeking control technique that is used to steer the system to a region of the unknown maximum retention. The objective of extremum seeking is to bring a dynamical system to a setpoint that maximizes a relevant objective function. In the current context, it is employed as a real-time optimizing controller that can be coupled to any existing control techniques such as model predictive control and distribution shaping control. While conventional control ensures that the requested setpoints are tracked effectively, the extremum-seeking controller provides a means to compute, in real-time, setpoints that are optimal with respect to a user-specified cost function.

In the past few years, Krstic et al. (Krstic, 2000; Krstic and Wang, 2000; Wang et al., 1998) have presented several schemes for extremum-seeking control of nonlinear systems. Their framework allows the use of black-box objective functions, with the restriction that the objective value to be minimized is measured online. Although this technique has been proven useful for some applications, the lack of guaranteed transient performance of the black-box schemes remains a significant drawback in its application. The lack of transient performance of perturbation-based extremum-seeking control, such as the one presented here, remains on open area of research. Some recent work indicates that it may possible to enhance the transient performance of such techniques (Chioua et al., 2007).

In contrast, the extremum-seeking framework proposed by Guay and Zhang (Guay and Zhang, 2003) assumes that the objective function is explicitly known as a function of the system states and uncertain parameters from the system dynamic equations. Parametric uncertainties make the online reconstruction of the true cost impossible such that only an estimated value based on parameter estimates is available. The control objective is to simultaneously identify and regulate the system to the lowest cost operating point, which depends on the uncertain parameters. The main advantage of this approach is that one can guarantee some degree of transient performance while achieving the optimization objectives when a reasonable functional approximation of the objective function is available.

In cases where the objective function is available for measurement, one may consider the identification of a functional approximation of the objective expressed as a function of the decision variables using an adaptive estimation approach. This functional approximation is used in the context of the approach proposed by Guay and Zhang (2003) to achieve the extremum-seeking objectives. Although the black-box approach may prove to be a simpler alternative, the application of the approach of Guay and Zhang (2003) provides additional guarantees of transient performance and robustness of process uncertainties. A similar approach has been considered by Guay et al. (2005) for the optimization of nonisothermal CSTRs and by Guay et al. (2004) for the optimization of bioreactors.

The adaptive extremum-seeking technique proposed in this manuscript merges the model-based approach, as proposed by Guay and Zhang (2003), with the perturbation-based techniques (see Krstic, 2000; Krstic and Wang, 2000). The key of model-based approach is to make use of the process model and a functional form of the unmeasured cost function. On the other hand, perturbation based methods require the measurement of the cost function but do not necessitate a process model. In this paper, the cost function is assumed to be measured and the adaptive extremum-seeking mechanism is used to obtain a model of the cost function. Based on this estimated model, and following the adaptive-extremum-seeking control approach, the gradient of the cost function is estimated. The RTO objective is achieved by ensuring that the system is driven to a local optimum of cost function where the gradient vanishes. In comparison to the perbutation-based methods, the proposed method is not limited to local convergence results. It relies on the approximation of the cost function over a specific compact set of parameter values. The convergence of the parameters over this set ensures that the cost function is accurately approximated and hence the gradient. Convergence of the process to its optimum only relies on the property of the underlying dynamical system providing that sufficient excitation is injected to the system to allow estimation of the nonlinearity of the cost function. Thus, the result provides a non-local (semiglobal) convergence of the process to its unknown optimum. It generalizes existing results in perturbation-based methods such as Tan et al. (2006). The proposed technique generalizes the results presented by Guay et al. (2004) and Guay et al. (2005). The theoretical analysis provides a precise characterization of the performance of the extremum-seeking control system. The technique provides an approximation of the cost function and its gradient via the estimation of an output bias term. The results also demonstrate that the technique is directly applicable to the multivariate case.

The paper is organized as follows. In the next section, we present the model of the head-box system that is considered in this study. The model is based on a dynamic mass balance of the various micro-particulate species present in the retention process. An empirical model of the retention is then presented. The adaptive extremum-seeking control scheme is then presented. A simulation study is used to demonstrate the applicability of the technique.

HEAD-BOX SYSTEM DYNAMICAL MODEL

Let us consider the dynamics of the head-box system in the wet-end of a paper machine. The formation of a sheet of paper occurs as a result of filtration and thickening effects that depend mostly on the size of the particles (see Figure 1). The water solution from the suspension flows through the mesh and is recovered on the other side.

[FIGURE 1 OMITTED]

A portion of the fines (which are fragments of fibres), fibres, and filler is not retained and flows through to the water solution. The fraction of the particles retained in the sheet relative to the particles lost to the water solution is called the retention. retention o = quantity retained in paper

retention [??] quantity retained in paper / quantity introduced in the water suspension (1)

The filtration mechanism is primarily active for larger particles (>50 [micro]m). For the smaller particles (<10 [micro]m), this mechanism accounts for less than 5%. Therefore the small particles such as fines and filler are primarily retained by a thickening mechanism and mechanical entrapment. The water suspension gathered in the vicinity of the location where the flow of suspension enters the forming sheet is rich in fibres and other particles that are then recycled for the preparation of the water suspension.

To improve the retention, various polymer-based retention agents are added to the suspension. A high molecular weight neutral polymer (polyethylene oxide, PEO) is adsorbed to the fibres, fines, and filler particles. A cofactor is added that adsorbs to the polymer. In solution, the addition of the cofactor promotes flocculation of the particles that are therefore less likely to be lost to the broke.

Although retention cannot be measured directly, it can be easily inferred from the fraction of solids in the sheet and in the water suspension. Consistencies in the flow can be measured online by using optical methods. This measurement provides an inferential measurement of the retention. It is assumed that the optimization of the measured quantity leads to the optimization of the retention.

Let us now concentrate on the dynamical model of the head-box network. This requires both hydrodynamic considerations and kinetic information. A more detailed study of the adsorption mechanisms of the polymer and the cofactor with the solid particles in the water suspension is necessary to provide a complete model description of the physical phenomena.

The head-box network is represented schematically in Figure 2. In practice, the polymer and cofactor are added at different locations. In the model it is assumed that polymer, cofactor, and filler are added at the same location. It is further assumed that all reactions take place in the primary vessel. The following notation is used throughout this paper. Concentrations, denoted by [C.sub.i,j], are expressed in units of mass per unit volume. Indices i and j indicate the stream number, according to Figure 2, and the species of interest, respectively. Volumetric flowrates are denoted by [F.sub.i]. Mass flowrates are denoted by [M.sub.i] or [[??].sub.i]. The total volume of each unit is denoted by [V.sub.k].

The primary vessel and secondary vessel are assumed to be perfectly mixed. The sheet acts as a separation process whose dynamics are neglected. The two outlet flows are the water suspension flow ([F.sub.4]) and the flow of the forming paper sheet ([F.sub.3]). Since the concentration of the polymer and of the cofactor are typically very low, it is assumed that these components are retained in the paper.

The volumes of water in the primary and secondary vessels are assumed to be constant. There are therefore no dynamics associated with the volumes. The state variables of the system are the concentrations of fines, fibres, filler, polymer, and cofactor in each vessel.

[FIGURE 2 OMITTED]

The Primary Vessel

Let [V.sub.hb] be the volume of water in the primary vessel. Since [V.sub.h]b is constant, it follows that:

[F.sub.2] = [F.sub.mp] + [F.sub.pol] + [F.sub.1] + [F.sub.filler] + [F.sub.5] (2)

The primary vessel contains filler, fines, fibres, polymer, and cofactor. One can write the following material balances for each component:

d[C.sub.2,fibres] / dt = 1 / [V.sub.hb]( [M.sub.fibres] + [F.sub.5][C.sub.5,fibres] - [F.sub.2][C.sub.2,fibres]) (3)

d[C.sub.2,fibres] / dt = 1 / [V.sub.hb]( [M.sub.fines] + [F.sub.5][C.sub.5,fines] - [F.sub.2][C.sub.2,fines]) (4)

d[C.sub.2,filler] / dt = 1 / [V.sub.hb]( [M.sub.filler] + [F.sub.5][C.sub.5,filler] - [F.sub.2][C.sub.2,filler]) (5)

where [M.sub.i] is the total mass flow of species i where i = fibres, fines, filler, pol, mp. Consequently, we have

[M.sub.fibres] = [F.sub.1] [C.sub.1,fibres], [M.sub.fines] = [F.sub.1] [C.sub.1,fines], [M.sub.filler] = [F.sub.filler] [C.sub.filler] + [F.sub.1] [C.sub.1,fibres] [M.sub.pol] = [F.sub.pol] [C.sub.pol], [M.sub.mp] = [F.sub.mp] [C.sub.mp]

The quantities [C.sub.2,fibres], [C.sub.2,fines], and [C.sub.2,filler] represent the fibres, fines, and filler concentrations per unit volume of the flow [F.sub.2], respectively. Since the processes are assumed to be perfectly mixed, the concentrations described as above are also the concentrations inside the primary vessel.

For the polymer and cofactor, one has to take into account the rate of consumption by adsorption phenomena. Let [r.sub.pol] (*) be the rate of consumption of polymer in the primary vessel and let [C.sub.2,pol] be the concentration of polymer in the primary vessel. The material balance for the polymer can be then written as follows:

d[C.sub.2,pol] / dt = [M.sub.pol] - [F.sub.2] [C.sub.2,pol] - [r.sub.pol] (x)(6)

The rate of uptake of cofactor in the primary vessel is assumed to be proportional to the reaction rate of the polymer. Let [r.sub.mp](*) = [Kr.sub.pol](*) where K is some strictly positive constant (K > 0). The material balance for the cofactor is then given by

d[C.sub.2,mp] / dt = [M.sub.mp] - [F.sub.2] [C.sub.2,mp] - [Kr.sub.pol] (x) (7)

The secondary vessel

Since the presence of polymer and cofactor can be neglected in the broke, no chemical reaction takes place and the material balances reduce to the hydrodynamics of the process. Assuming constant volume, the overall material balance yields

[F.sub.4] = [F.sub.5] + [F.sub.6] (8)

The flowrate [F.sub.6] is a purge stream that is used to regulate the flowrate [F.sub.5] as a function of the inlet flowrates [F.sub.mp], [F.sub.pol], [F.sub.filler] and [F.sub.1] in order to keep [F.sub.2] constant. Consequently, we can assume that the consistencies in [F.sub.5] and [F.sub.6] are equal. Letting [V.sub.wp] be the water volume in the secondary vessel, the material balances for the fibres, fines, and filler is written as

d[C.sub.5,fibres]/dt = 1 [V.sub.wp] ([F.sub.4][C.sub.4,fibres] - [F.sub.6][C.sub.6,fibres] - [F.sub.5][C.sub.5,fibres]) (9)

d[C.sub.5,fines]/dt = 1/[V.sub.wp] ([F.sub.4][C.sub.4,fines] - [F.sub.6][C.sub.6,fines] - [F.sub.5][C.sub.5,fines]) (10)

d[C.sub.5,filler]/dt = 1/[V.sub.wp] ([F.sub.4][C.sub.4,filler] - [F.sub.6][C.sub.6,filler] - [F.sub.5][C.sub.5,filler]) (11)

Using the relation (8) and the hypothesis that the consistencies are the same in each outlet flows, the material balance equations reduce to

d[C.sub.5,fibres]/dt = [F.sub.4]/[V.sub.wp] ([C.sub.4,fibres] - [C.sub.5,fibres]) (12)

d[C.sub.5,fines]/dt = [F.sub.4]/[V.sub.wp] ([C.sub.4,fines] - [C.sub.5,fines]) (13)

d[C.sub.5,filler]/dt = [F.sub.4]/[V.sub.wp] ([C.sub.4,filler] - [C.sub.5,filler]) (14)

The flowrate [F.sub.5] is given as a function of all inlet flowrates and the flowrate [F.sub.2] following Equation (2).

The sheet of paper

As mentioned above, it is assumed that the dynamics of the forming paper sheet are negligible. The material balance equations for this unit reduce to a simple set of algebraic equations. These equations serve to link the inlet of the secondary vessel with the outlet of the primary vessel.

The sheet is considered as a separation unit that takes the flowrate [F.sub.2] and produces two outlet flows, [F.sub.3] and [F.sub.4] with different compositions. This unit is indeed modeled as a separator with an efficiency equal to the retention in the system. Applying the definition of retention given in Equation (1) for the fibres, fines, and filler, we can write

[R.sub.i] = [F.sub.3][C.sub.3,i]/[F.sub.2][C.sub.2,i] = 1 - [F.sub.4][C.sub.4,i]/[F.sub.2][C.sub.2,i], i = fibres, fines, filler (15)

The quantity of water remaining in the sheet is much smaller than the amount water that is collected in the secondary vessel. Consequently, we can write [F.sub.4] [much greater than] [F.sub.3] and, accordingly, it follows that [F.sub.4] [??] [F.sub.2]. Finally we state the material balance equations for the sheet as follows:

[C.sub.4,i] = [C.sub.2,i](1 - [R.sub.i]) i = fibres, fines, filler (16)

As mentioned above, the fibres are essentially retained by filtration. The fraction of fibres retained does not depend on the amount of particulate present in the primary vessel. As a result, we generally assume this fraction to be constant. The typical amount observed in practice is approximately 97%.

On the other hand, the fines and filler concentrations are primarily retained on the sheet following the addition of polymer. For these solid particles, the retention is therefore a function of the quantities entering the primary vessel. Since the mechanism of retention responsible for accumulation of fines and filler on the sheet are related to the size of the particles, it can be assumed that the retention of fines and filler in the sheet are the same: [R.sub.filler] = [R.sub.fines] = R.

Polymer Adsorption Kinetics

The mechanism of flocculation that results from the addition of polymer is extremely complex. The polymer molecule must first be adsorbed to the fines, fibres, and filler particles in the primary vessel. The polymer complex can then be adsorbed creating links between the various complexes. The cofactor adsorbs to the polymer complexes, promoting links between complexes.

This mechanism demonstrates that retention mechanisms are related to the adsorption of the polymer on the solid particles. In order to simplify the model, we do not directly consider phenomena that lead to the exchange or transfer of polymer molecules from one surface to the next. The model only considers adsorption of the polymer and of the cofactor from the solution to the surface of particles. Since the reaction involving the polymer is an adsorption to a surface, a Langmuir -Hinshelwood kinetic model is considered in line with many related studies (Asselman, 1999; Asselman and Garnier, 2000; van de Ven, 1994, 2005). The Langmuir-Hinshelwood kinetic model proposes that the polymer molecules do not interact with each other but only with a single free site on the surface. The adsorption is an equilibrium reaction given by

Polymer + [[GAMMA].sup.*] [??] [GAMMA]

where [[GAMMA].sup.*] represents a free active site on the surface and [GAMMA] is an occupied site.

The rate of adsorption is based on the mass-action law assuming that free active sites act as regular reactants. The number of occupied sites is proportional to the fraction of the surface covered by the polymer. If we let [[sigma].sub.i] be the surface of particulate i (= filler, fines, fibres) covered by the polymer per unit volume and if [S.sub.i] is the fraction of surface available for adsorption of the polymer per unit mass of particulate i, then the free surface in the primary vessel per unit volume is equal to [S.sub.i][C.sub.i] - [[sigma].sub.i]. Finally, the rate of adsorption of polymer on the particulate i is given by the expression:

[r.sub.pol[right arrow]i] = [k.sub.ads,i][C.sub.2,pol]([S.sub.i][C.sub.i] - [[sigma].sub.i]) - [k.sub.des,i][[sigma].sub.i] (17)

The notation [r.sub.pol[right arrow]i] refers to the rate of adsorption of the polymer on the particulate i. Consequently, it is necessary to introduce three new state variables for the primary vessel to account for the fraction of polymer that covers the active sites of the fines, fibres, and filler particles. The corresponding state equations are obtained by a material balance for each surface in the primary vessel. Assuming that none of the particle surfaces contain polymer in the inlet stream of the primary vessel, the following equations are obtained:

d[[sigma].sub.fibres]/dt = [k.sub.ads,fibres][C.sub.2,pol]([C.sub.2,fibres][S.sub.fibres] - [[sigma].sub.fibres])

-[k.sub.des,fibres][[sigma].sub.fibres] - [F.sub.2][[sigma].sub.fibres]/[V.sub.hb] (18)

d[[sigma].sub.fines]/dt = [k.sub.ads,fines][C.sub.2,pol]([C.sub.2,fines][S.sub.fines] - [[sigma].sub.fines])

-[k.sub.des,fines][[sigma].sub.fines] - [F.sub.2][[sigma].sub.fines]/[V.sub.hb] (19)

d[[sigma].sub.filler]/dt = [k.sub.ads,filler][C.sub.2,pol]([C.sub.2,filler][S.sub.filler] - [[sigma].sub.filler])

-[k.sub.des,filler][[sigma].sub.filler] - [F.sub.2][[sigma].sub.filler]/[V.sub.hb] (20)

The rate of consumption of the polymer is equal to the sum of the rates of adsorption to each particulate:

[r.sub.pol] = 1/[S.sub.pol] [[summation over (i)] ([k.sub.ads,i][C.sub.2,pol] ([C.sub.2,i][S.sub.i] - [[sigma].sub.i]) - [k.sub.des,i][[sigma].sub.i]) (21)

where [S.sub.pol] is the surface occupied by a unit mass of polymer.

Retention Model

The retention R is a function of the fraction of the surface of the particles that are covered by polymer molecules. Assuming that this fraction is at an equilibrium, the values of [[sigma].sub.fines] and [[sigma].sub.filler] become a function of [C.sub.2,fines] [C.sub.2,filler] and [C.sub.2,pol]. As a result we consider the following expression for the equilibrium value of the retention:

R = R ([C.sub.2,fines], [C.sub.2,filler], [C.sub.2,pol]) (22)

In general, there is no direct function available for the retention, but there exists some empirical expressions that have been obtained from empirical studies. In most experimental settings, the equilibrium retention value is expressed as a function of the mass inlet flowrates of each particles:

[bar.R] = [R.sub.emp] ([[bar.M].sub.filler], [[bar.M].sub.fines], [[bar.M].sub.pol], [[bar.M].sub.mp], [[bar.M].sub.fibres]) (23)

where [[bar.M].sub.i] represents the mass flowrate of species i. Following evidence gathered from an experimental study conducted by one of the authors, [R.sub.emp] is given by the following empirical relationship:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (24)

The variables [X.sub.i] are defined from the mass flowrates expressed in g/s as follows:

[X.sub.pulp] = [[bar.M].sub.fines] + [[bar.M].sub.fibres]/60, [X.sub.filler] = [[bar.M].sub.filler] 0.6([[bar.M].sub.fines] + [[bar.M].sub.fibres]),

Xpol = [[bar.M].sub.pol]/0.03, [X.sub.mp] = [[bar.M].sub.mp]/0.15

The parameters of the empirical function are given in Table 1. These parameters provide an adequate fit of the operating conditions of interest.

The dependence of the retention function on the quantities in the inlet flows to the primary vessel are based on the equilibrium conditions. In fact, one requires that, at the equilibrium, the following condition holds:

R ([[bar.C].sub.2,fines], [[bar.C].sub.2,filler], [[bar.C].sub.2,pol] = Remp ([[bar.M].sub.filler], [[bar.M].sub.fines], [[bar.M].sub.pol], [[bar.M].sub.[M.sub.mp]], [[bar.M].sub.fibres])

It can be easily shown that the equilibrium values of [[bar.C].sub.2,fines], [[bar.C].sub.2,filler], and [[bar.C].sub.2,pol] do not depend on [[bar.M].sub.mp] or [[bar.M].sub.fibres]. We assume that both [M.sub.mp] and [M.sub.fibres] are kept constant to some nominal values [M.sup.*.sub.mp] and [M.sup.*.sub.fibres] that are likely to hold during normal operation of the paper machine. The retention function that must be obtained is of the general form:

R ([[bar.C].sub.2,fines], [[bar.C].sub.2,filler], [[bar.C].sub.2,pol]) = [R'.sub.emp] ([[bar.M].sub.filler], [[bar.M].sub.fines], [[bar.M].sub.pol]) = Remp ([[bar.M].sub.filler], [[bar.M].sub.fines], [[bar.M].sub.pol], [[bar.M].sup.*.sub.mp], [[bar.M].sup.*.sub.fibres])

If one considers the equilibrium concentration of [C.sub.2,fines] and [C.sub.2,filler] as a function of the fixed inputs, it can be easily seen that the value of the equilibrium consistencies of the fines and filler particles do not require the knowledge of specific model constants. In fact, since an empirical relation is available, one can deduce the equilibrium value for the retention [bar.R] and compute the corresponding equilibrium values [[bar.C].sub.2,fines] and [[bar.C].sub.2,filler]. The approach considered here is to compute the effect of [M.sub.filler] and [M.sub.fines], at constant [M.sub.pol]. For each specific pair of values ([M.sub.filler], [M.sub.fines]), we compute the equilibrium retention using the empirical relation as well as the equilibrium consistencies of fines and filler particles. This approach provides a way to produce a plot of the retention as a function of [[bar.C].sub.2,fines] and of [[bar.C].sub.2,filler]. As shown in Figure 3, the dependence of retention on [[bar.C].sub.2,fines] and [[bar.C].sub.2,filler] can be well approximated by a plane. Note that this is only true for certain ranges of values of [M.sub.pol]. In this case, values in the range 0.55 [less than or equal to] [M.sub.pol] [less than or equal to] 0.65 g/min yield prediction errors less than 5%.

The procedure highlighted above was then repeated for various values of polymer flowrate. For each value, it was observed that the retention surface as a function of [[bar.C].sub.2,fines] and [[bar.C].sub.2,filler] was well approximated by a plane. The coefficients of the planar approximation

[bar.R] = [a.sub.0] + [a.sub.1] [[bar.C].sub.2,fines] + [a.sub.2] [[bar.C].sub.2,filler] (25)

were computed as a function of [M.sub.pol].

[FIGURE 3 OMITTED]

To check the validity of Equation (25), we substitute the equilibrium values of [C.sub.2,fines] and [C.sub.2,filler] in the retention expression to obtain:

[bar.R] = [a.sub.0] + [a.sub.1] [[bar.M].sub.fines]/[F.sub.2] - [[bar.F].sub.5](1 - [bar.R]) + [a.sub.2] [[bar.M].sub. filler]/[F.sub.2] - [[bar.F].sub.5](1 - [bar.R]) (26)

Solving this equation for [bar.R], we find the following relationship:

[F.sub.5] [[bar.R].sup.2] + R([F.sub.2] - [F.sub.5] - [a.sub.0][F.sub.5]) - [a.sub.0]([F.sub.2] - [F.sub.5]) - [a.sub.1] [[bar.M].sub.fines] - [a.sub.2] [[bar.M].sub.filler] = 0. (27)

This relationship can be compared to the empirical relation. This comparison is presented in the contour plot shown in Figure 4 for a given value of [M.sub.pol]. The dashed contours represent the contours of the actual function while the full lines are the contours obtained from the linear approximation expressed as a function of consistencies of fines and filler.

[FIGURE 4 OMITTED]

Reduced Order Model

It is well known that the absorption of the polymer on the fibres is negligible compared to that on the fines and filler. In the following, this will interpreted as, [[sigma].sub.fibres] [??] 0. (Strictly speaking, this assumption should be interpreted as [[sigma].sub.fibres] [[sigma].sub.fines] [??] 0 and [[sigma].sub.fibres] [[sigma].sub.filler] [??] 0). Moreover it is shown by Cho (2005) and van de Ven (2004) that the fibres consistency in the head-box has a low influence on the retention of the fines and of the filler as long as there are enough fibres for the paper sheet forming, and we shall therefore assume that the retention does not depend on the fibres consistency in the head-box. Finally the kinetics of the absorption of the polymer on the fines and the filler are very fast (reaction time betwen 0.1 and 1 s (Odberg et al., 1993; Asselman and Garnier, 2000)) and it can be therefore assumed that the surface of fines and of fibres covered by the polymer is at equilibrium

[[sigma].sub.i] = [k.sub.ads,I][C.sub.2,pol][C.sub.2,i][S.sub.i]/[k.sub.ads,i][C.sub.2,pol] + [k.sub.des,i] + [F.sub.2]/[V.sub.hb], i= fines, fibres (28)

We can therefore rewrite the dynamics of the process by considering the following set of six mass balance equations:

d[C.sub.2,fines]/dt = 1/[V.sub.hb] ([M.sub.fines] + [F.sub.5][C.sub.5,fines] - [F.sub.2][C.sub.2,fines]) (29)

d[C.sub.2,filler]/dt = 1/[V.sub.hb] ([M.sub.filler] + [F.sub.5][C.sub.5,filler] - [F.sub.2][C.sub.2,filler]) (30)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (31)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (32)

d[C.sub.5,fines]/dt = [F.sub.2]/[V.sub.wp] ([C.sub.2,fines](1 - R) - [C.sub.5,fines]) (33)

d[C.sub.5,filler]/dt = [F.sub.2]/[V.sub.wp] ([C.sub.2,filler](1 - R) - [C.sub.5,filler]) (34)

In the following section, we consider the development of an extremum-seeking control scheme. It is assumed that the input variables of interest are [M.sub.pol] and [M.sub.mp].

EXTREMUM-SEEKING CONTROL

General Problem and Assumptions

The process dynamics (22)(29)-(34) can be expressed as a nonlinear system of the form

[??] = f (x, u) (35)

y = h(x) (36)

For any fixed value of the input u, we assume that there exists an equilibrium value of y. The set of equilibrium values of y as a function of u is represented by R(u). That is, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] where y([t.sub.0], t; u, y([t.sub.0])) represents the value of y(t) at time t starting from initial condition y([t.sub.0]) at time [t.sub.0] under input u. Let us consider the following assumptions about the system dynamics.

Assumption 1: The output function R(u) is a convex function of u over a compact subset [S.sub.u] [subset] U and it admits of maximum over [S.sub.u].

Assumption 2: The dynamics (35) are input to state stable (ISS) with respect to u. Assumption 2 is easily checked for the dynamical system (29)-(34).

In this work, we choose to parameterize the unknown function R(u) via some universal approximation (e.g. a neural network) as

R(u(t)) = [theta] + [W.sup.T]S(u(t)) + [[mu].sub.l](t) (37)

with

S(u) = [[[s.sub.1](T(u)), [s.sub.2](T(u)), ..., [s.sub.l](T(u))].sup.T] (38)

One possible choice for the basis functions [s.sub.1](T(u)) is given as follows:

[s.sub.i](u) = 1/[square root of 2[pi][[sigma].sup.2]] exp [-(u - [[phi].sub.i])T (u - [[phi].sub.i])/[[sigma].sup.2.sub.i]] [(u - [[phi].sub.i]).sup.2]/[[sigma].sub.2]

for i=1,2, ..., l. The parameters of the basis function approximation (37) are assumed to take values in a known compact set [[OMEGA].sub.w]. The functional approximation of the equilibrium manifold is dependent on the centre of the receptive field, [[phi].sub.i], and the width of the Gaussian function [[sigma].sub.i]. Universal approximation results stated by Sanner and Slotine (1992) and Seshagiri and Khalil (2000) indicate that if l is chosen sufficiently large, then [W.sup.T]S(u) can approximate any continuous function to any desired accuracy on a compact set.

For the neural network approximation, we assume that the following holds.

Assumption 3: The neural network approximation errors satisfies [absolute value of [[mu].sub.l](t)] [less than or equal to] [[bar.[mu]].sub.l] with constant [[bar.[mu]].sub.l] > 0 over the compact set [[OMEGA].sub.w] x [S.sub.u]. We further assume that the approximation error is such that [absolute value of [partial derivative][[mu].sub.l](t)/[partial derivative]t] [less than or equal to] [bar.d[mu]] where [[bar.d[mu]] > 0 is positive constant.

Controller design in the univariate case

Let us first concentrate on the design of the adaptive extremum seeking when the input is the polymer flowrate:

u(t) = [M.sub.pol](t) (39)

The first step in the design of the controller is to provide a parameter estimation strategy. This step is crucial since it allows the reconstruction of the unknown nonlinearity R(u). Let [??] be the estimate of y and define the estimation error [??] = y - [??]. The time-derivative of y is given by

[??] = [partial derivative]R/[partial derivative]u [??] (40)

If one considers the basis function approximation Equation (37), then Equation (40) is written as follows:

[??] = [W.sup.T] [partial derivative]S/[partial derivative]u [??] + [partial derivative][[mu].sub.l]/[partial derivative]t (41)

Note that, by assumption, [partial derivative][[mu].sub.l]/[partial derivative]t exists and is bounded. This can always be ensured if the input is itself continuously differentiable and the approximation error [[mu].sub.l](t) is small.

Let [??] be the estimate of W. The estimate [??] is generated by the following estimation law:

[??] = [??] [partial derivative]S/[partial derivative]u [??] + c[(t).sup.T] [??] + kr(y - [??]) (42)

where c(t) [member of] [R.sup.l] is a vector-valued time varying function to be defined later, [??] is the parameter update law and [k.sub.r] > 0 is a positive constant. The error dynamics are given by:

[??] = [[??].sup.T] [partial derivative]S/[partial derivative]u [??] + [partial derivative][[mu].sub.l]/[partial derivative]t - c[(t).sup.T] [??] - [k.sub.r] [??] (43)

where [??] = W - [??] is the vector of parameter estimation errors.

The control objective is to maximize R(u) with respect to u. Since R(u) is convex and is known to exhibit a maximum over the compact set [S.sub.u], we only need to find the value of u, [u.sup.*], at which:

[partial derivative]R/[partial derivative]u ([u.sup.*]) = 0 (44)

Using the basis function approximation (37), Equation (44) is restated as follows:

[W.sup.T] [partial derivative]S/[partial derivative]u ([u.sub.a]) = 0 (45)

where [u.sub.a] is the approximation of the true solution [u.sup.*].

Since W is not known, we substitute its estimate [??] to define the quantity:

z = [[??].sup.T] [partial derivative]S/[partial derivative]u (46)

In the controller design, we consider a perturbed value of z given by:

[z.sub.s] = [[??].sup.T] [partial derivative]S/[partial derivative]u - d(t) (47)

where d(t) is a dither signal to be designed in the following.

As in our previous works, we adopt a Lyapunov based approach to design an extremum-seeking controller that ensures that u reaches a neighbourhood of [u.sup.*]. Let us define the filtered estimation error:

[eta] = [??] - c[(t).sup.T] [??] (48)

and the Lyapunov function candidate:

[V.sub.1] = 1/2 [[eta].sup.2] (49)

Its time derivative is given by the following expression:

[[??].sub.1] = [eta] ([[??].sup.T] [partial derivative]S/[partial derivative]u [??] + [partial derivative][[mu].sub.l]/ [partial derivative]t - [k.sub.r] [??] - [??][(t).sup.T] [??] (50)

Substituting [??] = [eta] + c[(t).sup.T] [??] , we have:

[[??].sub.1] = [eta] ([[??].sup.T] [partial derivative]S/[partial derivative]u [??] - [k.sub.r][eta] - [k.sub.r]c[(t).sup.T] [??] + [partial derivative][[mu].sub.l]/[partial derivative]t - [??][(t).sup.T] [??]) (51)

With this formulation, we can compute [??] to remove the dependency of [[??].sub.1] on [??] by considering the following dynamical equation:

[??] = -[k.sub.1]c(t) + [??] [partial derivative]S/[partial derivative]u (52)

where [k.sub.1] is a positive gain assigned by the user. As a result, we get

[[??].sub.1] = -[k.sub.1][[eta].sup.2] + [partial derivative][[mu].sub.l]/[partial derivative]t [eta] (53)

Using the upper bound on [absolute value of [partial derivative][[mu].sub.l]/[partial derivative]t] = [bar.d[mu]], we finally obtain

[[??].sub.1] [less than or equal to] -[k.sub.1] [absolute value of [eta]] ([absolute value of [eta]] - [bar.[partial derivative][mu]/[k.sub.1] (54)

Consequently, it follows from standard Lyapunov stability theory arguments that [eta] is bounded and enters the region

{[eta] | [absolute value of [eta]] < [bar.d[mu]]/[square root of [k.sub.1]]} as t [right arrow] [infinity].

By completing the squares, inequality (54) can be rewritten as follows:

[[??].sub.1] [less than or equal to] 1/2[k.sub.1] [[bar.d[mu]].sup.2] - ([k.sub.r - [k.sub.1]/2) [[eta].sup.2] (55)

for some positive constant [k.sub.1]. If we choose [k.sub.1] < 2[k.sub.r] and define the [[lambda].sub.1] = [k.sub.r] - [k.sub.1]/2, we have

[[??].sub.1] [less than or equal to] 1/2[k.sub.1] [[[bar.d[mu]].sup.2] - 2[[lambda].sub.1][V.sub.1] (56)

and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (57)

The optimization objective is expressed in terms of the Lyapunov function:

[V.sub.2] = 1/2 [z.sup.2.sub.s] (58)

Its time derivative is given by

[[??].sub.2] = [z.sub.s] ([[??].sup.T][[partial derivative].sup.2]S/[partial derivative][u.sup.2] [??] + [[??].sup.T] [partial derivative]S/[partial derivative]u - [??] (59)

The existence of a maximum for the approximation of R(u) imposes the constraint [[??].sup.T] [[partial derivative].sup.2]S/[partial derivative][u.sup.2] < 0 over [[OMEGA].sub.w] x [S.sub.u]. Thus, given that [[OMEGA].sub.w] is chosen accordingly, we can consider the dynamical control law:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (60)

Hence we get that:

[[??].sub.2] = -[k.sub.z][z.sup.2.sub.s] (61)

from which we deduce that [z.sub.s] [right arrow] 0 as t [right arrow] [infinity]. The convergence of [z.sub.S] and [eta] to zero guarantees that all signals of the closed-loop system are bounded. However, it does not imply that one recovers the optimum. This is clear since the convergence of [eta] to zero alone cannot guarantee that either [??] nor [??] converges to zero. A proper choice of the parameter update law guarantees that the optimization objective can be achieved. Such an update law is given as follows:

[??] = [[gamma].sub.w]c(t)[??] (62)

where [[gamma].sub.w] is strictly positive constant. By assumption, the parameters W belong to a convex set described by the known elementwise intervals given by [W.sub.i] [member of] [[a.sub.i], [b.sub.i]] where i=1,2..., l. To ensure that the parameter estimates, [[??].sub.i], remain within the intervals [[a.sub.i], [b.sub.i]], i=1, ...,l, we consider the following projection algorithm:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (63)

This algorithm guarantees that the parameter estimates [[??].sub.i] remain within the intervals [[a.sub.i] - [delta], [b.sub.i] + [delta]]. As a result, the tuning parameter [delta] must be a small, strictly positive constant. Using Equations (37) and (62), the parameter estimation error dynamics are given by

[??] = - [??] = -[[gamma].sub.c](t)(y - [??]) = -[[gamma].sub.c](t) ([eta] + [c.sup.T](t) [??]) = -[[gamma].sub.c](t)[eta] - [[gamma].sub.c](t)[c.sup.T](t)[??] (64)

Let us now consider the following lemma.

Lemma 1: Consider the dynamical system:

[??](t) = -[[phi].sup.T](t)[phi](t)z(t)

and assume that there exist strictly positive constants T and k such that

[[integral].sup.t+T.sub.t] [[phi].sup.T] ([tau])[phi]([tau])d[tau] [greater than or equal to] kI

Then the origin is exponential stable equilibrium of the system. This lemma and its proof can be found by Anderson et al. (1986), Theorem 2.2. As a result, it follows that [??](t) can be shown to converge asymptotically to the origin if the dynamics of c(t) fulfills the assumption of Lemma 1. The vector c(t) must fulfill the persistency of excitation condition which states that there exist strictly positive constants T and k such that:

[[integral].sup.t+T.sub.t] c([tau])[c.sup.T]([tau])d[tau] [greater than or equal to] kI (65)

i.e. that the dither signal d(t) must be chosen to provide sufficient excitation. We now consider the following Lyapunov function candidate:

[V.sub.3](t) = 1/2[[gamma].sub.w] [[??].sup.T] [??] (66)

Its time derivative is given by

[[??].sub.3] = -[[??].sup.T]c(t)[c.sup.T](t) [??] - [[??].sup.T]c(t)[eta] (67)

We complete the squares and introduce a strictly positive constant [k.sup.*.sub.2] ([k.sup.*.sub.2] > 0) to obtain the following inequality:

[[??].sub.3] [less than or equal to] 1/2[k.sup.*.sub.2] [[eta].sup.2] - (1 - [k.sup.*.sub.2]/2) [[??].sup.T] c(t)[c.sup.T] (t) [??] (68)

The constant [k.sup.*.sub.2] must be chosen such that [k.sup.*.sub.2] < 2. In the remainder, we use the constant [k.sub.2] rather than [k.sup.*.sub.2] with [k.sub.2]/2 = (1 - [k.sup.*.sub.2]/2]) > 0. From Lemma 1, we obtain the following persistency of excitation condition:

[[??].sup.T]c(t)[c.sup.T](t) [??] [greater than or equal to] [[alpha].sub.1] [[parallel][??][parallel].sup.2] (69)

where [[alpha].sub.1] is a strictly positive constant. Consequently it follows that:

[[??].sub.3] [less than or equal to] -[[alpha].sub.1][k.sub.2]/2 [[parallel][??][parallel].sup.2] + 1/2[k.sup.*.sub.2] [[eta].sup.2] = -[[gamma].sub.w][[alpha].sub.1][k.sub.2][V.sub.3] + 1/2[k.sup.*.sub.2] [[eta].sup.2] (70)

Upon substitution of Equation (57), we obtain the following inequality:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (71)

By integrating Equation (71), the following inequality is obtained:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (72)

or

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (73)

The result confirms that [??] convergences exponentially to a small neighbourhood of the origin which we state as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (74)

where [[alpha].sub.1] and [[beta].sub.3] are strictly positive constants.

The last step in the design is to provide a parameter estimation algorithm for [theta] that was introduced in the original Equation (37). It is worth noting that this variable does not affect the optimization routine. However, its knowledge provides a better approximation of the equilibrium surface for y as a function of u. The adaptation law is simply chosen as follows:

[??] = -[k.sub.[theta]] ([[??].sup.T]S(u) + [??] - y) (75)

The parameter estimation error dynamics for [theta] are given by

[??] = [??] = [k.sub.[theta]] ([[??].sup.T]S(u) + [??] -[W.sup.T]S(u) - [theta] - [[mu].sub.l](t)) = [k.sub.[theta]] [??]S(u) - [k.sub.[theta]] [??] - [k.sub.[theta]] [[mu].sub.l](t) (76)

where [??] = [theta] - [??]. We consider the following candidate Lyapunov function:

[V.sub.4](t) = 1/2[k.sub.[theta]] [[??].sup.2] (77)

Its time derivative is equal to

[[??].sub.4] = [??] (-[??] - [[??].sup.T]S(u) - [[mu].sub.l](t))

leading to

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (78)

where [k.sub.3] and [k.sub.4] are strictly positive constants that are introduced in the completion of squares. These constants follow the relation 1 > [k.sub.3]/2 + [k.sub.4]/2 . In what follows, we set

[[lambda].sub.3] = [k.sub.[theta]] (1 - [k.sub.3]/2 - [k.sub.4]/2) (79)

From Equation (78), we deduce

[[??].sub.4] [less than or equal to] - [[lambda].sub.3]/[k.sub.[theta]] [[??].sup.2] + [[mu].sub.l][(t).sup.2]/2[k.sub.3] + [(p[??].sup.T]S(u)).sup.2]/2[k.sub.4] (80)

We note that the function [[mu].sub.l](t) is bounded. Similarly, the signal u (t) is bounded and so is S (u). Let [bar.S] represent the bound on [parallel] S [parallel]. We can then write the following:

[[??].sub.4] [less than or equal to] -[[lambda].sub.3]/[k.sub.[theta]] [[??].sup.2] + [[bar.[mu]].sup.2]/2[k.sub.3] + 1/2[k.sub.4] [[parallel][??][parallel].sup.2] [[bar.S].sup.2] (81)

Substituting Equation (74), we obtain:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Integration of the last inequality yields the following expression:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (82)

where

* [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

* [[beta].sub.5] = 2[k.sub.[theta]] max [[[beta].sub.4], [absolute value of [[beta].sub.3] [[bar.S].sub.2]/2[k.sub.4](2[[lambda].sub.3] - 2[[lambda].sub.2])]]

* [[lambda].sub.4] = min[[[lambda].sub.3], [[lambda].sub.2]]

The last inequality proves that [??] approaches a small neighbourhood of the origin asymptotically. Based on the analysis, it is therefore possible to reconstruct the function R (u) in a neighbourhood of the optimal one. We state the result as follows.

Theorem 1: Consider the nonlinear system (35) (36) subject to assumptions (1), (2), and (3). If there exists a time T and a constant k such that the persistency of excitation condition (65) is satisfied, then the extremum-seeking controller (60), (42), (52), and (63) steers the control system to a neighbourhood of the extremum of the measured objective function y.

Multivariate Optimization

Since the polymer is combined with a cofactor in order to improve the retention, it seems interesting to consider the control design when both flowrates of polymer and cofactor are considered as control inputs, i.e. [u.sub.1](t)=[M.sub.pol](t) and [u.sub.2](t)=[M.sub.mp](t). In this section, we consider the case where the retention is treated as a function of [u.sub.1](t) and [u.sub.2](t) assumed to independent. The input u(t) is represented in vector form as u(t) = [[[u.sub.1](t) [u.sub.2](t)].sup.T[. The retention is given as the unknown function:

y = R([M.sub.pol],[M.sub.mp]) = R([u.sub.1], [u.sub.2]) = R(u) (83)

This function is approximated using a basis function approximation. The basis function vector is represented by S (u) which is defined in two-dimensional space as follows:

R(u(t)) = [W.sup.T]S(u(t)) + [theta] + [[mu].sub.l](t) (84)

with weights W [member of] [R.sup.l], [theta] [member of] R being a constant. We consider a set of Gaussian basis functions given by:

[s.sub.i](u) = 1/[square root of 2[pi][[sigma].sup.2]] exp -[(u - [[phi].sub.i]).sup.T] (u - [[phi].sub.i])/[[sigma].sup.2]) - [(u - [[phi].sub.i]).sup.T](u - [[phi].sub.i])/[[sigma].sup.2]

The centres of the Gaussian functions, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], are evenly located over a rectangular grid in the space of input variables [u.sub.1] and [u.sub.2]. It is assumed that the approximation error [[mu].sub.l] (t) and its derivative are bounded.

For the control design we first consider the dynamics of the retention expression. It is given by

[??] = [partial derivative]R/[partial derivative]u [??] = [W.sup.T[ [partial derivative]S/[partial derivative]u [??] + [partial derivative][[mu].sub.l]/[partial derivative]t (85)

where [partial derivative]S/[partial derivative]u is the 2 x 2 Jacobian matrix of S (u). We pose the estimation of y as [??] which is obtained for the following expression equation:

[??] = [[??].sup.T] [partial derivative]S/[partial derivative]u [??] + c(t) [??] + [k.sub.r](y - [??]) (86)

where c(t) [member of] [R.sup.l] and [k.sub.r] is a strictly positive constant. The goal of the control law is to find the input value u (t) that maximizes the retention y. The optimal vector value [u.sup.*] of the input is the vector which is such that:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (87)

The optimum of the approximate surface is obtained at the input vector value [u.sub.a] such that:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (88)

(a) Convergence of [??]: Let us define the variable [eta] as follows:

[eta] = [??] - [c.sup.T](t) [??] (89)

where [??] = y - [??], [??] = W - [??]. We consider the Lyapunov function candidate:

[V.sub.1](t) = 1/2 [[eta].sup.2] (90)

Its time derivative is given by

[[??].sub.1](t) = [eta] ([[??].sup.T] (dS/du [??] - [k.sub.r]c(t) - [??](t)) + d[[mu].sub.l]/dt - [k.sub.r][eta] (91)

As a result we are led to assign the vector function c (t) as the solution of the differential equation:

[??](t) = dS/du [??] - [k.sub.r]c(t) (92)

The derivative of [V.sub.1] (t) becomes

[[??].sub.1] = -[k.sub.r][[eta].sup.2] + d[[mu].sub.l]/dt [eta] [less than or equal to] -[k.sub.r] [absolute value of [eta]] ([absolute value of [eta]] - [bar.d[[mu].sub.l]/[k.sub.r])

From LaSalle's invariance principle, we conclude that the function [eta](t) is bounded and converges exponentially to the set: {[eta] | ([absolute value of [eta]] < [bar.d[[mu].sub.l]/[k.sub.r]}

(b) Convergence to the optimum: Let us define the vector [z.sub.s] [member of] [R.sup.2] as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (93)

At the optimum of the approximate surface, we have:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Consequently we choose the candidate Lyapunov function:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (94)

whose time derivative is equal to

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (95)

This last expression allows one to choose the control law such that the optimum can be reached. In fact, if the control input is chosen as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (96)

Then Equation (95) becomes

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (97)

which proves that z[s.sub.1] and z[s.sub.2] tend to zero asymptotically and the input u (t) converges to the optimal value corresponding to the optimum of the approximate surface.

(c) Convergence of the weight parameters [??] : The adaptation law for the weights [??] is given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (98)

The update law is based on a projection algorithm which guarantees that parameter estimate [[??].sub.i] remains in the interval [[a.sub.i], [b.sub.i]]. The dynamics of the parameter estimation error [??] is then given by

[??] = -[??] = -[[gamma].sub.w]c(t)(y - [??]) = -[[gamma].sub.w]c(t) ([eta] + [c.sup.T](t) [??]) = -[[gamma].sub.w]c(t)[eta] - [[gamma].sub.w]c(t)[c.sup.T](t) [??] (99)

To guarantee the convergence we must ensure that the signal c(t) obeys to a persistency of excitation condition that requires that there exists a time [T.sub.0] and a strictly positive constant [beta] such that:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (100)

where I is the l by l identity matrix. The parameter [theta] is obtained by the adaptation law:

[??] = -[k.sub.[theta]] ([[??].sup.T]S(u) + [theta] - y) (101)

Using a standard adaptive control argument, one can easily show that the parameter estimates [??] and [??] converge to W and [theta]. As a result the corresponding input converges to the optimal value associated with the optimum value of the approximate surface and the goal of the control law is achieved.

SIMULATION RESULTS

In this section, the control law is tested by using the nominal dynamical model presented above. The empirical retention function is used both for the one variable case with polymer flowrate as the control input, and for the two variable case using the polymer and cofactor flowrates as the input variables. The excitation signal d (t) is chosen as follows:

[??] = -[k.sub.d]d(t) + a(t) (102)

where the signal a (t) ensuring persistent excitation is given by:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (103)

The centres of the Gaussian basis functions are chosen to be regularly spaced on the interval [[u.sub.min], [u.sub.max]] in the univariate case, and on the grid 4 x 4 in the subset of [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] in the multivariate case. The frequencies are chosen to span a range in the low and high frequency range. Typically, this choice should reflect the dynamics of the system and the total number of parameters estimated. In this case, the choice of 10 different frequencies was deemed sufficient to ensure sufficient excitation. The choice of the amplitudes, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] dictate the extent of the perturbations used for the dither signal. In (Adetola and Guay, 2007), it is shown how one can assign amplitudes of the signals to guarantee satisfaction of the persistency of excitation condition while minimizing the magnitude of the dither signal.

The values of the controller tuning parameters used in the simulation are the following ones in the one variable simulations:

[k.sub.[theta]] = 30, [k.sub.d] = 80, [k.sub.r] = 10, [k.sub.z] = 20, [[gamma].sub.w] = 20000

while the following values have been considered in the two variable simulations:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The results of the first case study are shown in Figures 5 and 6. The piecewise constant red line highlights the position of the unknown optimum steady-state conditions. In order to illustrate the parameter convergence, one can use the model parameters to reconstruct the retention function. The comparison is shown in Figure 6.

The adaptive extremum-seeking controller performances in the two variable case are illustrated in Figures 7 to 11. The first two figures illustrate the convergence of the retention to the optimum. In order to illustrate that the function approximation provides a good estimate of the retention in a neighborhood of the optimum, the contour lines of the actual function and its approximation [??] = [??]S(u) + [??] are given in Figure 9. The figure shows that the contour levels of the approximate function (dotted lines) approximates the actual contour levels of the actual retention function (full lines).

In the both case studies, we have investigated the performance of the control law subject to operation condition changes. In this case, we assume that an abrupt change occurs at t=12 min. The results shown in Figures 10 and 11 confirm the effectiveness of the control law to maintain the dynamical system at the maximum equilibrium retention.

[FIGURE 5 OMITTED]

[FIGURE 6 OMITTED]

[FIGURE 7 OMITTED]

[FIGURE 8 OMITTED]

[FIGURE 9 OMITTED]

[FIGURE 10 OMITTED]

[FIGURE 11 OMITTED]

CONCLUSIONS

In this study, we have demonstrated the applicability of the extremum-seeking controller for the maximization of the retention in the wet-end of a paper machine.

First, we have attempted to provide a simple dynamic model of the hydrodynamics and the kinetics of the retention of particles that occur in the wet-end as a sheet of paper is formed. Despite the fact that the exact mechanisms are somewhat poorly, we consider a case where an empirical expression relating the retention and polymer and cofactor mass flowrates in the primary vessel. We also show one could provide interpolation methods that can generate locally validated models of the retention in terms of the system state variables.

Using the model developed, we considered the design of RTO control system that can be used to maximize the retention despite the absence of suitable models. The control law is based on a uniquely steady-state optimization of the system which neglects all dynamics. It was shown that the control law can effectively steer the dynamical control system to a neighbourhood of its optimal value. In this study, the standard extremum-seeking algorithm for single variable problems is generalized to the multivariate case where it is shown that convergence to the optimum can be achieved when a suitable persistency of excitation condition is met. Simulation results have confirmed that this control is valid for typical applications.

ACKNOWLEDGEMENTS

This paper presents research results of the Belgian Network DYSCO (Dynamical Systems, Control, and Optimization), funded by the Interuniversity Attraction Poles Programme, initiated by the Belgian State, Science Policy Office. The scientific responsibility rests with its authors.

END NOTES

For solids, consistency is the mass per unit volume of water

Manuscript received January 5, 2008; revised manuscript received June 9, 2008; accepted for publication June 21, 2008.

REFERENCES

Adetola, V. and M. Guay, "Parameter Convergence in Adaptive Extremum-Seeking Control," Automatica 43(1), 105-110 (2007).

Anderson, B., R. Bitmead, C. Johnson, P. Kokotovic, R. Kosut, I. Mareels, L. Praly and B. Ridle, "Stability of Adaptive Systems--Passivity and Averaging Analysis," MIT press, Cambridge (1986).

Asselman, T., "Heterofloculation of Wood Fibres and Fines Induced by Polymers and Microparticules," PhD Thesis, McGill University, Montreal (1999).

Asselman, T. and G. Garnier, "Dynamics of Polymer-Induced Hetero-Floculation of Wood Fibres and Fines," Colloids surf. A Physicochem. Eng. Asp. 147, 297-306 (2000).

Chioua, M., B. Srinivasan, M. Guay and M. Perrier, "Dependence of the Error in the Optimal Solution of Perturbation-Based Extremum Seeking Methods on the Excitation Frequency," Can. J. Chem. Eng. 85(4), 447-453 (2007).

Cho, B., "Control Strategies for Retention and Formation on a Paper Machine using a Microparticulate Retention Aid System," PhD Thesis, McGill University, Montreal (2005).

Guay, M. and T. Zhang, "Adaptive Extremum Seeking Control of Nonlinear Dynamic Systems with Parametric Uncertainties," Automatica 39, 1283-1293 (2003).

Guay, M., D. Dochain and M. Perrier, "Adaptive Extremum Seeking Control of Continuous Stirred Tank Bioreactors with Unknown Growth Kinetics," Automatica 40, 881-888 (2004).

Guay, M., D. Dochain and M. Perrier, "Adaptive Extremum Seeking Control of Nonisothermal," CSTR. Chem. Eng. Sci. 60(13), 3671-3681 (2005).

Krstic, M., "Performance Improvement and Limitation in Extremum Seeking Control," Syst. Control Lett. 39(5), 313-326 (2000).

Krstic, M. and H. Wang, "Stability of Extremum Seeking Feedback for General Dynamic Systems," Automatica 36(4), 595-601 (2000).

Odberg, L., H. Tanaka and A. Swerin, "Kinetic Aspects of the Adsorption of Polymers on Cellulosic Fibres," Nord. Pulp Pap. Pres. J. 1, 159-166 (1993).

Sanner, R. M. and J. J. E. Slotine, "Gaussian Networks for Direct Adaptive Control," IEEE Trans. Neural Netw. 3(6), 837-863 (1992).

Seshagiri, S. and H. K. Khalil, "Output Feedback Control of Nonlinear Systems Using RBF Neural Networks," IEEE Trans. Neural Netw. 11(1), 69-79 (2000).

Tan, Y., D. Nesic and I. Mareels, "On Non-Local Stability Properties of Extremum-Seeking Control," Automatica 42(6), 889-903 (2006).

van de Ven, T., "Kinetic aspects of polymer and polyelectrolyte adsorption on surfaces," Adv. Colloid Interface Sci. 48, 121-140 (1994).

van de Ven, T., M. Qasaimeh and J. Paris, "Peo-Induced Flocculation of Fines: Effects of Peo Dissolution Conditions and Shear History," Colloids surface A Physicochem. Eng. Asp. 248, 151-156 (2004).

van de Ven, T., "Association-Induced Polymer Bridging by Poly(ethylene oxide)-Cofactor Flocculation Systems," Adv. Colloid Interface Sci. 114-115, 147-157 (2005).

Wang, H., S. Yeung and M. Krstic, "Experimental Application of Extremum Seeking on an Axial Flow Compressor," Proceedings of the American Control Conference, Philadelphia (1998), pp. 1989-1993.

Audrey Favache, (1) Denis Dochain, (1) * Michel Perrier (2) and Martin Guay (3)

(1.) CESAME, Universite Catholique de Louvain, Louvain-la-Neuve, Belgium

(2.) Departement de Genie Chimique, Ecole Polytechnique de Montreal, Montreal, PQ, Canada

(3.) Department of Chemical Engineering, Queen's University, Kingston, Ontario, Canada K7L 3N6, Canada

* Author to whom correspondence may be addressed. E-mail address: denis.dochain@uclouvain.be

Table 1. Parameter values of the empirical retention model Parameter Numerical value [a.sub.pulp] 2.72 [a.sub.filler] 27.47 [a.sub.pol] 30.91 [b.sub.pol] -1.92 [b.sub.mp] -0.4044 [a.sub.mp] 0.8940 [b.sub.pulp,filler] 0 [b.sub.pulp,pol] -6.36 [b.sub.pulp,mp] -2.66 [b.sub.filler,pol] -81.75 [b.sub.filler,mp] 0 [b.sub.pol,mp] 2.09 c -9.65

Printer friendly Cite/link Email Feedback | |

Author: | Favache, Audrey; Dochain, Denis; Perrier, Michel; Guay, Martin |
---|---|

Publication: | Canadian Journal of Chemical Engineering |

Date: | Oct 1, 2008 |

Words: | 10018 |

Previous Article: | Parameter and state estimation in nonlinear stochastic continuous-time dynamic models with unknown disturbance intensity. |

Next Article: | Preamble for special issue honouring John F. MacGregor. |