Paper Machine Modeling at Norske Skog Saugbrugs : A Mechanistic Approach

In this paper a mechanistic model of a paper machine is presented. The model is developed as a foundation for the control of three selected variables, the basis weight, the paper ash content and the white water total concentration. The model is of high order and reduced order models are developed and fitted to experimental mill data. The fitted models are validated with historical operational data.


Introduction
A paper machine is a complex process due to its multivariable nature and mixture of physical, chemical and mechanical sub-processes.Several researchers consider modeling of this process to be an impossible task (see e.g.(Roberts 1996b, page 8)), and no denying: an all-including model would not be possible given the present knowledge.The approach taken in this paper is one in which we focus on a mechanistic model 1 which will be used in an MPC (Model Predictive Control) application.There are three (manipulated) inputs and three outputs in the model.Several more inputs are present in the model and these will be considered as "measured disturbances".The model is simplified to make it more suitable for control purposes.It is beyond the scope of this paper to present a model which in all aspects have the correct physical structure, however it is important that the model captures the essential dynamic Paper C: Paper Machine Modeling at Norske Skog Saugbrugs: ... 123 behavior of the process and that it is applicable over a wider range of operating conditions than would be expected from an empiric model2 .The manipulated inputs and outputs of the model are as shown in Figure 1.The modeling is part of a larger project for stabilization of the wet end of paper machine 6 (PM6) at Norske Skog Saugbrugs, Norway (Hauge, Ergon, Forsland, Slora & Lie 2000).Norske Skog Saugbrugs is the world's second largest manufacturer of uncoated super calendered magazine paper (Norske Skog 2000), and the mill incorporates three paper machines.PM6 is built in the early 1990's and is the largest and most modern papermachine at the Saugbrugs mill.Empirical modeling or system identification of paper machines are reported in several papers and books.Some of these focus on so-called cross-directional (CD) modeling (i.e. a model for the profile across the paper web), e.g.(Featherstone, VanAntwerp & Braatz 2000), (Campbell 1997) and (Heaven, Manness, Vu & Vyse 1996), while others focus on the machine-direction (i.e.changes in average values across the web), e.g.(Hauge et al. 2000), (Hauge & Lie 2000), (Menani, Koivo, Huhtelin & Kuusisto 1998) and (Noreus & Saltin 1998).
The reported works on mechanistic modeling of paper machines are in most cases constrained to smaller parts of the paper machine.However, (Rao, Xia & Ying 1994), (Larsson & Olsson 1996) and (Hagberg & Isaksson 1993) consider a larger part of the paper machine, e.g. the wet end and the wire, press, and dryer sections, although the chemistry involved in papermaking is not considered at all.As far as we know, the only mechanistic models of a larger part of a paper machine which includes chemical modelling is found in (Shirt 1997), (Hauge et al. 2000) and (Hauge & Lie 2000).In Shirt's work both chemical aspects, which include adsorption and flocculation, and physical aspects, which include drainage on the wire, refining, tanks, headbox, wire section, etc., are part of the overall model, although transportation delays in pipelines are neglected and not all aspects are presented in detail.The mechanistic model in (Hauge et al. 2000) and (Hauge & Lie 2000) is in those papers not presented in detail, giving only an introduction to the equations used to describe the paper machine.In  neither of these latter references are the mechanistic models validated properly with real time data.The contributions in this paper are to bring more details on the model, to report on further refinements of the model, and to bring results from the model fitting and validation using experimental data.
A simplified overview of PM6 is given in Figure 2. Note that the developed models cover the process from the thick stock pump to the reel.A thick stock (lower left area of Figure 2, including the mixing-and machine-chests) model is developed in (Slora 2001), and this model is implemented in PM6 and provides estimates of totaland filler concentrations in the thick stock.
The paper is organized as follows: In Section 2 the paper machine sub-processes are discussed in detail and suggestions for equations describing them are given.In Section 3 we describe how simplified mechanistic models, more suitable for control purposes, may be obtained.In Section 4 we report on fitting and validation of the simplified models using real time process data.We improve the simplified mechanistic model by extending it with a first order empiric model in Section 5.In this section we also identify a Kalman filter and validate the model.Finally some conclusions are given in Section 6.

A Comprehensive Mechanistic Approach
The model described in this section consists of 28 ordinary di erential equations (ODE), 100 partial di erential equations (PDE) and hundreds of algebraic equations.For implementation in Matlab, each PDE is discretized by the method of lines (MOL) (see e.g.(Schiesser 1991)).

Chests and tanks
Chests and tanks are modelled as ideally stirred volumes, i.e. where is the density of the mass flow and q j [m 3 = s] is the volumetric flow j into the volume.We will get back to the "ideally stirred" assumption when discussing the deculator and the white water tank, and the generation term when discussing the retention aid and flocculation.
where v is the velocity of the mass inside the pipeline and x is the variable corresponding to the direction along the pipeline.Thus, the concentration is here a function of both time t and space x.When the reaction rates are small such that the advection term dominates, these models are notoriously di cult to discretize using the method of lines (MOL).With constant velocity v and r i = 0, these models can be interpreted as time delays.For implementation, the partial di erential equations are discretized into five ordinary di erential equations each, i.e. a pipeline is approximated by dividing it into five ideally stirred volumes.The choice of discretization level is a trade-o between factors such as accuracy, complexity and numerical properties.With an increasing number of volumes, the model is more accurate but also more complex and the sti ness of the overall system is increased.Keeping in mind that the model is developed for control purposes, we found that the input-output properties of the overall model, with a discretization level of five volumes, was very close to higher levels of discretization.The choice of five volumes was also convenient as a starting point with respect to model complexity.The trade-o can be studied from the responses in Figure 3, where a step change (from 0 to 0:1) in the initial concentration is applied to the pipeline between the screens and the headbox.The pipeline is 40 m long, it has a cross section area of 0:7 m 2 and a mass flow of 2500 kg = s.A density of ½ = 1000 kg = m 3 is assumed.If the pipeline is a pure time delay then the step change would appear at the outlet at t = L=(w=(A½)) = 11:2 s, where L is the length of the pipeline.

