Real-Time Optimization of Pulp Mill Operations with Wood Moisture Content Variation

In tropical countries, such as Thailand, the variation of tree moisture content can be significant based on seasonal variations in rainfall. Pulp mill operation optimization accounting for wood moisture variation was used to determine optimal operation conditions and minimize production cost. The optimization models were built using empirical modeling techniques with simulated data from the IDEAS software package. Three case studies were performed. First, a base case of nominal annual operation at a fixed production rate was used to calculate production cost that varies with wood moisture content. The second case is annual optimization where production was allowed to vary monthly over an annual cycle to minimize production cost. For the third case, real-time optimization (RTO) was used to determine optimal production rate with the wood moisture content varying every 3 days. The rolling horizon approach was used to schedule production to keep inventory levels within bounds and with a penalty applied to deviations from the annual expected values of inventory. The advantage of RTO in accounting for moisture content variation was confirmed by annual production costs results simulated for 20 years. These results statistically demonstrated that the overall cost was reduced compared to the second case of monthly production targets.


Introduction
The pulp and paper industry is an important and growing sector in the forest industry of Thailand for which the majority of fiber inputs are produced domestically. The pulp and paper sector is energy and raw materials intensive and faces several challenges that include rising feedstock and energy prices and increased concern for environmental impact. Reducing cost is the main objective for the operation of a pulp and paper mill and there is a drive toward efficiency to maximize the profit rate of the whole plant [1]. Optimization models can significantly improve pulp and paper mill operation and many of them focus on overall amount of operating consumptions as raw material, water, and energy [2]. In this paper, a process integration model was developed to improve operation and minimize operating cost, and sub-models of different mill operations were developed to support optimization of plant performance [3].
In process operations involving natural inputs, there is a variability that affects physical properties of raw material, product quality, and operating cost; therefore, on-line optimization can be an important contributor to overall plant efficiency. Specifically, real-time optimization (RTO) is an approach used in various industrial continuous processes and can be applied to the pulp and paper industries. RTO is a methodology to support decisions based on optimizing an economic objective and plant operating their operation have significant data requirements and involve a number of assumptions. The IDEAS software program of ANDRITZ Inc. is a pulp mill simulator used by the pulp manufacturing industry to simulate their processes. In this work, the IDEAS software program was used to simulate the steady-state pulping process and generate the data necessary to fit empirical models that could be used in optimization. operation have significant data requirements and involve a number of assumptions. The IDEAS software program of ANDRITZ Inc. is a pulp mill simulator used by the pulp manufacturing industry to simulate their processes. In this work, the IDEAS software program was used to simulate the steady-state pulping process and generate the data necessary to fit empirical models that could be used in optimization. Various levels of empirical modeling are used to model the whole pulping process. We term each of these levels as white-, black-, and grey-box models. A white-box model for pulp modeling consists of mass and energy balances of sub-models with parameters to represent the performance of equipment [13], such as washing efficiency that is used to calculate water flow for washing and black liquor content in evaporation and to calculate amount of steam to evaporate. This is still an empirical model because the parameters are estimated from data or assumed to take values found in previous studies. This is the type of model built in IDEAS and will serve as the starting point for other empirical models that are cheaper to evaluate and more explicitly link the degrees of freedom for optimization with the variables in the objective function and constraints.
A black-box model predicts output variable from input variable values and uses no details from within the process it is modeling. In this work, the black-box technique is used to build a sub-process model of pulp mill unit operations. While the building of mathematical models of pulping unit operations, a polynomial function was used in modeling that includes three types of terms, linear, quadratic, and cubic, as shown in Equation (1). The final model structure, the β values, is determined from the observed values of variables in process, output variable (Y), and input variables (x) and is a subset of the terms in Equation (1). The specific variables for different operations are shown in Table  1 where coefficient value (β) of each term is from regression.  Various levels of empirical modeling are used to model the whole pulping process. We term each of these levels as white-, black-, and grey-box models. A white-box model for pulp modeling consists of mass and energy balances of sub-models with parameters to represent the performance of equipment [13], such as washing efficiency that is used to calculate water flow for washing and black liquor content in evaporation and to calculate amount of steam to evaporate. This is still an empirical model because the parameters are estimated from data or assumed to take values found in previous studies. This is the type of model built in IDEAS and will serve as the starting point for other empirical models that are cheaper to evaluate and more explicitly link the degrees of freedom for optimization with the variables in the objective function and constraints.
A black-box model predicts output variable from input variable values and uses no details from within the process it is modeling. In this work, the black-box technique is used to build a sub-process model of pulp mill unit operations. While the building of mathematical models of pulping unit operations, a polynomial function was used in modeling that includes three types of terms, linear, quadratic, and cubic, as shown in Equation (1). The final model structure, the β values, is determined from the observed values of variables in process, output variable (Y), and input variables (x) and is a subset of the terms in Equation (1). The specific variables for different operations are shown in Table 1 where coefficient value (β) of each term is from regression. The model selection selects one model from the list of models in which the terms are significant. We use a model selection criterion to find the best fitting model [14]. Hirotugu Akaike proposed the Akaike Information Criterion (AIC) method to select the best approximate model [15]. The general formulation of AIC is given in Equation (2).
where K refers to number of estimate parameters in the model, the optimal fitted model is identified by the minimum value of AIC. Terms in Equation (1) were introduced sequentially from linear to cubic and the minimum AIC selected. In yearly operation, the season affects operating cost because the moisture content of wood changes monthly. This affects the amount of chip, chemical, electrical, and water consumption to produce a ton of pulp. A white-box model is used to calculate amounts of required input variables, which are raw material costs and operating costs of pulping by using the IDEAS simulation that depends on daily production rate (P). In addition, the moisture content (MC) was varied as an input to the white-box model and then this study used an empirical black-box model for creating daily operating cost (C total ) in function of moisture content and production rate; the model is shown in Appendix A. The operating cost is directly related to production rate and is maximized at the highest moisture content of the wood. Therefore, the daily operating cost equation is function of production rate and moisture content of wood represented in Equation (3).
It was necessary to include a bi-linear term to represent the interaction between the moisture content and the production rate. The equation was tested valid for the range of production rates [1000, 2000] air dry metric tons(admt)/day with moisture content of 40-55% and returned errors less than 1%. Extrapolations beyond these ranges for the specific pulp mill modeled should be done with caution.

