1D and 2D parameter estimation of a rotary drum filter based on RLS and IV methods

ABSTRACT This work deals with a physical one- and two-dimensional (1D and 2D) parameters estimation of a filtration process of slurry, the second stage of phosphoric acid manufacture. This study focuses on recursive least square and instrumental variable techniques applied to the (1D) and (2D) models. The model of the rotary drum filter is based on different physical laws involved in the filtration phase in order to get a simulator of the filtration process. Besides, many physical parameters rise in the system model and effect enormously the efficiency which should be modelled with precision such as permeability, porosity and viscosity. We use a constructive realization procedure for (2D) systems which may lead to a Fornasini–Marchesini local state-space model to describe the dynamic of the system states.


Introduction
In recent years, membrane filtration technology and particularly low pressure membrane technology have been widely applied in liquid aspiration fields more and more [1].
For some time, it is well known that membrane separation processes are considered very beneficial, advanced and efficient technologies to use in factories. These methods are used to separate suspended particles. The filtration technique is classified according to the size of the pores. So, there are four types of filtration: macrofiltration, microfiltration, ultrafiltration and nanofiltration [2]. Filtration processes can be classified into two categories: • Discontinuous: The press filter and the Nutshe filter are most common. • Continuous: The rotary drum filter and the pass band filter are most common.
In our case, we are interested in the rotary drum filter which ensures the separation between the phosphoric acid and the phosphogypsum. Few works are dealt with the modelling of the filtration phases due to its complexity and the number of parameters interfering in the filtration process. The problem is therefore to model the system to be studied so that the model is as faithful as possible to the actual behaviour of the system.
Modelling filter is essentially based on the Darcy's law treated by [3] that applies to an homogeneous and isotropic porous environment traversed by a flow at low speed in order to calculate the rate of phosphoric acid produced. Darcy's law is processed also by [4], a similar strategy has been discussed by [5]. Then a relationship frequently quoted was proposed by [6] and modified by [7].
In this research, we describe the dynamic behaviour of a rotary drum filter by an homogeneous twodimensional (2D) model. Fornasini and Marchesini [8] focus on the importance of (2D) systems. (2D) systems have attracted attention because of their wide range of practical projects and significant significance. (2D) system theory can be used efficiently in many fields such as digital filtering and batch processes [9,10]. For instance, investigations on (2D) systems have been related to state-space model realization and stability [11], estimation and observers [12][13][14] and filtering [15][16][17].
More recently, many other interesting contexts where (2D) systems prove to be the appropriate setting for carrying on a thorough and successful analysis have been illuminated [18,19]. Indeed, in the (1D) context, there has been a long stream of research, which originated in the 70s and flourished in the 80s [20,21], but still represents a very important research topic. During the last years, research focused on the (2D) system theory and it have been interested in most of the classical topics already investigated in the (1D) setting [22][23][24][25]. Generally, the main interesting point of the (2D) systems is its dependence on two variables which are the time and the space. However, there are many (2D) systems such that none of their variables is time or spatial [26][27][28][29]. Moreover, one of the fundamental issues in the (2D) systems theory is the realization of a given transfer function or transfer matrix by a certain kind of the (2D) local state-space model, typically by Roesser model or Fornasini-Marchesini second (FM-II) model.
Our contribution is mainly based on the modelling of a rotary drum filter using physical laws and particularly concerned with the realization problem for a given (2D) MIMO system by the (FM-II) model. As well as a (2D) least squares identification is applied and it has a great importance in a huge range of applications [30]. However, the special structure of the (2D) data set gives rise to the development of cost-effective algorithms for the sequential determination of parameters [31,32].
This paper is organized as follows: Section 2 provides the characteristics of the rotary drum filter subject of steady with its different running modes. Section 3 describes the (1D) nonlinear modelling of the filter using several laws. Section 4 presents the linearization of the (1D) system. Section 5 describes the implementation of (FM-II) model on our system and the algorithms used. Section 6 represents parameter estimation with the recursive least square (RLS) method in the case (1D) and (2D ). Finally, concluding remarks are made in Section 7.