Fibers and fillers
The finished paper typically consists of 65% (wood-) fibers, 30% filler particles and 5% water.The filler particles are added to improve certain properties of the paper, such as brightness and smoothness, and also often to reduce the production costs.At PM6 two di erent types of filler particles are used depending on the requirements from the end-user.One of these type of filler particles is used occasionally when particularly high brightness is required.
The fibers that enter the process come in a variety of dimensions, and may be crudely classified as fibers and fines where strict definitions of fines appear in the literature (Britt & Unbehend 1976).The fines are generated in the refining process by the shearing action of the refiner bars upon the fiber cell walls (Roberts 1996a).

The thick stock
The thick stock area is the lower left area of Figure 2. Cellulose, thermomechanical pulp (TMP), broke (repulped fibers and fillers from sheet breaks and edge trimmings), recovered stock from the disk filters (thickened mixture of cellulose, white water, and more) and filler particles are the main additives, and they are blended in the mixingand machine chests.These tanks have relatively large volumes for smoothing rapid changes in the additive flows.The stock is transported to the paper machine area (Figure 2, except for lower left area) by the thick stock pump.A mechanistic model of the thick stock area, estimating the total and filler consistencies in the flow to the paper machine, is developed and implemented at PM6 (Slora 2001).
The components in the flow from the thick stock pump are assumed to be fibers, fines and filler particles.The total-3 and filler-concentrations are estimated by the thick stock model, and the concentration of fines is assumed to be where ® fines is a constant and C total and C filler are total-and filler-concentrations in the thick stock.

Retention aid and flocculation
The filler particles and fines are in general too small to be trapped on the wire (a fine meshed woven cloth) and one adds retention aid to help them flocculate, thus increasing the possibility of mechanical entrapment.The retention aid is also added for other reasons, such as improving the drainage from the sheet, but these e ects will not be studied here; see e.g.(Roberts 1996b) for a general introduction to retention aids and flocculation.
Fibers and filler particles are mostly negatively charged, and at PM6 a two component cationic (positively charged) retention aid system is used.The two components have quite di erent charge densities and molecular weights.A low molecular weight, high charge density polymer is added first, mainly to fix or neutralize highly anionic (negatively charged) impurities in the stock but also to improve retention as illustrated in Figure 4. We will assume that the flocculation due to this polymer is negligible as is also experienced on the paper machine when the addition of this polymer is stopped.The addition of this polymer is now used in a control loop for stabilizing the charge (or cationic demand) of the stock at PM6.
The second polymer that is added has a low charge density and high molecular weight.This results in adsorption onto fibers and filler particles leaving "tails" which other fibers and filler particles may adsorb onto (see Figure 5).This is termed "bridging flocculation" and is assumed to be the main contributor to the flocculation in the process.The addition of this retention aid is located just after the screens as seen on Figure 2. On a more microscopic level one may go on to describe the adsorption and flocculation of various components (e.g.adsorption of polymers onto fibers, fines and filler particles, and the flocculation of fibers with fibers, fibers with fines, filler particles with filler particles, etc.) as is found in e.g.(Van de Ven 1993) and (Shirt 1997).This, however requires many (ordinary and partial) di erential equations, and a simpler solution was sought here.An equation for the concentration of flocculated component i that provided an overall good fit for the model, and which also was relatively simple was @C floc;i This equation describes concentration of flocculated components in pipelines (following the convention of Equation 3), with obvious extensions to equations for chests and tanks as outlined in Equation 2. Component i is either filler particles or fines (we are not interested in flocculation of fibres as we assume later on that they will be retained on the wire in any case), C floc;i is the concentration of flocculated component i, C non-floc;i is the concentration of non-flocculated component i, k i;ret and k i;fiber are flocculation constants and C ret and C fiber are concentrations of retention aid and fibers, respectively.The reaction rate is empiric and should not be regarded as an attempt to explain the complicated process of flocculation.This is obvious considering that as long as there are fibers, the equation predicts a positive reaction rate even with no filler particles present in the system.
In addition to the flocculation equation ( 5) itself, there are several other issues which must be addressed.For example that the equation fails to account for several major disturbances on the flocculation.The temperature and charge are two of these disturbances and they are stabilized at PM6 by separate control loops.Another parameter which has tremendous e ect on the flocculation is the pH (see e.g.(Horn & Linhart 1996)), and at PM6 the pH is observed for longer periods of time showing no dramatic variations.
The retention aid is added at the screens, and only seconds later it arrives at the wire section where obviously some of it must be recirculated through the wire cloth.We will get back to this issue in Section 2.11.The recirculated polymer will undergo a deactivation process in which it di uses into the pores of the fibres (Koethe & Scott 1993) thus loosing some ability to cause flocculation (Pelssers, Cohen Stuart & Fleer 1989).In addition, the recirculated flocs must pass through a number of pumps and other process equipment (such as the hydrocyclones and the screens) with high shear rates, causing break-up of flocculated particles4 (Gregory 1988).These facts caused (Shirt 1997), in his dynamic model, to assume that all particles which are recirculated through the white water system loose any active high molecular weight (low charge density) polymer coverage.In our model we allow for some flocculation to take place in the white water system, thus giving some initial concentrations (larger than zero) of flocs at the outlet of the screens where the retention aid is added.