Inventory Balance Approach
Multi-period production is important degree of freedom for production and operation management of businesses in which the period depends on machine's performance, time, customer's orders, etc. Effective inventory flow management is one factor for profitability due to enterprises requiring enough inventory to satisfy customer demands [16]. On the other hand, too much inventory increases the cost of inventory carrying [17].
In improving manufacturing, the pulp industry is faced with the challenge of producing products at the right time and quantity at minimum costs. Inventory balance approach is a network flow analysis to optimize production quantity for periods in planning process [18]. In any period, the supply for process is the inventory from the prior period (I t−1 ) plus the production in the period P t . This supply can be used to respond to the demand in the period D t or held in inventory as I t . The simplest form of inventory balance is beginning inventory + production = ending inventory + goods sold, as shown in Equation (4).
A series of these equations is used to represent the evolution of inventory over time in response to demands and production decisions. The production decisions are optimized to meet demands while minimizing inventory holding costs.

Real-Time Optimization with Rolling Horizon
The purpose of the optimization in this work is to determine the optimal production rate with minimum total operating cost, under wood moisture content variation, when the variation is dependent on seasonal changes. A discrete time RTO technique is used in optimization employing a rolling horizon approach.

Real-Time Optimization
The plant profitability can be enhanced by optimization of operating conditions and is particularly important in plants with high throughputs of material and energy. RTO is a technique to enable business decisions in real time based on planning models. It involves the on-line calculation of optimal process setpoints and allows the profits from the process to be maximized or costs to be minimized while conforming to operating constraints [19,20]. Typically, the economic cost model of a plant involves costs of raw materials, price of products, and operating cost, which are functions of operating conditions in RTO. The optimization goal in this study is to minimize operating cost of pulping process.
The RTO was used to optimize the process operating conditions, and process set points were updated to the control system. Figure 2 represents the structure of the RTO system. Disturbances such as feedstock variations and other environmental variables related to economics are fed into the closed loop and the RTO model is updated. The supply setpoint is recalculated at each iteration until the minimum deviation between model variables and plant measurements occurs and this establishes the optimal setpoints for operation. The disturbance used in this study is the moisture content of wood, which cannot be measured but varies seasonally. The wood moisture content value is initialized as the average moisture content of the wood, which is an input variable for the pulping process optimization. In the optimization, the setpoints of overall operation that can be manipulated are the water shower in washing, evaporation steam, makeup chemical as Na 2 SO 4 and CaCO 3 , water shower in mud filter, and white liquor flow (WL). The measured outputs of the process, specifically the production rate and liquor flows, will be different given the variation of the wood moisture content from its average value. A function was created that used these two outputs to estimate the moisture content, given the other values of the setpoints. These setpoints were re-optimized using the estimated moisture content of wood using the optimization formulation. This process was iterated until the results from the optimization and simulation become equal.

