Parameter Estimation of a Pulp Digester Model with Derivative-free Optimization Strategies

The work concerns the parameter estimation in the context of the mechanistic modelling of a pulp digester. The problem is cast as a box bounded nonlinear global optimization problem in order to minimize the mismatch between the model outputs with the experimental data observed at a real pulp and paper plant. MCSFilter and Simulated Annealing global optimization methods were used to solve the optimization problem. While the former took longer to converge to the global minimum, the latter terminated faster at a significantly higher value of the objective function and, thus, failed to find the global solution.


INTRODUCTION
Pulp and paper industry plays an important role in world economy. The pulp production line is quite complex as it involves thousands of variables that have to be monitored and controlled. The key unit of the pulp mill is the digester, a pressurized upright cylindrical reactor whose height can reach 60 m [1]. Its very sui generis design and physical/chemical characteristics of the process require sophisticated control strategies in order to ensure safe and economically viable operation. Since these strategies rely on mathematical models that predict the unit behaviour, their accuracy determine the steady-state and dynamic performance of the closed-loop system.
The digester, represented schematically in Figure 1, presents a set of peculiarities that make its modelling a quite challenging task [2,3]. It is a heterogeneous moving bed reactor that operates in continuous mode and where  the organic wood components suffer complex chemical reactions under the combined action of inorganic reactants at adequate temperature. Simultaneously, and since different phases co-exist in the vessel, significant mass transfer phenomena also occur. The digester is a distributed parameters system, being the model a set of nonlinear partial differential equations (PDE) [3] that describes the behaviour of the wood chips as well as of the entrapped and free liquid in terms of composition and temperature. The digester is equipped with several sets of screens in order to extract part of the free liquid at different axial locations and a corresponding system of re-injections. Moreover, these recirculation mechanisms are integrated with a multiple-feed system of fresh inorganic reactants. Besides promoting the mass and heat transfer between the chips and the surrounding liquid, the recirculations allow to heat and makeup for inorganic reactants depletion due to chemical reaction. The digester's design is such that co-and counter-current regimens coexist inside the vessel (in the upper and lower halves, respectively). From a mathematical point of view, the described characteristics render the model discontinuous [4,5] and non-convex. Furthermore, the functions describing the chemical reactions are highly non-linear and also discontinuous.
Notwithstanding the fact that mechanistic models are conceptually more complex, they are able to capture the process behaviour in a wider range of operating conditions as compared to the so-called black-box models. The model used herein expands the continuous digester dynamic model of Fernandes [3] in the following ways: (i) the digester parts between the injections and the corresponding extraction points of the co-current part were modelled considering complete radial homogeneity of the characteristics instead of the virtual conical surface eventually induced by the liquid displacement in the recirculation zones; (ii) the discontinuities exhibited by the kinetic equations (but not the ones associated with the free liquid concentrations and temperature) were smoothened by the technique of Abbo and Sloan [6]. These differences call for a new set of parameters in order to predict accurately the real digester behaviour. However, it is expectable that the new set of parameters is not very different from that encountered by Fernandes [3] as the differences in the models are light.
The objective of this work is to carry out parameter estimation in order to match the updated model predictions to an industrial dataset that represent typical operating conditions. The parameter estimation problem is solved via the resolution of an optimization problem that takes into account the specificities of the problem. First, due to the complexity of the mathematical equations and the presence of discontinuities, the use of derivative free optimization method is preferable. Second, it is desirable to find the parameter set that corresponds to the global minimum of the plant-model mismatch. Based on these considerations, two global and derivative-free optimization methods were used: the Simulated Annealing (SA) and the MCSFilter. While the former uses a heuristics for accepting trial points with higher values of the objective function [7], the second is based on a multistart strategy with a coordinate search filter method as the local procedure [8].