The white water tank
The white water tank is modeled as a perfectly stirred tank (Equation 2).The validity of this assumption may be questioned as the main input flow (the dewatering from the wire section) enters on top of the tank and the main output flow (the flow to the first stage of the hydrocyclones) leaves at the bottom of the tank, and there are no mixing arrangement present in this tank.The tank is however only moderately tall and the main input and output flows are quite large and have high velocities, so the turbulence inside the tank is assumed to blend the contents well.Another approach is taken in (Nissinen 1999) where the upper part of the tank is assumed perfectly mixed, and the lower part a plug flow.

The hydrocyclones
There are a total of five hydrocyclone stages and a tank between the fourth and fifth stage, as shown in Figure 6.Each stage consists of several units where the accept flow goes through the top of the cyclones and the reject flows out at the bottom.
In the model we have neglected the time delays, or volumes, of the hydrocyclone units themselves, because these volumes are small compared to the pipeline volumes between the units.The pipelines are modeled as shown in Equation 3, and the cyclone units split the incoming flows into an accept flow and a reject flow.The concentrations in the accept flow will be lower than in the reject flow because gravitational forces are used to separate the two outlet flows and heavier particles will have a tendency to be rejected.
Equations like are used for finding the accepted total mass flows in stage i ( i {1; :::; 5}), and are used for finding the concentration of component j in the accept flow at stage i.Component j is recirculated filler particles, "fresh" filler particles, fiber, fines, retention aid, flocculated fines or flocculated filler.Observations after step changes in the addition of filler particles indicate that recirculated and newly added filler behaves quite di erently in the hydrocyclones, and this is the reason for treating them separately.
Due to high shear rates and deactivation (as explained in Section 2.5), it is assumed that no flocculation takes place in the fourth and fifth stage, and also that the polymer entering the fourth stage is completely inactive.The tank between the fourth and fifth stage is modeled by ordinary di erential equations similar to Equation 2, but without the flocculation term, and only for the components filler, fiber and fines.

The deculator
The deculator is a relatively small two chamber tank whose purpose is to remove air bubbles from the stock.The "right side" chamber (refer to Figure 2) has the largest volume and an overflow to the "left side" chamber keeps the level constant.
The level on the "left" side is controlled and assumed constant.Equations similar to Equation 2 are used for both chambers.The assumption of ideally stirred masses in each chamber is probably good due to large masses entering and the small volumes involved.In addition, the "right" side chamber has an input flow (recirculated from the headbox) entering at the bottom which should make the assumption even better.

The screens
There are two screens in parallel, in addition to a reject system for recirculation of usable fibers and filler particles back into the white water tank.In the model the screens are seen as one unit, splitting the flow in a reject and an accept flow.The equations used are similar to those used for the hydrocyclones, thus ignoring the volume of the screens which are small.The reject system is not part of the model due to the reject flow and concentrations being small.The reject flow is about 2% of the inject flow and the concentrations are found by laboratory measurements to be no larger than the inject and accept concentrations.

The headbox
The headbox is modeled as a pure mass flow splitting unit.The two output flows, one recirculation flow to the deculator and one flow onto the wire, are assumed to have the same concentrations as the input flow.The mass flow to the wire section is calculated by which can be derived from Bernoulli's equation, where w H W [kg = s] is the mass flow onto the wire, ® I [kg 1=2 = m 5=2 ] is a lumped constant (dependent on the geometry of the slice opening amongst others), h I and `I are the height and length of the slice opening respectively and P I [Pa] is the pressure inside the headbox.Note that the length of the slice opening is a constant while the height of the slice opening is used for controlling the cross directional (CD) profile.Thus the height varies across the slice opening and the average is used in the equation above.More detailed models of headboxes can be found in e.g.(Rao et al. 1994) and (Tuladhar, Davies, Yim & Woods 1997).