Moisture Content of Wood Variation
Production management becomes more complex when considering different uncertainty sources in process operations. Types of uncertainty sources can be classified as external sources, such as uncertainty in demand, prices, and availability of resources, and internal sources, such as fluctuations in process parameters and other sources, such as errors of measurements [21]. Pulp mill operations have several uncertainty sources, particularly the moisture content of wood varies in a seasonal pattern with an overlying random variation from wood lot to lot. This affects production capacity of a pulp mill. In this work, the moisture content of wood is represented by a pre-determined average that varies month to month and an actual moisture content that is estimated by solving an optimization problem. The average moisture content monthly depends on wet and dry season in Thailand and ranges from 40% to 55%. However, in day-to-day production, actual moisture content can vary from the average. The optimal production rate is determined in real time due to moisture content variation. In this work, the actual moisture content was assumed to follow a normal distribution about the mean where the standard deviation of the distribution was specified as five percentage points from the average moisture content of each month. The moisture content variation framework is illustrated in Figure 3. The assumed actual moisture content is updated in every three days, which results in 10 periods in a month, as depicted in Figure 3.

Moisture Content of Wood Variation
Production management becomes more complex when considering different uncertainty sources in process operations. Types of uncertainty sources can be classified as external sources, such as uncertainty in demand, prices, and availability of resources, and internal sources, such as fluctuations in process parameters and other sources, such as errors of measurements [21]. Pulp mill operations have several uncertainty sources, particularly the moisture content of wood varies in a seasonal pattern with an overlying random variation from wood lot to lot. This affects production capacity of a pulp mill. In this work, the moisture content of wood is represented by a pre-determined average that varies month to month and an actual moisture content that is estimated by solving an optimization problem. The average moisture content monthly depends on wet and dry season in Thailand and ranges from 40% to 55%. However, in day-to-day production, actual moisture content can vary from the average. The optimal production rate is determined in real time due to moisture content variation. In this work, the actual moisture content was assumed to follow a normal distribution about the mean where the standard deviation of the distribution was specified as five percentage points from the average moisture content of each month. The moisture content variation framework is illustrated in Figure 3. The assumed actual moisture content is updated in every three days, which results in 10 periods in a month, as depicted in Figure 3.

Moisture Content of Wood Variation
Production management becomes more complex when considering different uncertainty sources in process operations. Types of uncertainty sources can be classified as external sources, such as uncertainty in demand, prices, and availability of resources, and internal sources, such as fluctuations in process parameters and other sources, such as errors of measurements [21]. Pulp mill operations have several uncertainty sources, particularly the moisture content of wood varies in a seasonal pattern with an overlying random variation from wood lot to lot. This affects production capacity of a pulp mill. In this work, the moisture content of wood is represented by a pre-determined average that varies month to month and an actual moisture content that is estimated by solving an optimization problem. The average moisture content monthly depends on wet and dry season in Thailand and ranges from 40% to 55%. However, in day-to-day production, actual moisture content can vary from the average. The optimal production rate is determined in real time due to moisture content variation. In this work, the actual moisture content was assumed to follow a normal distribution about the mean where the standard deviation of the distribution was specified as five percentage points from the average moisture content of each month. The moisture content variation framework is illustrated in Figure 3. The assumed actual moisture content is updated in every three days, which results in 10 periods in a month, as depicted in Figure 3.