The rotary drum filter
There are mainly two vacuum-operated devices that are rotary drum filters and band filters. They have the same applications but the band filters deal thicker slurries (50% solids). Among the continuous filters which are considered, we find the rotary drum filter. The rotary drum filter is constituted by two coaxial cylindrical drum and a outer drum carries a filter canvas.

Characteristics
Our filter has a huge toric surface which has a value close to 262 m 2 . This filter ( Figure 1) is composed by 36 flat cells covered by 36 membranes with 36 tubes. Each tube is connected to a cell, a vacuum box and a driving motor which is responsible for the transition of the slurry by five sectors which are: the pre-sector, the strong acid sector, the acid medium sector, the weak acid sector and the sector of wash canvases. Figure 2 shows a sketch [33] of the filter. The outputs of the filter are the strong acid which goes to storage, the gypsum which goes to outer area and the medium acid which returns to the reactor [34]. Tables 1 and 2 gathers the different characteristics of the filter.
The cake (slurry) reaches a certain thickness and by the end of the cycle, the gypsum is removed from the filter by scraping system. A cycle of wash and spin is often assistant. These filters provide greater investment but they have a lower cost of operation: they are therefore suitable for large productions.   In continuous filtration, a cake stand (or filter media) moves slowly to the suspension. When the cake is thick enough (or clogging of the mass is large enough), we shall remove the cake with, where appropriate, an intermediate steps of washing and dewatering of the cake.

Running modes
The membrane filtration can be used in two main operations: frontal filtration (dead end) and tangential flow   figure 3). These two modes are very important and correspond to two technologies and two completely different approaches to filtration. In our case, we have a frontal filtration because the slurry is injected in a perpendicular way in the filter. For this mode of running, the fluid flows vertical to the membrane by causing a shear which allows the accumulation of material. All the material retained accumulates on the membrane. Nevertheless, this type of filtration increases a little the occurrence of clogging phenomenon. The frontal filtration is the best solution but it requires a washing step to avoid the phenomenon of clogging. In general, the filtration can be operated at constant pressure or constant flow. In our case, we have a constant pressure. The pressure gradient may be generated by the operation of a pump that circulates the liquid above the membrane. When the pressure of the vacuum pump decreases, the thickness of the cake increases along the axis of the membrane [35]. Due to the coriolis effect, the thickness of the cake has tendency to move to the outer edge which decreases the height of the edge. This phenomenon is neglected in our case because the speed of the rotation of the filter is so low.

(1D) nonlinear modelling of the filter
The literature distinguishes two types of approach to model filtration: macroscopic models based on the material conservation equations and a filtration rate determined empirically as well as a microscopic models describing the grain transportation process. Concerning macroscopic models, these types are not based on a physical modelling of particle deposition mechanisms. The parameters of these models do not necessarily have a clear physical meaning but its results are usually close to the experimental results.
Developing a model requires first to write mathematical equations that relate state variables descriptive of physical process. Compared to other procedures, fixed cultures, the filtration system has the complementary particularity to include an additional model of physical filtration, which increases the difficulty of resolution and integration of system equations. For activated sludge processes like the filter, the problem of the complexity of the models make them control possibilities in practice very limited cases [36].
The model we are looking for is a model for a control objective. So we will look for the physical equations that contain the state variables X and control variables U. The state variables are the height of the slurry in the filter which denoted h and the volume of the filtrate acid denoted V . The two actuators involved in the process are the pump providing slurry which is regulated to a prescribed volumetric flow q b , and the motor for rotating the filter which is regulated to the described angular speed w f .