The wire, press, and dryer sections
Following (Shirt 1997) we assume that (long) fibers, flocculated fines and flocculated filler particles are retained on the wire 5 .In addition we allow for some filler particles to be mechanically entrapped, so that the component mass flows entrapped on the wire are where w W;i [kg = s] are the component mass flows retained on the wire (i {filler, fines, fiber}), C H;j are the concentrations in the headbox (j {flocculated filler, non-flocculated filler, flocculated fines, fiber}), and ® W;filler is the fraction of nonflocculated filler mechanically entrapped on the wire.The significance of mechanical entrapment seems to be somewhat controversial in the literature.For example (Van de Ven 1984) found (theoretically) that mechanical entrapment was low, while (Bown 1996) reports that mechanical entrapment can be a dominant mechanism.
Most of the water is drained from the wire to the white water tank, and this is modeled by where w W W W T [kg = s] and w H W [kg = s] are total mass flows from the wire to the white water tank and from the headbox to the wire respectively, while ® W is the fraction (close to one) of mass flow out on the wire which is recirculated to the white water tank.
The concentrations of filler particles, fines and fibers in the mass flow from the wire section to the white water tank are easily calculated using Equations 8-12 and the concentrations in the headbox (the concentration of fibers in this flow is equal to zero).As explained in Section 2.5 we also allow for some recirculation of retention aid, and the concentration of retention aid in the flow from the wire to the white water tank is calculated as where C W W W T;ret.aid and C H;ret.aid are concentrations of retention aid in the flow from the wire to the white water tank and in the headbox respectively, and ® ret.aid is the fraction of retention aid being recirculated.The parameter ® ret.aid may be viewed as a lumped parameter considering that we do not account for deactivation and destruction of polymers by shear rates elsewhere in the model.An exception is in the hydrocyclones, where we assume that polymer entering the fourth stage is completely inactive.The transportation of the solids from the wire section, through the pressing and dryer section, and onto the reel is modeled by an advection equation where w W R;i [kg = s] is the mass flow of component i (i =filler or solids (filler, fiber and fines added)) from the wire section to the reel, and v [m = s] is the paper machine velocity (near the reel).We do not model the filtration process or drainage process in any detail, and we only focus on the solids on the wire.

The output equations
The outputs are as shown in Figure 1, and the equations connecting the outputs to the internal states of the dynamic model are where y basis weight (t) [g = m 2 ], y paper ash (t) [%], and y W W tot. conc.(t) [%] are the basis weight, paper ash content, and white water total concentration respectively.The basis weight and paper ash are measured by a scanning device between the dryer section and the reel.® edge trim is a constant which adjusts for edge trimmings of the paper, is the mass flow of component i (i = filler and/or solids (filler, fiber and fines added)) at the scanning device (L paper indicates the length of the paper sheet from the wire section to the scanning device), v(t) [m = s] is the paper machine velocity at the scanning device, b R [m] is the width of the paper sheet at the scanning device and f(t) is the measured moisture content in the paper at the scanning device.C W W W T;filler (t) and C W W W T;fines (t) are the concentrations of filler and fines in the flow from the wire section to the white water tank.
3 Model reduction

Finite dimensional models
Assume that we have an input-output model where x R n , u R r and y R m .The dimension of the model is n, and we want to reduce the order of the model.We illustrate the principle using a linear model: Assume that matrix A has a full set of eigenvectors such that there exist a non-singular matrix M and a diagonal matrix such that AM = M .We can then introduce a change of variables and let z = M 1 x such that we find Thus, we have a decoupled system: . . .
The system is stable if i : <¸i 0. For simplicity, assume that all eigenvalues ¸i are real, and that we have ordered the eigenvalues such that ¸1 ¸2 • • • ¸n 0. If now for some k we have ¸k ¿ ¸k+1 , it is common to introduce the approximation dz i =dt 0 for i {1; : : : ; k}.Thus we have reduced the order of the model to n k, and it is now: . . .
In addition, we have the algebraic equations A more formal analysis of this approximation can be performed using singular perturbation techniques, see e.g.(Hoppensteadt 2000).The technique of singular perturbation is also applicable to nonlinear systems such as dx=dt = f (x; u).However, for nonlinear systems, it is much more di cult to find a transformation of the state such that the transformed system is decoupled.Such a transformation would ideally have the form: where z R n , and Á is invertible yielding x = Á 1 (z; u).In that case, we would search for a transformation Á such that where f (z; u) has the following structure: Assuming that such a transformation is available, the method of singular perturbation is readily applicable to the transformed system.
In practice, we can not expect to find such a decoupling transformation for nonlinear systems.A partial solution could be to linearize the model at an operating point, and then let Á (x; u) = M 1 x where M is the matrix of non-singular eigenvectors at this operating point.This would yield an approximate decoupling in the neighborhood of the operating point.
A simple case where we in fact have the desired decoupling, is in cascaded processes.This is relatively common, but feedback will then easily ruin the decoupling.
Even if there is no such natural decoupling, we may have an almost decoupled system.We may thus still try to use singular perturbation techniques on a system and thus set dx i =dt 0, even if this approximation is not strictly correct.
The final word regarding what eigenvalue/mode of the system we can allow ourselves to set at steady state, will depend on the use of the model.If the model is used for control purposes, we will have to select the model reduction such that the input-output response is not severely changed.If we are interested in the accuracy of some state, however, we will have to make sure that this particular state is not severely changed.

Infinite dimensional models
Assume that we have an infinite dimensional model (i.e., the model includes PDEs) which we discretize using the Method of Lines (MOL) such that we arrive at a model where dim x depends on the discretization level N .This case can be treated just as in the previous subsection.The particular problem with infinite dimensional models is that as part of the model reduction, we have to decide the discretization level.If there are several PDEs in the model, these can in principle have individual discretization levels.
When using the model for control, a natural way to solve the problem of finding the appropriate discretization level is to vary N such that the input-output behavior is not severely changed from the infinitely (high) dimensional case.If we are particularly concerned with the approximating capabilities of an internal state, we should select the discretization level by making N as small as possible while retaining the accuracy of the state in question.