Rolling Horizon Approach
The rolling horizon approach has been represented as an effective decomposition technique for planning and scheduling problems and broadly used in control methods. The rolling horizon method decomposes sequentially the problem by dividing the time horizon and solving subproblems by fixing the discrete variables of previous time periods [22].
In this work, a rolling horizon approach is used to schedule production, with constraints on inventory levels, to reject disturbances introduced by the moisture content. The rolling horizon period duration is 3 days. An estimate of the moisture content is calculated based on the deviation of the measured production rate from the set point. The production set point is then adjusted to reflect whether the moisture content is higher or lower than the monthly average. In general, production rate will be increased if the moisture content is below the average but the overall optimization will balance inventory holding with production costs over an extended horizon.

Results and Discussions
"fmincon", a MATLAB function, is used as the function minimization routine to minimize production cost and determine optimal pulp production rate. The results are divided into three different optimizations at different time scales and with different levels of information.

Annual Optimization
The optimal monthly production rate of each month in a given year is goal of annual optimization where the average moisture content of wood is known for each month and varies seasonally. The pulp demand (D) for each month was specified as 1600 admt/day but could be varied as well. In normal annual production, production cost is calculated by empirical Equation (3) with constant production rates at 48,000 admt/month. The annual production cost is M$150.19, which is related to the moisture content where the maximum production rate is in the months of low moisture content; these results are represented in Figure 4.

Rolling Horizon Approach
The rolling horizon approach has been represented as an effective decomposition technique for planning and scheduling problems and broadly used in control methods. The rolling horizon method decomposes sequentially the problem by dividing the time horizon and solving subproblems by fixing the discrete variables of previous time periods [22].
In this work, a rolling horizon approach is used to schedule production, with constraints on inventory levels, to reject disturbances introduced by the moisture content. The rolling horizon period duration is 3 days. An estimate of the moisture content is calculated based on the deviation of the measured production rate from the set point. The production set point is then adjusted to reflect whether the moisture content is higher or lower than the monthly average. In general, production rate will be increased if the moisture content is below the average but the overall optimization will balance inventory holding with production costs over an extended horizon.

Results and Discussions
"fmincon", a MATLAB function, is used as the function minimization routine to minimize production cost and determine optimal pulp production rate. The results are divided into three different optimizations at different time scales and with different levels of information.

Annual Optimization
The optimal monthly production rate of each month in a given year is goal of annual optimization where the average moisture content of wood is known for each month and varies seasonally. The pulp demand (D) for each month was specified as 1600 admt/day but could be varied as well. In normal annual production, production cost is calculated by empirical Equation (3) with constant production rates at 48,000 admt/month. The annual production cost is M$150.19, which is related to the moisture content where the maximum production rate is in the months of low moisture content; these results are represented in Figure 4. The objective function is to minimize annual production cost as the sum of monthly costs Cm, where the monthly production cost is a function of monthly average moisture content MCm and The objective function is to minimize annual production cost as the sum of monthly costs C m , where the monthly production cost is a function of monthly average moisture content MC m and production rate P m . An inventory balance is enforced where I m-1 and I m refer to the inventory levels, above the safety stock, of the previous and current month. This leads to the following optimization formulation: Minimize: Subject to: Optimization results are illustrated in Table 2. The optimal production and inventory levels in low moisture content months are greater than the high moisture months of May and June. Production cost results of annual optimization can be reduced by $0.44M when compared to operation at a constant production rate of 48,000 admt, a 0.3% reduction in cost.

RTO with Moisture Content Variation
The RTO uses the updated actual moisture content of wood that is considered to vary in 3-day subperiods within a month. This variation affects pulp production rates and costs. Inventory constraints are imposed using a rolling horizon approach with actual moisture content variations in which inventory and production rate results from optimization are updated and used in next period.
The objective function minimizes production cost in 10 periods over a given month and includes a penalty term to reduce errors from inventory changing with moisture content variation where coefficient (α) of this term is weighting inventory cost against production cost; in this study, it was specified as 0.1%. This percentage value is a tuning parameter of the algorithm, but its value reflects the relative size of the inventory compared to the production rate and the cost of storage relative to the cost of production. Constraints of the RTO are different from annual optimization, inventory constraints with moisture content variation are added, based on a rolling horizon formulation; the details are as follow: Minimize: where ∆I m is the deviation of the inventory rate in the period of moisture content variation, N is number of periods that vary moisture content, and I 1 is inventory rate at current period for optimization. The optimization problem consists of 10 periods where the first period has the estimated moisture content of the wood for that period which, for purposes of the results in this paper, is drawn from a Gaussian distribution with a mean equal to the monthly moisture content and with a standard deviation of 5% of the mean. The later periods' moisture contents are fixed at their average MC value based on the month in which they fall. In the simulation below, the month start optimization of this work is June because the May inventory is neither maximum nor minimum level and also there is a different moisture content of wood between these months. The results for other months are similar and shown in Appendix A. RTO results for the month of June are represented in Table 3. Although low moisture content periods have higher production rates in general, sometimes the production rate is less than high moisture content sub-periods, such as period 2 and 10. In overall production cost, RTO is slightly more effective than the optimized annual case since it has costs lower by $10,000. Reduction of production cost is not always guaranteed every month due to some months having low inventory rate from the previous period and, therefore, must have an increased production rate to remain feasible. Therefore, performance of the RTO should be considered on an annual basis as shown in Table 4.  Tables 3 and 4 demonstrate that the RTO reduces the monthly and annual cost relative to the optimized monthly and annual approach, but the gain is small and could be due to random variations in the 3-day period moisture content relative to the average for this particular sample of the random variable.