The height of the cake h
According to the mass balance, the quantity of the slurry in the filter is equal to the difference between the quantity of slurry injected in the filter and the quantity of acid gotten from the slurry: where M b is the mass of the slurry injected into the filter, M filter is the mass of the slurry after filtration and M A is the mass of the phosphoric acid filtrated. The mass of the slurry M b injected from the reactor to the filter is modelled by the volume of the slurry V b injected into the filter: where ρ b is the volumetric density of the slurry and V b is the volume of slurry filled the nacelle. The mass of the slurry M filter after filtration is modelled by the height of the cake h and the surface S of the nacelle: The mass of the filtrated acid M A is modelled by the volume of this acid V : where ρ A is the volumetric density of the phosphoric acid P 2 O 5 .
Replacing (2), (3) and (4) into (1), we get If we divide (5) by ρ b and we suppose that ρ A /ρ b = m ρ , we find the following formula: Then the filter has the shape of a ring with external radius R 1 and internal radius R 2 as shown in Figure 4. So the surface S of the nacelle is modelled by the angle θ: If we introduce (7) into (6), we obtain If we derivate (8) with respect to time, we obtain (9) In fact, the angular deviation with respect to time dθ/dt is equivalent to the angular speed of the filter w f , and the slurry volume deviation with respect to time dV b /dt is equivalent to the volumetric flow of the slurry q b , so we obtain from (9) the following equation: (10) We take dV/dt = q AF . From (10), we can get the height deviation with respect to time dh/dt: (11) Finally, we get

Volume flow rate of the phosphoric acid q AF
To describe the flow of a fluid through a porous medium, we use the Navier Stokes law and the Darcy's law. The Darcy's law is a relation between the flow of a fluid and the pressure. Permeability is a physical characteristic that represents the ease of a material to allow the transfer of fluid through a connected network. The Darcy's law makes it possible to connect a flow rate to a pressure gradient applied to the fluid thanks to a characteristic parameter of the crossed place: the permeability k. The mass conservation equation coupled to the fluid motion equation yields the incompressible Navier Stokes equations and the main unknown parameters in this equation are the mass density, the velocity, the pressure and the temperature, but this list could be longer depending on the case being studied. The purpose of this part is to describe the corresponding fundamental solution for the equations that model an incompressible fluid in the exterior of a rigid body that is rotating at a constant velocity and rotation [37]. The Navier Stokes equation is: with ϑ is the velocity field, P is the pressure, ρ is the fluid density, μ is the kinematic viscosity, t is the time and is Laplace operator. = 1 j=1 ∂ 2 /∂z 2 j stands for the Laplace operator.
If we multiply (13) by S, we obtain If we take ϑ · S = q AF , so we find Now we should calculate the following two terms: To find these two terms, Darcy's law is used which has the following formula: The strong acid flow (28%P 2 O 5 ) is characterized by a ϑ, according to the Darcy's law, is proportional to the pressure P and the permeability K and inversely proportional to the dynamic viscosity μ of the fluid.
First, we will express the pressure drop across the filter. The term is derived from the Darcy's law.
If we derivate the previous equation, we obtain If we derivate (18), we get (19) Replacing (18) and (19) in (16), we find The overall pressure P in the slurry mixture is the sum of the pressure of the solid particles P s and the pressure of the fluid P l such as P = P s + P l .
The differential of the pressure with respect to the z axis: The expressions of P s and P l are given by the following relations: where Q e is the speed of strong acid (ms −1 ); ρ s is the solid density (dimensionless); ρ l is the liquid density (dimensionless); μ l is the liquid viscosity (Pas.s ); E is the ratio of lateral and normal stress (dimensionless ).
So, the overall pressure is the sum of the liquid and solid pressure: If we replace ∂P(z)/∂z by its expression, we obtain finally, we get