A simplified paper machine model
In (Hauge & Lie 2000) it was shown, by simulation, that the full scale model described in Section 2 can be reduced to a much simpler and lower order model without a ecting its properties to any large extent.The implemented full scale model is of order 528, and it was shown that a 38th order model have basically the same input-output properties.The simplifications were by and large carried out by lumping pipeline volumes into existing tanks, discretizing the remaining pipeline volume (the pipeline between the screens and the headbox) into one volume and discretizing the wire, press, and dryer sections into one "volume".The response from the full scale model and the simplified model, when filtered PRBSs (Pseudo Random Binary Signal) are applied to the inputs, is shown in Figure 7.
Further model reductions were obtained by e.g.comparing the responses of simplified models of various order, fitted to process data.The simplifications were carried out with a step by step approach in which the model was carefully studied after each phase of simplifications and model fitting.The model reductions culminated in a third order model which will be presented next.
In Figure 8, only elements relevant for the simplified model are shown.There are three lumped volumes, and these are the white water tank, the reject tank and the deculator ("right" side).Only two components are accounted for in the simplified model, and these are filler particles and fibers, thus no flocculated filler or fines concentrations are calculated throughout the white water system.We will go briefly  The thick stock The fiber and filler concentrations are estimated.In addition we calculate the concentration of "fresh" filler in the thick stock flow C T S;fresh filler = ® T S;fresh filler • C T S;filler . ( Here, only a share ® T S;fresh filler of the filler is assumed "fresh", while all the filler added at the outlet of the white water tank is assumed "fresh".The reason for this is that the filler from the thick stock is assumed better mixed when entering the hydrocyclone arrangement. The white water tank The white water tank is modeled as a perfectly stirred tank with one flow entering (from the wire) and one flow leaving (to the hydrocyclones).
The two components in the tank are filler particles and fiber, and the filler particles are modeled by an ODE (ordinary di erential equation), while the fibers are at steady state.
The hydrocyclones The hydrocyclones consist of two stages and a reject tank between the two stages, as shown in Figure 8.The equations used to describe the hydrocyclones are equal to those presented in section 2.7, although fines, retention aid, flocculated fines, and flocculated filler particles are not part of the model.In addition the ODE for fiber is at steady state, thus only the ODE for filler particles in the tank remain.
The deculator The ODEs in the "left" side are at steady state, and only the ODE for fiber in the "right" side remain.
The headbox The screens are neglected and the flow from the "right" side of the deculator goes directly to the headbox.Between the deculator and the headbox the retention aid is added, and the equation for flocculation of filler particles, as described in section 2.5, is at steady state.Since fines are no longer a component in the model, no equations for flocculated fines exist.The flow which is recirculated from the headbox to the deculator is removed from the model.
The wire We assume that a part ® W;fiber of the fibers are retained on the wire where we follow the notational convention of section 2.11.Furthermore, we assume that flocculated filler particles are retained on the wire, in addition to some mechanical entrapment of non-flocculated filler particles The output equations The output equations are equal to those presented in section 2.12, except that the transportation delay through the wire, press, and dryer sections are neglected, and no fines are accounted for

Summary of simplified third order model
The model equations are where w T S is the flow of thick stock, w filler is the flow of filler which is added at the outlet of the white water tank, w ret.aid is the flow of retention aid added at the outlet of the screens, and the outputs are as explained in section 2.12.
The measured disturbances, which are accounted for in the model, are where C T S;total and C T S;filler are total and filler thick stock concentrations, S 1.stage pump is the speed of the pump between the white water tank and the first stage of the hydrocyclones, v is the machine velocity, P H is the pressure inside the headbox, h slice opening is the height of the slice opening and f is the paper moisture percentage.
The parameter vector µ consists of various more or less unknown parameters, which we tune to fit the model to process data.Several other parameters exist, some which are known, and some which are set at fixed values due to identifiability considerations.The parameter vector then, is where ® i;j;accept is the share of accepted component j (fresh filler, filler or fiber) at the i'th hydrocyclone stage, ® W;fiber and ® W;filler are shares of fiber and filler mechanically entrapped on the wire, ® T S;fresh filler is the share of fresh filler in the filler flow from the thick stock, V i are the volumes of the white water tank, reject tank and deculator ("right" side), and k j are flocculation constants for filler particles.

Implementation issues
The algebraic equations in the model are calculated in an order similar to the physical appearance of the variables in the process (e.g.algebraic equations associated with the first stage of the hydrocyclones are computed before equations associated with the reject tank and the second stage).The advantage of calculating them in this order is that the model file is well arranged, and changes in the model can be easily implemented and tested.The disadvantage is that due to several equations being mutually dependent, we can not compute all of them as they appear in the model equations.In the model file this is solved using shifted values for some variables, i.e.
The method used here corresponds to the fixed-point iteration method (see e.g.(Gerald & Wheatley 1994)), using only one iteration.A simple explicit Euler method is used for discretizing the model in Equation 35.One reason for using the explicit Euler method, contrary to e.g. a Runge-Kutta method, is that it can be used in a straight forward manner, even though we use a fixed-point iteration method with one iteration to calculate some of the algebraic equations.
It is possible, by substitution, to eliminate the algebraic equations from the model.We have compared validation results (as in Section 4), using a model where the algebraic equations are eliminated, and a model where the fixed-point iteration method using one iteration is used.The two methods gave practically the same validation results and the model outputs were close to indistinguishable.The reason for this is that most of the variables in the process change slowly and thus the error of using a delayed value is small.The model file, when the algebraic equations are eliminated, consists of a few very large and complex equations.Changing the structure of the model would not be possible using this file, and in the remainder of this paper we use the model file where fixed-point iteration is used.

Criterion and minimization algorithm
The least squares criterion is used for measuring the model fit where e is a vector of errors and Q is a diagonal weighting matrix.The minimization algorithm's single task is to solve where μ is the estimated parameter vector.
The function lsqnonlin in the Optimization toolbox (version 2) in Matlab (version 6) is used for solving the minimization problem.The function relies on the Levenberg-Marquardt algorithm in its search for the optimal parameter values (The MathWorks, Inc. 2000).
There are at least two alternatives when deciding how the errors e should be calculated.In the prediction error method (PEM) and in subspace methods one calculates the prediction error where y(t) is the measured output at time t, and ŷ(t|t 1) is the predicted output at time t based on past input-output data, i.e. a one-step-ahead prediction.In this case the error vector would be e.g. with where m is the number of outputs and N is the number of samples in the data set.
Another approach is to simulate the system, with only the initial conditions given.The error is then where ŷ(t|0) is the model output at time t given only the initial conditions.The error vector for output i is then

7)
Traditional system identification is carried out by using the one-step-ahead method, however in our case we wish to emphasize the need for a model with good long term prediction abilities.The reason for this is that the model will be used for model predictive control (MPC).Then, it seems natural to use the simulation approach in the parameter estimation algorithm.The simulation approach results in a deterministic model, and it will be necessary to identify or model the noise as well.We will return to this issue soon.