PROBLEM FORMULATION
The set of PDE that form the model is discretized along the digester's axial coordinate taking into consideration the geometric specificities of the unit and applying the finite differences approximation. This procedure originates 750 ordinary differential equations (ODE),Ẋ = f (X, U, P). For the purpose of this work, the steady-state condition of the digester is assumed, that is, f (X, U, P) = 0. Finally, a set of output variables, y, are computed from the state variables X by linear functions.
The quality of the pulp is characterized by a number of factors. Of those, the most important is the kappa number that represents the amount of lignin that remains in the wood and is a direct measure of the quality of the pulp. In practice, it should have values between 14.5 and 15.5 for kraft pulp, because below this range the resulting pulp has low fiber strength that affects the paper production and because values above this range lead to a increase in the consumption of reactants in the bleaching stages. The economic performance of the process is measured via the yield that has a typical range between 53% and 54%. Both kappa number and the yield provide an indication of the degradation of the organic matter. Higher reaction rates correspond to further degradation and thus lower values of yield and kappa number. Besides the kinetic parameters, there is a variable that influences significantly the kinetic rates which is the concentration of the alkali (one of the reactants present). Its spacial evolution along the digester is inferred by the measurable concentration at the extractions, that are represented as C E,z where z represents the extraction zone in question. For the exposed reasons, a subset of seven of the output variables was selected to evaluate the model performance, namely the kappa number, the yield, and the alkali concentrations of the extracted flows at all recirculation screens, that is, y = κ η C E,C5 C E,C6 C E,ITC C E,C8 . Among all the model parameters, a subset P of five parameters associated with the chemical reaction was selected for optimization: the stoichiometric coefficients for the alkali consumption in the reaction and three of the pre-exponential factors of the reaction rates, that is, P = ζ L ζ C,H λ H,i λ H,p λ L,p1 . The optimization objective is to determine the values for this set of parameters that allow to minimize the deviations between the model predictions and their corresponding values based on the real practice with the industrial digester. The parameters to be found should lay within pre-defined ranges in order to warranty their physical meaning.
The degree of deviation between predicted and real values can be quantified by the sum of the square difference of the previously selected output variables of the process, y, which are directly dependent on the state variables X.
Therefore, mathematically the optimization problem consists of: where J is the objective function, P is the vector containing the parameters to be estimated, ν i is the weighting factor of the contribution of each variable, y exp,i represents the experimental (industrial) values for output variable i, y i (X) designates the prediction value for output variable i while the subscripts LB and UB stand for lower and upper bounds, respectively.

RESULTS AND DISCUSSION
The parameter estimation problem (1) is solved using the SA and the MCSFilter methods. The weighting factors in the cost function (1a) assigned to each output variable are given in Table 2 and were chosen taking into account the relative importance of each of the output variables. The initial guess of the model parameters, their corresponding lower and upper bounds in (1b), and the obtained solutions are given in Table 1. The initial guess for the model parameters coincide with the same initial guess of the previous study of Fernandes [3]. The lower and upper bounds were set to avoid excessive deviations that could jeopardize the kinetic model used. Table 2 shows the industrially relevant output variable values experimental and predicted with the set of the obtained minimizers. The most pronounced differences between the two optimal solutions are in the parameters ζ L and ζ C,H (see Table 1). These two parameters have a strong direct influence on the consumption of the reactant alkali and therefore affect the wood degradation rate. This has a significant impact on the output variables. On the other hand, for the remaining parameters the difference between the two estimates is small. These parameters affect essentially the rate of degradation of the organic compounds, which is accounted at the bottom of the digester by the kappa number and yield variables. From Table 2 it follows that the solutions obtained by the MCSFilter method slightly outperforms the one of the SA method in what concerns the most important output variables (kappa number, yield and the concentrations of the cooking recirculations).
The moderately better solution provided by the MSCFilter requires a higher computational cost. Nevertheless, at this point and for the purposes of solving the problem once this does not constitute a limitation in terms of applicability of the method. Globally, the parameter estimation from the MCSFilter provides slightly better results.

CONCLUSIONS
Two derivative-free optimization strategies are applied to estimate a set of parameters of a steady-state model of a pulp digester: the Simulated Annealing and the MCSFilter. In what concerns the optimal objective function, the MCSFilter moderately outperforms the SA method. However, the MCSFilter method needs considerably more computational time. Since no derivative information is required, the application of these optimization techniques is of great practical interest, in particular when the computational time is not a critical aspect. As future work, a more detailed characterization of the model parameter space will be carried out to investigate how the starting point affects the performance of the optimization framework.