The (1D) nonlinear model
The nonlinear model obtained from (12) and (24) is as follows: For controlling such systems, it is quite necessary to take into account the nonlinear phenomena. The simulation of the nonlinear system gives the two following outputs: • The height of the cake, h (cm).
• The volume flow rate of the strong acid, q AF (m 3 · h −1 ).
We note that the thickness of the cake increases over time until it reaches 6 cm. The volume flow of the filtrate increases to a value close to 74 m 3 · h −1 . These results are close to the real system.
with h 0 is a small variation of the height of the cake. So, we obtain We suppose that: So the linear model of the rotary drum filter is given by the two following equations: The equations of the model contains the product of a state variable and a control variable. In order to linearize the model, we must use the Jacobian linearization. We will select an operating point at steady-state P(X,Ū) such as f (X(t),Ū(t)) = 0. Moreover, if we suppose that U = 0, so the operating point is: The model is as follows: where F x , B u are the Jacobian matrices of the partial derivatives of f (X, U) with respect to X and U, the state and the input variable of the system, respectively, which measured at point P(X,Ū). The two matrices F x and B u change according to the operating point. At steady-state, we havē So, the final Jacobian matrices of the state-space are Validation of the model is the process of determining the degree to which a simulation model and its associated data are an accurate representation of the real world from the perspective of the intended uses of the model. In order to investigate the filter and its changing states, our system was modelled and simulated by "MATLAB" which is among the most developed software for the simulation of production systems. We suppose that the rotary filter on the TCG factory is a multi-input multi-output nonlinear system (See Figure 5), which has as inputs U1 flow rate of the slurry and U2 speed of the filter, and outputs Y1 height of the cake and Y2 the volume of the filtrate.
Also we have The Figure 6 has the output Y 1 = Y r , which corresponds to the thickness of the cake as well as the output of the model Y m are almost constant. The relative error between the output of the real system Y r and the output of the model Y m belongs to ±20% (see Figure 8) The mean squared error is MSE = 0.4%. The Figure 7 describes the output Y 2 = Y r , which corresponds to the flow rate of the filtrate and the output of the model Y m are almost constant. The relative error between the output of the real system Y r and the model Y m belongs to ±0, 47%, see Figure 9. The mean squared error is MSE = 0.313%. We note that there is a gap between the real output and the output of the model, this difference is due to the variation of the  physical parameters of the system during the filtration phase.

2D modelling of the rotary drum filter based on (FM-II) method
Consider the linear 2D system described by the following (FM-II) model: where X(d, k) ∈ R n , U(d, k) ∈ R l and Y(d, k) ∈ R m are, respectively, the local state, input and output vectors, d is the time t and k is the angle 1/θ and A 1 , A 2 , B 1 , B 2 and C are real system matrices of suitable sizes. The system is also conventionally denoted by (A 1 , A 2 , B 1 , B 2 , C).  The following lemma gives a sufficient condition for the asymptotic stability of the (2D) system in term of an LMI.