Experiment design and preprocessing of data
For a linear system, the concept of persistent excitation (see e.g.(Ljung 1999) and (Söderström & Stoica 1989)) provides an adequate characterization of the input signal needed to identify the model.The order of the excitation is dependent on the power spectrum of the signal only, and is independent of its shape (e.g.amplitude).This is not the case for a non-linear system; because such systems are amplitude dependent, i.e. the response to an input sequence u(t) may be qualitatively very di erent from that for a • u(t).One has the opportunity to design optimal experiments (Goodwin & Payne 1977) where ˆ is the optimal experimental condition, e.g.shape of inputs and sampling times, z is the set of feasible values for , and J( ) is a scalar criterion.However, this approach is seldom practicable for mechanistic models due to a chain of assumptions which must be made, e.g.choosing nominal parameter values.Often, experiment design for non-linear systems is based on a few rules of thumb, e.g. to use an input sequence of several amplitude levels (Pearson & Ogunnaike 1997), and to excite the relevant frequency bands.On a paper machine an open loop experiment is carried out with high risk of poor paper quality or even sheet breaks.A solution to this problem is to perform closed loop experiments, i.e. in this case experiments where the basis weight, paper ash and wire tray (or white water) total consistency controllers are in automatic mode.There is a vast amount of published material on closed loop system identification, and various approaches and algorithms are treated in more detail in e.g.(Ljung 1999), (Söderström & Stoica 1989) and (Forssell 1999).Our approach is "the direct approach" (Ljung 1999) in which we use the process inputs u and outputs y in the same way as for open loop identification, ignoring the feedback mechanisms and the reference signals.We specify changes in the setpoints, thus forcing the inputs to perturb the process.For example a rough approximation of the filtered PRBS signal is possible by changing the setpoints of the mass flows according to a PRBS scheme and let the valve and pump controllers work out the correct openings and velocities.Such an experiment plan is shown in Figure 9.There is no need to introduce several amplitude levels in the plan, since the process inputs and outputs are far from typical binary signals.The resulting inputs and outputs are shown in Figures 10 and 11 respectively.Filtered data are used when the deterministic model was identified, while the rawdata are used when the stochastic part of the model was identified.The filtering was carried out by a second order Butterworth filter, and as is seen in Figure 11, a cubic interpolation routine was applied to the paper ash data in a region near the 125th minute due to erroneous measurements.

Model fitting and validation of deterministic model
For comparison, we compute the value of a root mean square error (RMSE) criterion where N is the number of observations, y i (t) is the measured value of output i at time t, and ŷi (t) is the predicted or simulated value of output i at time t.
Optimal fitting of the mechanistic model to experimental data was carried out as described in section 4.1.The fitted model output is shown in Figure 12, along with     the measured output.Note that only a deterministic model, or ballistic model, is used in this simulation.The basis weight and paper ash dynamics seem well captured by the model, while the wire tray concentration fit is less good.In addition to scaling of the parameters in Equation 39, the outputs in the criterion (Equation 41) where also scaled such that more weight was put on the wire tray concentration.Adding even more weight on the wire tray concentration did not improve the model fit for this particular output.
RMSE values for the fitted third order mechanistic model is shown below: Basis weight Paper ash Wire tray cons.0:30 0:33 0:022 For validation of the mechanistic model, three data sets containing operational data were used.The first data set was collected during March 8-11, 2001.Figure 13 shows the validation results, when the mechanistic model is simulated (ballistic) with the measured inputs.Although corrected for bias, the model fit for basis weight and paper ash seem reasonably good considering that the time span is more than 90 hours.The initial oscillations in the basis weight may be caused by large oscillations in the estimated thick stock consistencies at this time.
RMSE values for validation of the third order mechanistic model with operational data collected during March 8-11, 2001, are:  Basis weight Paper ash Wire tray cons.1:17 0:66 0:037