Confirming the Importance of Moisture Content of Wood Variation
To confirm that the moisture content of wood variation is a significant factor for RTO, RTO results over a 20-year horizon were generated in which the annual production cost found from the rolling horizon RTO is compared with results of the annual optimization case. The comparison results are shown in Figure 5. This demonstrates that, for every year, the production costs of RTO are smaller.  Tables 3 and 4 demonstrate that the RTO reduces the monthly and annual cost relative to the optimized monthly and annual approach, but the gain is small and could be due to random variations in the 3-day period moisture content relative to the average for this particular sample of the random variable.

Confirming the Importance of Moisture Content of Wood Variation
To confirm that the moisture content of wood variation is a significant factor for RTO, RTO results over a 20-year horizon were generated in which the annual production cost found from the rolling horizon RTO is compared with results of the annual optimization case. The comparison results are shown in Figure 5. This demonstrates that, for every year, the production costs of RTO are smaller. To provide more quantitative statistical confirmation, we propose two hypotheses:

Hypothesis 2 (H2). Mean of annual production cost of RTO is less than annual algorithm
where H1 and H2 refer to null and alternative hypotheses. The t-test hypothesis result in Table 5 can identify that the mean of production cost of RTO is less than annual algorithm with 95% confidence because the t-stat is not inside the region of the null hypothesis and the p-value is less than 0.05.

Conclusions
In this work, a surrogate model was developed for pulp mill modeling and optimization in which the production cost variation with the wood moisture content and production rate is empirically determined using data from an IDEAS simulation of pulping operations. The model was used for optimization to minimize production cost and determine optimal production rate at different time scales.
The monthly variation in wood moisture content arises from seasonal patterns of rainfall and growth and can be predicted on average. The production rate is increased during low moisture content periods to reduce the production cost in high moisture content periods with inventory used to buffer the system operation. The annual production cost can reduce to 0.3% when compared with case without inventory balance.
For real-time optimization, production is scheduled using a rolling horizon approach where the actual wood moisture content varies every 3 days during the month. RTO with moisture content variation reduces the annual production cost each year over a twenty-year horizon compared to an annual optimization approach. The results support the conclusion that the moisture content of wood can be a significant factor in the economic operation of pulp mills in Thailand.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A.
Appendix A. 1

. Pulp Process Model
To generate the black-box model used in this paper, a polynomial superstructure and least squares fitting was used along with the Aike Information Criterion (AIC) method to determine the final structure of model that balanced the accuracy of the fitting with the model complexity as measured by the number of terms in the polynomial. In the polynomial superstructure, all input variables of each operation were considered to predict output variable. Models were constructed by retaining subsets of terms of variables and combinations of variables in quadratic or ternary form. A p-value (less than 0.05) criterion was used to assess whether the model terms were considered a fit, where the least square errors of the model results from the simulation (IDEAS) results were used as the objective function of the fitting procedure. The minimum AIC was used to select the best model. The optimal model includes significant terms, lowest AIC, and percentage error with simulated result comparison. An example of selecting the best model of daily operating cost is shown in Table A1. As result, model 3 is the optimal model to predict daily operating cost (C total ) as shown in Equation (A1), the same as Equation (3)  The result of RTO with moisture content variation in one year, as shown in Table A2.