Lemma 5.3 ([39]):
The (2D) system is asymptotically stable if there exist matrices P 0, Q 0, R = R T = 0, 0, F ∈ R 2n·n and G ∈ R n·n such that: The different stages of the (2D) modelling are: (1) The transfer matrix of the previous model is: where ψ i is the column vector defined as: and diag is a diagonal matrix. (2) Write D(z 1 , z 2 ) and N(z 1 , z 2 ) in the form of: where: We have F(z 1 , z 2 ) = N(z 1 ,z 2 ) where D HT and N HT are real matrices with sizes conformable to ψ.
(3) Construct the matrices A i0 , B i ; i = 1, 2 such that: The realization is finally obtained as: By thoroughly investigating the structural properties of (FM-II) model, using the relations (34), (35) and (36), a realization can be obtained.
ψ has to satisfy the following conditions: (1) Condition 1: ψ i contains all the power products occurring in the polynomial entries of the i th column of F(z 1 , z 2 ). (2) Condition 2: ψ i contains z 1 or z 2 , or both z 1 and z 2 .
The above realization procedure produces a state-space description (A 1 , A 2 , B 1 , B 2 , C) which have the following properties: Starting from the initial ψ i (i = 1, . . . , l), the following algorithm is given below for constructing the final ψ i (i = 1, . . . , l) which has to satisfy the conditions of (2) and (3). Once ψ 1 , . . . , ψ l are constructed by Algorithm 1, ψ = diag{ψ · · · ψ}, can be readily obtained. ψ is now of size n × l, where n = l i=1 n i . Note that n will be smaller than n' if some power products are absent in the polynomial entries of F(z 1 , z 2 ). This will be illustrated by a non-trivial example in the next section. Next, express D(z 1 , z 2 ) and N(z 1 , z 2 ) in the form D(z 1 , z 2 ) = D HT ψ, N(z 1 , z 2 ) = N HT ψ, with D i j ∈ R 1n j and D i j ∈ R 1n j are row vectors whose entries are the coefficients of the (i, j)− indexed polynomial in D(z 1 , z 2 ) and N(z 1 , z 2 ), respectively. The system matrices A 1 , A 2 , B 1 , B 2 and C can then be constructed by the Algorithm 2.
Calculation of the transfer matrix M The state vector of the rotary drum filter has the following form: If the answer is yes, return to Step 2, and if the answer is no, go to Step 7. If j = n i , proceed to Step 6. 6: r = r + 1. Check whether there exists an entry . If the answer is yes, insert either z r−s 1 ψ i (k) or z r−s 2 ψ i (k) at an appropriate position according to the descending order of power products in ψ for s = 1, ..., r − 1. Let n i = n i + (r − 1) and return to Step 3. In the case that the answer is no, check whether r ≺ deg ψ(j) 2 . If yes, repeat Step 6. Otherwise, proceed to Step 7. 7: Check whether steps 1 and 2 and that , into ψ i according to the descending order of power products. n i = n i+1 . 8: return to Step 3. with: The command vector of the rotary drum filter has the following form: (38) with: and it can be shown, in the same way of [40], that Throughout this part, we denote by p the Laplace operator.
At first, we calculate the determinant of F x : Using the previous determinant, we calculate the adjugate matrix of F x to obtain the final transfer matrix M(t, 1/θ ). with We have p = z 1 and 1/θ = z 2 so, From the transfer matrix, we can formulate the following matrices: Thus the column degrees of columns 1 and 2 in F(z 1 , z 2 ) are k 1 = 1 and k 2 = 2, and the power products occurring in the polynomial entries of the columns with non zero coefficients are {z 2 , z 1 } and {z 2 , z 1 , z 2 1 , z 2 z 1 }, respectively. The initial column vectors ψ 1 and ψ 2 satisfying condition (1) are: ψ 1 = (z 1 z 2 z 2 z 1 ) T ; ψ 2 = (z 1 z 2 z 2 1 z 2 z 1 ) T Applying now the 1, ψ 1 and ψ 2 should satisfy conditions (2) and (3).
So, for ψ 1 a new term z 2 1 has to be inserted into ψ 1 , so, ψ 1 = (z 1 z 2 z 2 1 z 2 z 1 ) T . This new ψ 1 now satisfies conditions (2) and (3); therefore, we have For ψ constructed above, it is ready to calculate that

RLS for the (1D) system
RLS algorithms for parameter estimation are usually described in relation to the batch-type least squares algorithms, minimizing a quadratic error criterion. This error is the difference between the estimated and the observed values of the process output [41]. Consider the following linear process with slowly timevarying parameters, described by the following difference equation: The signal y k is the process output and x k is the input. For a notation most practical a parameter vector is defined: t = (b 0 , b 1 , . . . , b M , −a 1 , −a 2 , . . . , −a N ) and a signal vector:u t = (x k , . . . , x k−M , y k−1 , . . . , y k−N ) u k , ∈ R p , p = M + N + 1 (44) can be rewritten as follow: The parameter vector of the filter is estimated by an m-vectorˆ . An error criterion is established: where η is an exponential forgetting factor η ∈]0, 1] and 0 contains initial values ofˆ 0 . Q ∈ R m×m is a positive definite matrix. Note that usually in the normal RLS algorithms the first term on the right-hand side of equation (51) is omitted By setting ∂J/∂ˆ = 0 it is readily calculated that In this equation, a matrix P k is introduced where which is a recursive algorithm; P 0 = Q P 0 = Q using P k defined in (48), (50) can be written as follows:ˆ with P 0 = Q, it is seen thatˆ 0 = 0 . Substituting (50) yieldsˆ .ˆ 0 = 0