A Hybrid Extension
A hybrid model will be introduced in this section.The hybrid model consists of the third order mechanistic model and a first order empiric model.The empiric model is identified with the error signals in Equation 47 as outputs, u and various measured disturbance signals as inputs.Denote the output from the mechanistic model by ŷmec and the output from the empiric model by ŷemp , then the hybrid model output is ŷhyb = ŷmec ŷemp .The empiric model may be viewed as originating from neglected and unknown dynamics in the third order mechanistic model.The neglected dynamics typically arise from the model reductions carried out, but also from such sources as sensors and actuators.Unknown dynamics can e.g.be a filter in a measuring device with proprietary software, some physical unit not known to the person modelling the system or a general lack of understanding of the physical process.We will start this section by fitting and validating a deterministic hybrid model.The validation will be carried out with operational data sets spanning several days.It is of course not realistic or the intention to use the models to predict several days ahead, so we will then identify a stochastic sub-model, and validate the combined deterministic and stochastic model.This validation will be carried out close to realistic conditions, validating the one-step ahead prediction abilities and the long-term prediction abilities.

Model fitting and validation of deterministic hybrid model
We use the experimental data set from February 28th, 2001, to identify the empiric model.The criterion and functions used are similar to those used with the mechanistic model (see Section 4.1).The model structure is The two sources of measured disturbances are the estimated thick stock concentrations.Although it would be preferable to add other disturbance sources, this is not possible using the February 28th, 2001, dataset due to lack of excitation from other sources.An alternative could be to use operational data spanning several days, such that more measured disturbances were excited.We will return to this issue soon, but point out that using operational data to identify e.g. the time constant of the process, which in this case is a simple transformation of A, in this case fails because the process itself is not properly excited.With operational data and a first order system, we found A 1, which is an integrator.Validation of the model identified from operational data was not successful.The basis weight and paper ash validation results are better than for the pure mechanistic model, although this is not the case for the wire tray consistency where the result is poorer.
We mentioned earlier that only two measured disturbances were used in the empiric model, due to lack of excitation from other sources.This is true for the short    experimentation data set used to identify the models, but in the operational data sets which span several days, many measured disturbances vary quite a lot.We take advantage of this and identify the E and F matrices in the state space model of Equation 50 a new, and in addition augmenting the d vector such that d R 6 .The other system matrices are not altered.The measured disturbance vector used in the empiric model is equal to that of Equation 38, except that the last element (f ) is not part of the empiric model.Figure 16

Identification and validation of "quasi extended" Kalman filter
Identification An extended Kalman filter normally updates the Kalman filter gain matrix at each sample, based on noise covariance matrices and a linearization of the model (see e.g.(Gelb 1974) and (Ergon 2001)).For simplicity we skip the linearization of the model and the identification or tuning of the covariance matrices, and identify the Kalman filter gain matrix directly from data as shown in Figure 19.Thus, the Kalman filter is not updated and therefore the name "quasi extended".Note that we do not use filtered data when we identify the Kalman filter.problem is to "freeze" the corresponding inputs at their present values (the values at the time of the sheet break).At PM6, the basis weight and paper ash measurements are lost during sheet breaks, while we still measure the wire tray concentration.Thus, in Figure 22 and 23 we have validated the model, simulating sheet breaks of various lengths.The sheet breaks last from 30 minutes to 2 hours, and they take place during normal operation or during grade changes.The simulation is carried out such that at sheet breaks we use the identified constant gain Kalman filter matrix, and the innovation signal for basis weight and paper ash is "frozen" at the mean value of their ten values prior to the sheet break.The innovation signal for wire tray concentration is calculated at each sample as before.The paper ash seems to be predicted reasonably good during sheet breaks, while the result for the basis weight is more mixed considering the drift in the 40-41st hour in Figure 22.An alternative, and perhaps more common, method for updating the states with missing measurements, is to set the innovation signal to zero for the lost measurements.The states will then only be updated through the available measurement.

One-step ahead validation
Validation of prediction ability during grade changes Prior to grade changes, the operators must give information to the control system about the time of change, new basis weight, paper ash, and other variables.With this information the long term prediction ability of the model must be reasonably good.grade changes.The prediction horizon ranges from 2 to 3 hours.The innovation signal for all three model outputs is "frozen" at the mean value of their ten values prior to the long-term prediction.

Conclusions
In this paper we have presented a high order mechanistic model of paper machine 6 (PM6) at Norske Skog Saugbrugs in Halden, Norway.The model is simplified making it more suitable for control purposes, and a third order mechanistic model is outlined.The third order model is fitted to experimental data and validated with historical operational data.In the experimental data set only a few of the measured disturbances are properly excited, thus the data set is not ideal for fitting of the model.Fitting of the model to operational data sets, in which the measured disturbances were properly excited, were tested but failed due to lack of excitation of the manipulated inputs.One may consider merging several operational data sets, or merging data during grade changes such that both manipulated inputs and measured disturbances are properly exited.However, this approach may also fail because the process itself probably is time varying, and unmodelled disturbances, such as e.g.variations in the raw material, cause drifts and trends which are not accounted for in the model.The fitting and validation reveals deficiencies in the model and perhaps in the experimentation, although it is not clear whether one can eliminate these deficiencies without increasing the order of the model.Increased order models which were fitted to experimental data performed better than the lower order models, thus the choice of a third order model is a trade-o between e.g.complexity and accuracy.
A first order empiric model is added to the mechanistic model to capture neglected and unknown dynamics in the process.The resulting fourth order hybrid model gives much better validation results than the pure mechanistic model for the basis weight, while not a clear better/worse answer for the paper ash and wire tray concentration.The basis weight is perhaps the most important quality variable of the three outputs, and thus we choose to use the hybrid model in the sequel.
A "quasi extended" Kalman filter, in which the Kalman filter gain matrix is constant, is then identified from operational data.Validation of the hybrid model with Kalman filter on operational data seems to give good results.Finally the model is validated with respect to prediction ability during sheet breaks, and prediction ability during grade changes.The prediction ability during sheet breaks seems good for paper ash, but is more uncertain for the basis weight which have a tendency to drift away in some cases.Without improving the prediction ability during sheet breaks it seems that one might as well freeze the inputs at the value prior to the sheet break.However, one may consider using the predictions for operator support if e.g.changes in machine velocity and other variables occur during a sheet break.The prediction ability during grade changes seems reasonably good for the basis weight and paper ash, however for the wire tray concentration, the prediction ability seems poorer.
Let us finally point at a few topics which could be of interest for further research: • Refining the model, especially focusing on improving the predictions of the wire tray concentration.
• Identify and validate covariance matrices, linearize the model, for implementation of extended Kalman filter.
• Online estimation of key parameters, i.e. an augmented Kalman filter.
• Compare with more traditional handling of measurement loss.

Figure 1 :
Figure 1: Manipulated inputs, and outputs of the model.

Figure 3 :
Figure 3: Step responses at the outlet of a pipeline (40 m length, 0:7 m 2 cross sectional area, 2500 kg = s mass flow).Discretization carried out with various numbers of ideally stirred volumes.

Figure 4 :Figure 5 :
Figure 4: (1) The polymer adsorbs onto anionic impurities in the stock.(2) Flocculation caused by the polymer.Note the tight adsorption of the particles and the polymers due to the charge and weight characteristica of the polymer.

Figure 7 :
Figure 7: The responses of the full scale model (solid lines) and of the 38 th order simplified model (dashed lines).

Figure 8 :
Figure 8: Sketch of PM6, according to the third order mechanistic model.

ẋ = f 1
(x; u; d; z; µ) (35) z = f 2 (x; u; d; z; µ) y = g(x; u; d; z; µ),where z R 31 is a set of algebraic equations.The states are Re;filler is the concentration of filler in the reject tank, C W W T;filler is the concentration of filler in the white water tank, and C Dr;fiber is the concentration of fiber in the deculator ("right" side)

Figure 9 :
Figure 9: Experiment plan for the PM6 process operators.The plan shows the changes in setpoints carried out for the data set collected at February 28, 2001.

Figure 12 :
Figure 12: Third order mechanistic model fitted to experimental data.Experimental data collected at February 28th 2001.

Figure 13 :
Figure 13: Validation of third order mechanistic model by simulation.Operational data collected during March 8-11, 2001.Simulated data are bias corrected.

Figure 14 :
Figure 14: Fourth order hybrid model fitted to experimental data.Experimental data collected at February 28th 2001.

Figure 14 Figure 15
Figure14shows the fitted deterministic hybrid model.The improvement compared to the pure mechanistic model seems quite large.RMSE values for fitted fourth order hybrid model are shown below:

Figure 15 :
Figure 15: Validation of fourth order hybrid model by simulation.Operational data collected during March 8-11, 2001.Simulated data are bias corrected.

Figure 16 :
Figure 16: Fitted fourth order hybrid model.Only those elements in the empiric model of Equation 50, corresponding to measured disturbances (elements in E and F matrices), were tuned.Elements in other system matrices and parameters in mechanistic model are as identified with experimental data set collected at February 28, 2001.Operational data collected during March 8-11, 2001.Simulated data are bias corrected.
shows the fitted model and measured outputs.RMSE values for fitted fourth order hybrid model with operational data collected duringMarch 8-11, 2001, are:    Basis weight Paper ash Wire tray cons.sets were used for validation of the new hybrid model.The first validation data set is collected during May 11-16, 2001, and the bias corrected results can be seen in Figure17.There appear to be problems with the basis weight and especially the wire tray concentration outputs, however compared to the model outputs when using the mechanistic model or the hybrid model with d R 2 , the re-

Figure 19 :
Figure 19: Identification of "quasi extended" Kalman filter for fourth order hybrid model.Operational data collected during March 08-11, 2001.Predictions are bias corrected.
Figure20and Figure21show validation results using the hybrid model with "quasi" extended Kalman filter.All in all, the Kalman filter seems to work properly for the two validation data sets.However, in Figure21,

Figure 22 :Figure 23 :
Figure 22: Three periods of simulated sheet breaks.The first and second column show simulated sheet breaks during normal operation, while the third column show a simulated sheet break just prior to a grade change.Operational data collected during May 11-16, 2001.
Figure 21: Validation of "quasi extended" Kalman filter for fourth order hybrid model.Operational data collected during May 19-23, 2001.Predictions are bias corrected.