Simulation results for the (1D) linear system
At the end of the modelling, we obtained a transfer matrix consisting of four transfer functions M 1 , M 2 , M 3 and M 4 . These transfer functions are arranged as in Figure 10. From the filter-state representation, we got a transfer matrix which is as follows: with The performance of the RLS method was tested by computer simulation. Now, we want to look for the estimation vectorsθ 1 ,θ 2 ,θ 3 andθ 4 which correspond respectively to the transfer functions M 1 (p), M 2 (p),   Tables 3, 4, 5 and 6: Using these tables, we notice thatθ and θ are close, which validates the performance of the (1D) RLS identification method.

RLS for the (2D) system
This section is concerned with a basic spatial RLS algorithm for a second Fornasini-Marchesini (2D) model. Either, x(d, k) is the input and y(d, k) is the output of the system.
L y and L x denote the support region for the ARX model. Given an input (2D) signal x(d, k) and a desired responsex(d, k), (h, k) ∈ [0, 0].[N 1 , N 2 ], the optimal least square (2D ) filter is obtained by minimizing the cost function.
where e(t 1 , t 2 ) =x(t 1 , t 2 , k) − y(t 1 , t 2 ) is the instantaneous error and η ∈ [0, 1] is the so called forgetting factor. In this section, (2D) support region of general shapes are considered for both L y and L x . More precisely, L consists of a union of intervals Using the previous notation, equation (49) takes the form Let us define the row-wise data vectors and the system parameters that are attached to the support regions L y and L x . We use the following notation k y 2 = 1 (1) Output data corresponding to L y : (2) Input data corresponding to L x : Then (59) can be expressed as a linear regression The regressor vector (h, k) is defined in terms of input-output data where Y L y (d, k) = y l y (i) i=k y 1 ···s The parameter vector carries the system coefficients The least square algorithm is as follows: 2D parameter update (d, k + 1) = (d, k) + w(d, k + 1)ε(d, k + 1), ε(d, k + 1) = e(d, k + 1)/β(d, k + 1), 2D Kalman gain update w(d, k + 1) = −η −1 P(d, k) (d, k + 1), P(d, k + 1) = η −1 P(d, k) + w(d, k + 1)w t (d, k + 1) β(d, k + 1) .

Simulation results for the (2D) linear system
We consider thatθ (2D) are the estimated vectors, respectively, of the transfer functions M 1 , M 2 , M 3 and M 4 . We apply the (2D) identification algorithm in order to obtain the following tables: From the simulation results presented in Tables 7, 8, 9 and 10, we note that the estimated parameters tend to the calculated parameters with a tolerance of ±10%. The parameters in (2D) model are closer to calculated values than those given by a (1D) model. Since the convergence of the parameters is done, we note that the (2D) (RLS) identification method is more efficient than the (1D) identification. For example for the (1D) model we have a 1 = −0.2331 whether for the (2D) model a (2D) 1 = −0.2530. So, we note that a (2D) 1 is closer to the calculated value equal to −0.2501. So the convergence of parameters is performed and the (RLS) identification method is so effective in the (2D ) case.    Table 9. Estimation vector ofθ

Conclusion
In this paper, (2D) model was proposed to describe the dynamic behaviour of a rotary filter, which is dedicated to the aspiration of various types of acid: strong, medium and weak. Therefore, its study requires a deep research in order to know all the parameters that have an influence on the normal operation of the filter. Using the parameters that characterize the "AOUSTIN-UCEGO" filter, we get to find a two dimensions model for the filter. A constructive state-space realization procedure has been used for (2D) systems which may produce (FM-II) local state-space model. An RLS identification is made in both cases (1D) and (2D) and from which it is noticed that the results of simulation are near to the real values. The RLS identification for the (2D) system is more efficient since the parameters in this case depend on the temporal and spatial variables.