Model Predictive Control and Optimization for Papermaking Processes

1. Introduction

Papermaking is a large-scale two-dimensional process. It has to be monitored and controlled continuously in order to ensure that the qualities of paper products stay within their specifications. There are two types of control problems involved in papermaking processes: machine directional (MD) control and cross directional (CD) control. Machine direction refers to the direction in which paper sheet travels and cross direction refers to the direction perpendicular to machine direction. The objectives of MD control and CD control are to minimize the variation of the sheet quality measurements in machine direction and cross direction, respectively. This chapter considers the design and applications of model predictive control (MPC) for papermaking MD and CD processes.

MPC, also known as moving horizon control (MHC), originated in the late seventies and has developed considerably in the past two decades (Bemporad and Morari 2004; Froisy 1994; Garcia et al. 1998; Morari & Lee 1999; Rawlings 1999; Chu 2006). It can explicitly incorporate the process’ physical constraints in the controller design and formulate the controller design problem into an optimization problem. MPC has become the most widely accepted advanced control scheme in industries. There are over 3000 commercial MPC implementations in different areas, including petro-chemicals, food processing, automotives, aerospace, and pulp and paper (Qin and Badgwell 2000; Qin and Badgwell 2003).

Honeywell introduced MPC for MD controls in 1994; this is likely the first time MPC technology was applied to MD controls (Backström and Baker, 2008). Increasingly, paper producers are adopting MPC as a standard approach for advanced MD controls.

MD control of paper machines requires regulation of a number of quality variables, such as paper dry weight, moisture, ash content, caliper, etc. All of these variables may be coupled to the process manipulated variables (MV’s), including thick stock flow, steam section pressures, filler flow, machine speed, and disturbance variables (DV’s) such as slice lip adjustments, thick stock consistency, broke recycle, and others. Paper machine MD control is truly a multivariable control problem.

In addition to regulation of the quality variables during normal operation, a modern advanced control system for a paper machine may be expected to provide dynamic economic optimization on the machine to reduce energy costs and eliminate waste of raw materials. For machines that produce more than one grade of paper, it is desired to have an automatic grade change feature that will create and track controlled variable (CV) and MV trajectories to quickly and safely transfer production from one grade to the next. Basic MD-MPC, economic optimization, and automatic grade change are discussed in this chapter.

MPC for CD control was introduced by Honeywell in 2001 (Backström et al. 2001). Today, MPC has become the trend of advanced CD control applications. Somesuccessful MPC applications for CD control have been reported in (Backström et al. 2001, Backström et al. 2002; Chu 2010a; Gheorghe 2009).

In papermaking processes, it is desired to control the CD profile of quality variables such as dry weight, moisture, thickness, etc. These properties are measured by scanning sensors that traverse back and forth across the paper sheet, taking as many as 2000 or more samples per sheet property across the machine. There may be several scanners installed at different points along the paper machine and so there may be multiple CD profiles for each quality variable.

The CD profiles are controlled using a number of CD actuator arrays. These arrays span the paper machine width and may contain up to 300 individual actuators. Common CD actuators arrays allow for local adjustment, across the machine, of: slice lip opening, headbox dilution, rewet water sprays, and induction heating of the rolls. As with the CD measurements, there may be multiple CD actuator arrays of each type available for control. By changing the setpoints of the individual CD actuators within an array, one can adjust the local profile of the CD measurements.

The CD process is a multiple-input-multiple-output (MIMO) system. It shows strong input and output off-diagonal coupling properties. One CD actuator array can have impact on multiple downstream CD measurement profiles. Conversely, one CD measurement profile can be affected by multiple upstream CD actuator arrays. Therefore, the CD control problem consists of attempting to minimize the variation of multiple CD measurement profiles by simultaneously optimizing the setpoints of all individual CD actuators (Duncan 1989).

MPC is a natural choice for paper machine CD control because it can systematically handle the coupling between multiple actuator and multiple measurement arrays, and also incorporate actuator physical constraints into the controller design. However, different from standard MPC problems, the most challenging part of the cross directional MPC (CD-MPC) is the size of the problem. The CD-MPC problem can involve up to 600 MVs, 6000 CVs, and 3000 hard constraints. Also, the new setpoints of MVs are required as often as every 10 to 20 seconds. This chapter discusses the details of the design for an efficient large-scale CD-MPC controller.

This chapter has 5 sections. Section 2 provides an overview of the papermaking process highlighting both the MD and CD aspects. Section 3 focuses on modelling, control and optimization for MD processes. Section 4 focuses on modelling, control and optimization for CD processes. Both Sections 3 and 4 give industrial examples of MPC applications. Finally, Section 5 draws conclusions and provides some perspective on the future of MD-MPC and CD-MPC.

2. Overview of papermaking processes

A flat sheet of paper is a network consisting of cellulose fibres bound to one another. A paper machine transforms a slurry of water and wood cellulose fibres into this type of network. The whole papermaking process can be regarded as a water-removal system: the consistency of fibre solutions, called stock by papermakers, increases from around 1% at the beginning of a paper machine (the headbox) to around 95% at the end (the reel).

2.1. Brief description of papermaking processes

In general a paper machine can be divided into four sections: forming section, press section, drying section, and calendering section. In the forming section, the stock flow enters the headbox to be distributed evenly across a continuously running fabric felt called the wire. The newly formed sheet is carried by the wire along the Fourdrinier table, which has a set of drainage elements that promote water removal by various gravity and suction mechanisms. These elements include suction boxes, couch rolls, foils, etc. The solid consistency of the paper web can reach 20% by the time the web leaves the forming section and enters the press section. Figure 1 illustrates the configuration of a Fourdrinier-type paper machine.

Figure 1.

The configuration of a Fourdrinier-type paper machine

The press section may be considered as an extension of the water-removal process that was started on the wire in the forming section. Typically, it consists of 1 – 3 rolling press nips. When the paper web passes through these nips, the pressing roll squeezes water out and consolidates the web formation at the same time. In the press section, both the surface smoothness and the web strength are improved. As higher web strength is achieved in the press section, better runability will be observed in the drying section. A paper machine is typically operated at a very high speed. The fastest machine speed may be as high as 2,200 meters per minute.

The drying section includes multiple drying cylinders which are heated by high temperature and high pressure steam. The heat is transferred from steam onto the paper surface through these rotating steel cylinders. The heat flow increases the paper surface temperature to the point where water starts evaporating and escaping from the paper web. The drying section is the most energy consuming part of paper manufacturing. Before the paper enters the drying section, the solid consistency is around 50%. After the drying section, the consistency can reach 95%, which corresponds to a finished product moisture specification.

The last section of the paper machine is called the calendering section. Calendering is a terminology referring to pressing with a roll. The surface and the interior properties of the paper web are modified when it passes through one or more calendering nips. Typically the calendering nip consists of one or multiple soft/hard or hard/hard roll pairs. The hard roll presses the paper web against the other roll, and deforms the paper web plastically. By this means, the calender roll surface is replicated onto the paper web. Depending on the type of paper being produced, the primary objective of calendering may be to produce a smooth paper surface (for printing), or to improve the uniformity of CD properties, such as paper caliper (thickness).

More details of paper machine design and operation are given in (Smook 2002; Gavelin 1998).

2.2. Paper qualitymeasurement

A paper machine can have one or more measurement scanners. The quality measurement sensors are mounted on the scanner head which travels back and forth across the paper web to provide online quality measurements. The most common paper machine quality measurements include dry weight, moisture, and caliper. Dry weight indicates the solid weight per unit area of a sheet of paper. For different types of products, the value of dry weight can vary from 10 grams per square meter (gsm), in the case of paper tissue, to 400 gsm, in the case of heavy paper board. Moisture content is another critical quality property of the finished paper product. It indicates the mass percentage of water contained in a sheet of paper. Moisture content is a key factor determining the strength of the finished product. Typical moisture targets range from 5% to 9%. Caliper is the measure of the thickness of a sheet of paper. It is a key factor determining the gloss and printability of the finished product. The caliper targets are in the range from 70$\text{μm}$ to 300 $\text{μm}$ depending on the production grade. In general the online measurements for dry weight, moisture and caliper are available and used for both the MD and CD feedback controller designs.

As the scanners travel across a moving sheet, the real data collected actually comes from a zig-zag trajectory (See Figure 2). These data contain both CD and MD variation. A reliable MD/CD separation scheme is the prerequisite for MD and CD controldesigns. Since the MD/CD separation is a separate topic, the rest of this chapter assumes that the pure MD/CD measurements have been obtained prior to the MD/CD controller development. The scanner measurements are denoted by x(i, t),$\text{i}=1,\cdots ,\text{n}$ indexes the n measurements taken across the sheet each scan (CD measurement index), and t is the time stamp of each scan (MD measurement index).$\overline{x}(\text{t})$ is the MD measurement given by

$$\overline{x}\left(t\right)=\frac{1}{n}{\displaystyle {\sum }_{i=1}^{n}x\left(i,t\right)}$$E1

Figure 2.

The zig-zag scanner trajectories

2.3. Brief description of MD control

The objective of MD control is to minimize the variation of the sheet quality measurements in machine direction.

A number of actuators are available for control of the MD variables. Stock flow to the headbox is regulated by the stock flow valve or variable speed pump. As stock flow increases, the amount of fibre flowing into the forming section increases and dry weight and caliper increase. At the same time, there is also more water coming through the machine and moisture will increase. So, changes in the stock flow affect dry weight, moisture and caliper. The steam pressure in the cylinders of the drying section may be adjusted. As the steam pressure in the cylinders increases, so does the temperature in the cylinders, and more heat is transferred to the paper. In this way, steam pressure affects moisture. Typically, the dryer cylinders are divided into groups, and the steam pressure for each of these dryer sections may be adjusted independently. Machine speed affects dry weight and caliper, as increasing the machine speed stretches the paper web thinner, giving it less mass per unit area. Machine speed also affects moisture as both the drying properties of the paper and the residence time in the dryer change. Clearly, paper machine MD control is a multivariable control problem.

2.4. Brief description of CD control

The objective of CD control is to achieve uniform paper qualities in the cross direction, i.e., to minimize the variation of CD profiles. The CD variation, can be formulated as two times of the standard deviation of the CD profile,

$$2{\sigma }_{CD}\left(t\right)=2^{\ast }{\left(\frac{1}{n-1}{\displaystyle {\sum }_{i=1}^n{\left(x\left(i,t\right)-\overline{x}\left(t\right)\right)}^2}\right)}^{\frac{1}{2}}$$E2

Often the term ‘CD spread’ is used interchangeably with 2σ_CD.

CD actuators are used to regulate CD profiles and improve the uniformity of paper quality properties in the cross directioni.e., reduce the value of 2σ_CD.

The most common dry weight CD actuators are both located at the headbox.The headbox slice opening is a full-width orifice or nozzle that can be adjusted at points across the width of the paper machine. This allowsfor differences in the local stock flow onto the wire across the machine. The consistency profiler changes the consistency of local stock flow by injecting dilution water and altering the local concentration of pulp fibre across the headbox. Headbox slice and consistency profiler are primarily designed for dry weight control, but they have the effects on both moisture and caliper profiles. Figure 1 indicates the location of headbox dry weight actuators.

The most common moisture actuators are the steam box and water spray. The steam box applies high temperature steam to the surface of the moving paper web. As the latent heat in the steam is released and heats up the paper web, it lowers the web viscosity and eases dewatering in the press section. The water spray regulates the moisture profiles according to a different mechanism. It deploys a fine water spray to the paper surface through a set of nozzles across the machine width to re-moisturize the paper web. Similar to the dry weight actuators, moisture actuators are designed for moisture profile regulation but they may have effects on the caliper profile. The steam box is typically installed in the press section and the water spray is located in the drying section. Figure 1 indicates the physical locations of moisture actuators.

The most common caliper actuators are hot air showers and induction heaters. Both types of actuators provide surface heating for calendering rolls. The hot shower uses the high temperature steam or air; the induction heater uses the high frequency alternating current. By heating up the calender roll, caliper actuators alter the local diameter of the calender roll and subsequently increase the local pressure applied to the paper web. The physical location of the caliper actuators can be also found in Figure 1.

3. Modelling, control and optimization of papermaking MD processes

Control of the MD process is typically a regulation problem where the paper quality variables need to be held within specified quality limits. At the same time, the constant pressure to increase operational efficiency demands that the paper production uses the least amount of raw materials and energy required to meet quality goals.

When multiple grades of paper are produced on a single paper machine, control must also be able to provide quick transitions of the quality variables along smooth trajectories. When the differences between the paper grades are large, it may be necessary to enhance the control algorithm to account for process nonlinearities that become apparent over a larger span of operating points.

3.1. Modeling of papermaking MD processes

In this section, modelling of the MD process for MPC controller design is discussed. The additional modelling required for paper grade change control is discussed in section 3.4.

For effective MPC control of paper MD quality variables, it is necessary to build a matrix of linear models relating the process MV’s to the quality variables (the CV’s). A basic paper machine model matrix most often includes stock flow, steam pressure of multiple dryer sections, and machine speed as MVs, and paper weight (basis weight or dry weight), and moisture as CV’s. Many other MV’s and CV’s can be included in the model matrix depending on the complexity of the paper machine and the paper quality requirements. An example model matrix is given in Figure 3.

Figure 3.

A basic model matrix for CD-MPC. The models are step responses.

The process models used for MPC control are often developed in transfer function form, such as:

$$y\left(s\right)=g\left(s\right)u\left(s\right)+d\left(s\right)$$E3

And

$$y\left(s\right)=\left[\begin{array}{c}{y}_1\left(s\right)\\ ⋮\\ y_{N_y}\left(s\right)\end{array}\right],g\left(s\right)=\left[\begin{array}{c}{g}_{11}\left(s\right)\\ ⋮\\ g_{N_{y}1}\left(s\right)\end{array}\begin{array}{c}\cdots \\ \ddots \\ \cdots \end{array}\begin{array}{c}{g}_{1N_u}\left(s\right)\\ ⋮\\ g_{N_y{N}_u}\left(s\right)\end{array}\right],u\left(s\right)=\left[\begin{array}{c}{u}_1\left(s\right)\\ ⋮\\ u_{N_u}\left(s\right)\end{array}\right]$$E4

where$\text{y}(\text{s})\in {ℂ}^{{\text{N}}_{\text{y}}}$ is the Laplace transformation of the${\text{N}}_{\text{y}}$ MD quality measurements (such as dry weight, moisture, caliper, etc).$\text{u}(\text{s})\in {ℂ}^{{\text{N}}_{\text{u}}}$ is the Laplace transformation of the actuator setpoints (such as thick stock flow, dryer section steam pressure(s), filler flows, machine speed, etc). d(s)$\in {ℂ}^{{\text{N}}_{\text{y}}}$ is the Laplace transformation of the augmented process disturbance array. ${\text{g}}_{\text{ij}}(\text{s})\in ℂ$($\text{i}=1\dots {\text{N}}_{\text{y}}\text{andj}=1\dots {\text{N}}_{\text{u}})$are the transfer functions from the j^th actuator ${\text{u}}_{\text{j}}$ to the i^th MD quality measurement${\text{y}}_{\text{i}}$.

3.1.1. Model identification

Typically the MD process models that are used as the basis for the MD-MPC controller are identified from data obtained during simple process experiments. A series of step changes is made for each MV. There is a delay after each step long enough so that the full responses of all of CV’s can be observed. That is, the CV responses reach steady state before the next step change in the MV is made. These types of process experiments are known as bump tests.

Once a set of identification data has been obtained, various techniques may be used to generate a process model from this data. In the simplest case, plots of the bump tests are reviewed to graphically estimate a process gain, $k_{ij}$, dead time,$T_{d_{ij}}$ , and time constant, $T_{p_{ij}}$, yielding the process model:

$$y_i\left(s\right)=\frac{k_{ij}{e}^{-T_{d_{ij}}s}}{1+T_{p_{ij}}s}{u}_j\left(s\right)$$E5

More complex methods involve use of regression and search techniques to find both the optimum model parameters, and the optimum model structure (transfer function numerator and denominator orders). These techniques typically use minimization of squared model prediction errors as the objective:

$$J={\displaystyle \sum }_{k=1}^N{\left(\widehat{y_i}\left(k\right)-y_i\left(k\right)\right)}^2$$E6

Where $\widehat{y_i}\left(k\right)$ and$y_i\left(k\right)$ are respectively the predicted and actual values of the i^th CV at time k. (Ljung 1998) is the classic reference on system identification, and there are commercial software packages available that automate much of the system identification work.

3.2. MD-MPC design

Once all of the bump test, and system identification activities have been performed, the complete process model (3) is used directly in the model predictive controller. MPC solves an optimization problem at each control execution. One robust MPC problem formulation is:

$${\min }_{\text{Δu},{\text{y}}_{\text{r}}}\frac{1}{2}\text{W}{\left({\text{y}}_{\text{r}}-\text{SΔu}\right)}_2^2$$E7

Subject to:

$$y_l\le y_{\text{r}}\le y_{\text{h}}$$E8

$$u_{\text{l}}\le \text{u}\le {\text{u}}_h$$E9

$$\Delta u_l\le \Delta \text{u}\le \Delta {\text{u}}_h$$E10

The values of $y_r$ are the CV targets and SΔU are the predicted future values of the CV’s. S is the prediction matrix, containing all the information from the process model (3). W is a weighting matrix, and ${\Vert ·\Vert }_2^2$ is the two-norm squared operator.$y_l$ and $y_{\text{h}}$ are the low and high CV quality limits, $u_{\text{l}}$and ${\text{u}}_h$ are the low and high MV limits, and $\Delta u_l$ and $\Delta {\text{u}}_h$ are the low and high limits for MV moves. As discussed in the section below, this problem formulation, combined with techniques employed in its solution implicitly provide robustness characteristics in the controller design.

The technical details of the solution of the problem (7) are given in (Ward 1996); however, some notable aspects of the solution methodology and beneficial characteristics of the solution are given in the section below.

3.2.1. Model scaling and controller robustness

Robust numerical solution of the optimization problem (7) depends on condition number of the system gain matrix, G. Prior to performing the controller design, the gain matrix condition number is minimized by solving the problem:

$${\min }_{{\text{D}}_{\text{r},}{\text{D}}_{\text{c}}}\left\{\text{κ}\left({\text{D}}_{\text{r}}{\text{gD}}_{\text{c}}\right)\right\}$$E11

Where$\text{κ}$ is condition number, and D_r and D_c are diagonal transformation matrices. The scaled system gain matrix g_s is then:

$${\text{g}}_{\text{s}}={\text{D}}_{\text{r}}{\text{gD}}_{\text{c}}$$E12

${\text{g}}_{\text{s}}$is then used for all MPC computations.

It should be noted that the objective (7) does not explicitly penalize the MV moves Δu as a method to promote controller robustness. Instead, controller robustness is provided by the CV range formulation and singular value thresholding.

First, the CV range formulation refers to the inequalities given in the problem formulation (7). Under this formulation, if a CV is predicted to be within it range in the future, no MV action is taken. Since MV moves are not made unless absolutely necessary, this is a very robust policy.

Second, the solution of the problem involves an active set method that allows the constrained optimization problem to be converted into an unconstrained problem. A URV orthogonal decomposition (see Ward 1996) of the matrix characterizing the unconstrained problem is then employed to solve the unconstrained problem. Prior to the decomposition, singular values of the problem matrix that are less than a certain threshold are dropped, reducing the dimension of the problem, and ensuring that the controller does not attempt to control weakly controllable directions of the process.

3.3. Economic optimization

Energy consumption is a big concern for papermakers. Increasing profits by minimizing operating costs without sacrificing paper quality and runability is always a goal for them. In theory, if the number of MVs of a process is greater than the number of CVs plus the number of active constraints, the process has degrees of freedom allowing for steady-state optimization. Product value optimization can be systematically integrated with the MD-MPC control. One can then take the feed, product, and utility costs into account with the MD controller design.

For the economic optimization of a process, the following objective is to be minimized:

$$\begin{array}{l}J={\displaystyle \sum }_i\left(a_{y_i}{\left(y_i-y_{i0}\right)}^2+b_{y_i}{y}_i\right)\\ +{\displaystyle \sum }_j\left(a_{u_j}{\left(u_j-u_{j0}\right)}^2+b_{u_j}{u}_j\right)\end{array}$$E13

Here $y_{i0}$ and $u_{j0}$ are the desired steady state values of the process CV’s and MV’s, $a_{y_i}$and $a_{u_j}$ are the costs of quadratic deviation from the desired values, and $b_{y_i}$ and $b_{u_j}$ are the linear costs of the CV’s and MV’s. This objective is useful for paper machines, for example, by placing costs on the different energy sources used in drying.

The economic objective is combined with the MPC control problem objective to give an augmented problem formulation. The augmented problem is then solved using the same solution method as described above.

Economic optimization is a lower priority for paper machines than quality control. If the paper does not meet quality specifications, it cannot be sold, and any savings made from economic optimization are more than lost. Therefore, economic optimization only occurs when there are extra degrees of freedom for the controller. Economic optimization is not attempted unless all of the CV’s are predicted to remain within their quality specifications over the whole of the controller’s prediction horizon.

3.3.1. Mill implementation results

MPC including an economic optimization layer was implemented for a tissue machine. A diagram of the tissue machine is given in Figure 4. As can be seen in the diagram, tissue dry weight and moisture are measured at the reel; moisture is measured between the second through-air dryer (TAD2) and the Yankee dryer, and TAD1 exhaust pressure must also be monitored and controlled. These four variables are the CV’s in this example. A large number of MV’s are available to control this machine. Stock flow, TAD1 supply temperature, TAD1 dry end differential pressure, TAD1 gap pressure, TAD2 exhaust temperature, TAD2 dry end differential pressure, TAD2 gap pressure, Yankee hood temperature, and Yankee supply fan speed are all used as MV’s in the MPC. Machine speed and tickler refiner were added as DV’s. The MPC model matrix is shown in Figure 5.

MV	Energy Fuel	Units	Linear Obj. Coef.; Cost/eng unit
TAD1 Supply Temp	Gas	Deg F	0.680
TAD1 DE DP	Electricity	Inch H2O	47.267
TAD1 Gap Pres	Electricity	Inch H2O	-0.030
TAD2 Exh Temp	Gas	deg F	5.858
TAD2 DE DP	Electricity	Inch H2O	40.249
TAD1 Gap Pres	Electricity	Inch H2O	-16.415

Table 1.

Tissue machine MV’s with linear objective coefficients

A large number of the MV’s in this control problem have an impact on the paper moisture, both after TAD2, and at the reel; however each MV uses a different energy source and has different drying efficiency. Overall,since there are more MV’s than CV’s, and significant cost differences between the MV’s, there is an opportunity for economic optimization in this system. Table 1 shows the energy sources and different energy cost efficiencies (Linear Obj Coeff Cost/eng unit) associated with each MV. Economic optimization can be accomplished by including these variables in the linear part of the economic objective function given by (10).

Figure 4.

Diagram of a tissue machine with CV’s, MV’s, and DV’s for MPC.

Once the economic cost function was added to the MPC, a plant trial was made. Figures 6-10 show the results of this trial. In Figure 7 it can be seen that prior to turning on the economic optimizer (the period from 8:30 to 9:30) there was a relative cost of energy of 100. The optimizer was turned on at 9:30. Initially there were some wind-up problems in the plant DCS which were preventing the MPC from optimizing. Once these were cleared, at 10:44, the controller drove the process to the low cost operating point (from 10:44 to 12:30). The relative cost of energy at this operating point was 98.8. In order to better interpret these results, it is necessary to rank the costs of each MV on the common basis of Cost/% Moi. This is accomplished by dividing the linear objective coefficients given in Table 1 by their respective process gains. These are shown in Table 2, along with the MV high and low limits, and the optimization behaviour. Looking again at the Figures 8 and 9, it can be seen that the highest costing MV’s are driven to their minimum operating points, and the lowest costing MV’s are driven to their maximum operating points. The TAD1 dry end differential pressure is left as the MV that is within limits and actively controlling the paper moistures. Figure 10 shows that throughout this trial, the MV’s are optimized without causing any disturbance to the CV’s.

Figure 5.

The MPC model matrix for the tissue machine control and optimization example

Figure 6.

Natural gas costs and electricity costs during the trial

Figure 7.

Total costs during the trial

Table 2.

The MV cost rankings.

Figure 8.

Manipulated variables during the trial

Figure 9.

Manipulated variables during the trial

Figure 10.

Controlled variables during the trial

3.4. Grade change strategies

Grade change is a terminology in MD control. It refers to the process of transitioning a paper machine from producing one grade of paper product to another. One can achieve a grade change by gradually ramping up a set of MVs to drive the setpoints of CVs from one operating point to another. During a grade change, the paper product is often off-specification and not sellable. It is important to develop an automatic control scheme to coordinate the MV trajectories and minimize the grade change transition times and the off-spec product. An offline model predictive controller can be designed to produce CV and MV trajectories to meet these grade change criteria. MPC is well-suited to this problem because it explicitly considers MV and CV trajectories over a finite horizon. By coordinating the offline grade change controller (linear or nonlinear) and an online MD-MPC, one can derive a fast grade change that minimizes off-spec production. This section discusses the design of MPC controllers for linear and nonlinear grade changes.

Figure 11 gives a block diagram of the grade change controller incorporated into an MD control system. The grade change controller calculates the MV and CV trajectories to meet the grade change criteria. This occurs as a separate MPC calculation performed offline so that grade change specific process models can be used, and so that the MPC weightings can be adjusted until the MV and CV trajectories meet the design criteria. The MV trajectories are sent to the regulatory loop as a series of MV setpoint changes. The CV trajectories are sent as setpointchanges to the MD controller. If the grade change is performed with the MD controller in closed-loop, additional corrections to the MV setpoints are made to eliminate any deviation of the CV from its target trajectory.

Figure 11.

Block diagram of MD-MPC control enhanced with grade change capability

The MV and CV trajectories are generated in a two step procedure. First there is a target calculation step that generates the MV setpoints required to bring the CV’s to their target values for the new grade. Once the MV setpoints are generated, then there is a trajectory generation step where the MV and CV trajectories are designed to meet the specifications of the grade change.

The MV targets are generated from solving a set of nonlinear equations:

$$\begin{array}{l}{\text{y}}_{\text{dw}}^1-{\text{f}}_{\text{dw}}^1\left({\text{u}}_1,{\text{u}}_2,{\text{u}}_3,\dots ,{\text{C}}_1,{\text{C}}_2,{\text{C}}_3,\dots \right)=0,\\ {\text{y}}_{\text{dw}}^2-{\text{f}}_{\text{dw}}^2\left({\text{u}}_1,{\text{u}}_2,{\text{u}}_3,\dots ,{\text{C}}_1,{\text{C}}_2,{\text{C}}_3,\dots \right)=0,\\ ⋮\\ {\text{y}}_{\text{moi}}^1-{\text{f}}_{\text{moi}}^1\left({\text{u}}_1,{\text{u}}_2,{\text{u}}_3,\dots ,{\text{C}}_1,{\text{C}}_2,{\text{C}}_3,\dots \right)=0,\\ {\text{y}}_{\text{moi}}^2-{\text{f}}_{\text{moi}}^2\left({\text{u}}_1,{\text{u}}_2,{\text{u}}_3,\dots ,{\text{C}}_1,{\text{C}}_2,{\text{C}}_3,\dots \right)=0,\\ ⋮\end{array}$$E14

Here y_dw/y_moi represents the CV target for the new grade. The functions $\text{f}(\cdot )$ are the models of dry weight and moisture. The process MV’s are denoted u_i and model constants are denoted C_i. The superscripts indicate the same paper properties measured by different scanners. Since the number of MV’s and the number of CV’s is not necessarily equal, these equations may have one, multiple or no solutions. To allow for all of these cases, the problem is recast as:

$${\min }_{{\text{u}}_{\text{i}}}\text{F}\left({\text{u}}_1,{\text{u}}_2,\dots \right)$$E15

Subject to:

$$\text{G}\left({\text{u}}_1,{\text{u}}_2,\dots \right)\le 0,$$E16

$$\text{H}\left({\text{u}}_1,{\text{u}}_2,\dots \right)=0,$$E17

Where $\text{F}(\cdot )$ is a quadratic objective function formulated to find the minimum travel solution. $\text{H}(\cdot )$represents the equality constraints given above, and $\text{G}(\cdot )$ represents the physical limitations of the CVs and MVs (high, low, and rate of change limits).

Once the MV targets have been generated, the MV and CV trajectories are then designed. Figure 12 gives a schematic representation of the trajectory generation algorithm. The process models are linearized (if necessary) and then scaled and normalized for application in an MPC controller. Process constraints such as the MV and CV targets, and the MV high and low limits are also given to the MPC controller. Internal controller tuning parameters are then used to adjust the MV and CV trajectories to meet the grade change requirements.

Figure 12.

Diagram of MPC-based grade change trajectory generation.

3.4.1. Linear grade change

In a linear grade change, the MD process models that are used in the MD-MPC controller are alsoused as the models for determining the MD targets, and for designing the MD grade change trajectories.

3.4.2. Nonlinear grade change

In a nonlinear grade change, a first principles model may be used for the target and trajectory generation. For example, a simple dry weight model is:

$${\text{m}}_{\text{dry}}=\text{K}\frac{{\text{q}}_{\text{stock}}}{\text{v}}$$E18

Where ${\text{m}}_{\text{dry}}$ is the paper dry weight, ${\text{q}}_{\text{stock}}$is the thick stock flow, and v is machine speed. K is the expression of a number of process constants and values including fibre retention, consistency, and fibre density. (Chu et al. 2008) gives a more detailed treatment of this dry weight model.

(Persson 1998, Slätteke 2006, and Wilhelmsson 1995) are examples of first principles moisture models that may be used.

3.4.3. Mill implementation results

In this section, some results of MPC-based grade changes for a fine paper machine are given. The grade change is from a paper with a dry weight of 53 lb/3000ft²(86 g/m²) to a paper with a dry weight of 44 lb/3000ft²(72 g/m²). Both paper grades have the same reel moisture setpoint of 4.8%. For the grade change, stock flow, 6^th section steam pressure, and machine speed are manipulated.

Figures 13 and 14 show a grade change performed on the paper machine using linear process models, and keeping the regular MPC in closed-loop during the grade change. The grade change was completed in 10 minutes, which is a significant improvement over the 22 minutes required by the grade change package of the plant’s previous control system. In Figure 13, the CV trajectories are shown. Here it can be seen that although there is initially a small gap between the actual dry weight and the planned trajectory, the regular MPC takes action with the thick stock valve (as shown in Figure 14) to quickly bring dry weight back on target. The deviation in the reel moisture is more obvious. This might be expected as the moisture dynamics of the paper machine display more nonlinear behaviour for this range of operations. The steam trajectory in Figure 14 is ramping up at its maximum rate and yet the paper still becomes too wet during the initial part of the grade change. This indicates that the grade change package is aggressively pushing the system to achieve short grade change times.

Figure 13.

CV trajectories under closed-loop GC with linear models

Figure 14.

MV trajectories under closed-loop GC with linear models

Figures 15 and 16 show a grade change performed on a high fidelity simulation of the fine paper machine. This grade change uses a nonlinear process model, and the regular MPC is kept in closed-loop during the grade change. Here it can be seen that the duration of the grade change is reduced to 8 minutes. Part of the improvement comes from using stock flow setpoint instead of stock valve position, allowing improved dry weight control. Another improvement is that the planned trajectories allow for some deviation in the reel moisture that cannot be eliminated. Both dry weight and reel moisture follow their trajectories more closely. At the end of the grade change, the nonlinear grade change package is able to anticipate the need to reduce steam preventing the sheet from becoming dry.

Figure 15.

CV trajectories under closed-loop GC with nonlinear models

Figure 16.

MV trajectories under closed-loop GC with nonlinear models

4. Modelling, control and optimization of papermaking CD processes

To produce quality paper it is not enough that the average value of paper weight, moisture, caliper, etc across the width of the sheet remains on target. Paper properties must be uniform across the sheet. This is the purpose of CD control.

4.1. Modelling of papermaking CD processes

The papermaking CD process is a large scaled two-dimensional process. It involves multiple actuator arrays and multiple quality measurement arrays. The process shows very strong input-output off-diagonal coupling properties. An accurate CD model is the prerequisite for an effective CD-MPC controller. We begin by discussing a model structure for the CD process and the details of the model identification.

4.1.1. Atwo-dimensional linear system

The CD process can be modelled as a linear multiple actuator arrays and multiple measurement arrays system,

$$Y(s)=G(s)U(s)+D(s)$$E19

and

$$\begin{array}{l}Y(s)=\left[\begin{array}{c}{\text{y}}_1(\text{s})\\ ⋮\\ {\text{y}}_{{\text{N}}_{\text{y}}}(\text{s})\end{array}\right],\text{G}(\text{s})=\left[\begin{array}{ccc}{\text{G}}_{11}(\text{s})& \dots & {\text{G}}_{1{\text{N}}_{\text{u}}}(\text{s})\\ ⋮& \ddots & ⋮\\ {\text{G}}_{{\text{N}}_{\text{y}}1}(\text{s})& \dots & {\text{G}}_{{\text{N}}_{\text{y}}{\text{N}}_{\text{u}}}(\text{s})\end{array}\right],\\ U(s)=\left[\begin{array}{c}{u}_1(\text{s})\\ ⋮\\ u_{{\text{N}}_u}(\text{s})\end{array}\right]\end{array}$$E20

where $\text{Y}(\text{s})\in {ℂ}^{({\text{N}}_{\text{y}}\cdot \text{m})}$ is the Laplace transformation of the augmented CD measurement array. The element ${\text{y}}_{\text{i}}(\text{s})\in {ℂ}^{\text{m}}(\text{i}=1,\dots ,{\text{N}}_{\text{y}})$ is the Laplace transformation of the i^th individual CD measurement profile, such as dry weight, moisture, or caliper. ${\text{N}}_{\text{y}}$is the total number of the quality measurements, and m is the number of elements of individual measurement arrays. $$\text{U}(\text{s})\in {ℂ}^{({\displaystyle \sum }_{\text{j}=1}^{{\text{N}}_{\text{u}}}{\text{n}}_{\text{j}})}$$is the Laplace transformation of the augmented actuator setpoint array. The element ${\text{u}}_{\text{j}}(\text{s})\in {ℂ}^{{\text{n}}_{\text{j}}}(\text{j}=1\dots {\text{N}}_{\text{u}})$ is the Laplace transformation of the j^th individual CD actuator setpoint profile, such asthe headbox slice, water spray, steam box, or induction heater. ${\text{N}}_{\text{u}}$is the total number of actuator beams available as MV’s, and ${\text{n}}_{\text{j}}$ is the number of individual zones of the jth actuator beam. In general a CD system has the same number of elements for all CD measurement profiles, but different numbers of actuator beam setpoints. D(s)$\in {ℂ}^{({\text{N}}_{\text{y}}\cdot \text{m})}$ is the Laplace transformation of the augmented process disturbance array. It represents process output disturbances.

${\text{G}}_{\text{ij}}(\text{s})\in {ℂ}^{\text{m}\times {\text{n}}_{\text{j}}}$ ($\text{i}=1\dots {\text{N}}_{\text{y}}\text{andj}=1\dots {\text{N}}_{\text{u}})$in (15) is the transfer matrix of the sub-system from the j^th actuator beam ${\text{u}}_{\text{j}}$ to the i^th CD quality measurement${\text{y}}_{\text{i}}$. The model of this sub-system can be represented by a spatial static matrix ${\text{P}}_{\text{ij}}\in {ℝ}^{\text{m}\times {\text{n}}_{\text{j}}}$ with a temporal dynamic transfer function${\text{h}}_{\text{ij}}(\text{s})$. In practice, ${\text{h}}_{\text{ij}}(\text{s})$is simplified as a first-order plus dead time system. Therefore, ${\text{G}}_{\text{ij}}(\text{s})$is given by

$${\text{G}}_{\text{ij}}(\text{s})={\text{P}}_{\text{ij}}{\text{h}}_{\text{ij}}(\text{s})={\text{P}}_{\text{ij}}\frac{1}{1+{\text{T}}_{\text{p}}\text{s}}{\text{e}}^{-{\text{T}}_{\text{d}}\text{s}}$$E21

where ${\text{T}}_{\text{p}}$ is the time constant and ${\text{T}}_{\text{d}}$ is the time delay. The static spatial matrix ${\text{P}}_{\text{ij}}$ is a matrix with ${\text{n}}_{\text{j}}$ columns, i.e., ${\text{P}}_{\text{ij}}=\left[\begin{array}{cccc}{\text{p}}^1& {\text{p}}^2& \cdots & {\text{p}}^{{\text{n}}_{\text{j}}}\end{array}\right]$and its kth column ${\text{p}}^{\text{k}}$ represents the spatial response of the k^th individual actuator zone of the j^th actuator beam. As proposed in (Gorinevsky & Gheorghe 2003), ${\text{p}}^{\text{k}}$can be formulated by,

$${\text{p}}^{\text{k}}=\frac{\text{g}}{2}\{{\text{e}}^{-\frac{\text{α}{((\text{x}-{\text{x}}_{\text{k}})-\text{βω})}^2}{{\text{ω}}^2}}\text{cos}(\frac{\text{π}}{\text{ω}}((\text{x}-{\text{x}}_{\text{k}})-\text{βω}))+{\text{ e}}^{-\frac{\text{α}{((\text{x}-{\text{x}}_{\text{k}})+\text{βω})}^2}{{\text{ω}}^2}}\text{cos}(\frac{\text{π}}{\text{ω}}((\text{x}-{\text{x}}_{\text{k}})+\text{βω}))\}$$E22

where x is the coordinate of CD measurements (CD bins), g is the process gain, $\text{ω}$is the response width, $\text{α}$is the attenuation and $\text{β}$ is divergence. ${\text{x}}_{\text{k}}$is the CD alignment that indicates the spatial relationship between the centre of the k^th individual CD actuator and the center of the corresponding measurement responses. A fuzzy function may be used to model the CD alignment. Refer to (Gorinevsky & Gheorghe 2003) for the technical details.

Figure17 illustrates the structure of the spatial response matrix${\text{P}}_{\text{ij}}$. The colour map on the left shows the band-diagonal property of${\text{P}}_{\text{ij}}$ ; and the plot in the right shows the spatial response of the individual spatial actuator${\text{p}}^{\text{k}}$. It can be seen that each individual actuator affects not only its own spatial zone area, but also adjacent zone areas.

4.1.2. Model identification

Model identification of the papermaking CD process is the procedure to determine the values of the parameters in (16, 17), i.e., the dynamic model parameters${\text{θ}}_{\text{T}}=\{{\text{T}}_{\text{p}},{\text{T}}_{\text{d}}\}$, the spatial model parameters${\text{θ}}_{\text{CD}}=\{\text{g},\text{ω},\text{α},\text{β}\}$, and the alignment${\text{x}}_{\text{k}}$.An iterative identification algorithm has been proposed in (Gorinevsky & Gheorghe 2003). As with MD model identification, this algorithm is an open-loop model identification approach. Identification experiment data are first collected by performing open-loop bump tests.

Figure 17.

The illustration to spatial response matrix${\text{P}}_{\text{ij}}$.

Figure18 illustrates the logic flow of this algorithm. This nontrivial system identification approach first estimates the overall dynamic response and spatial response, and subsequently identifies the dynamic model parameter ${\text{θ}}_{\text{T}}$ and the spatial model parameter${\text{θ}}_{\text{CD}}$. $\hat{\text{h}}$in Figure 18 is the estimated finite impulse response (FIR) of the dynamic model $\text{h}(\text{s})$ in (16). $\hat{\text{p}}$in Figure 18 is the estimated steady state measurement profile, i.e., overall spatial response. For easier notation, we omit the indexes i and j here. The key concept of the algorithm is to optimize the model parameters iteratively. Refer to (Gorinevsky & Gheorghe 2003) for technical details of this algorithm, and (Gorinevsky & Heaven, 2001) for the theoretical proof of the algorithm convergence.

Figure 18.

The schematic of the iterative CD system identification algorithm

The algorithm described above has been implemented in a software package, named IntelliMap^TM, which has been widely used in pulp and paper industries. The tool executes the open-loop bump testsautomatically and, at the end of the experiments, provides a continuous-time transfer matrix model (defined in (14)). For convenience, the MPC controller design discussed in the next section will use the state space model representation. Conversion of the continuous-time transfer matrix model into the discrete-time state space model is trivial (Chen 1999) and is omitted here.

4.2. CD-MPC design

In this section, a state space realization of (14) is used for the MPC controller development,

$$\begin{array}{cc}\text{X}(\text{k}+1)& =\text{AX}(\text{k})+\text{BΔU}(\text{k})\\ \text{Y}(\text{k})& =\text{CX}(\text{k})+\text{D}(\text{k})\end{array}$$E23

$\text{X}(\text{k})\in {ℝ}^{(2\cdot {\text{N}}_{\text{y}}\cdot \text{m})}$, $\text{Y}(\text{k})\in {ℝ}^{({\text{N}}_{\text{y}}\cdot \text{m})}$, $\text{ΔU}(\text{k})\in {ℝ}^{({\displaystyle \sum }_{\text{j}=1}^{{\text{N}}_{\text{u}}}{\text{n}}_{\text{j}})}$and$\text{D}(\text{k})\in {ℝ}^{({\text{N}}_{\text{y}}\cdot \text{m})}$ are the augmented state, output, actuator move, and output disturbance arrays of the papermaking CD process with multiple CD actuator beams and multiple quality measurement arrays. {A, B, C} are the model matrices with compatible dimensions. Assume (A, B) is controllable and (A, C) is observable. In this section, the objective function of CD-MPC is developed first. Then the CD actuator constraints are incorporated in the objective function. Finally a fast QP solver is presented for solving the large scale constrained CD-MPC optimization problem. How to tune a CD-MPC controller is also covered in this section

4.2.1. Objective function of CD-MPC

The first step of MPC development is performing the system output prediction over a certain length of prediction horizon. From the state space model defined in (18), we can predict the future states,

$$\begin{array}{cc}\mathcal{X}(\text{k})& ={\text{A}}_{}\text{X}(\text{k})+{\text{B}}_{}\text{Δ}\mathcal{U}(\text{k})\end{array}$$E24

where $\mathcal{X}(\text{k})\in {ℝ}^{(2\cdot {\text{N}}_{\text{y}}\cdot \text{m}\cdot {\text{H}}_{\text{p}})}$ is the state prediction, $\text{Δ}\mathcal{U}(\text{k})\in {ℝ}^{({\text{H}}_{\text{u}\cdot }{\displaystyle \sum }_{\text{j}=1}^{{\text{N}}_{\text{u}}}{\text{n}}_{\text{j}})}$is the augmented actuator moves. ${\text{A}}_{}$and ${\text{B}}_{}$ are the state and input prediction matrices with the compatible dimensions. ${\text{H}}_{\text{p}}{\text{andH}}_{\text{u}}$are the output and input prediction horizons, respectively.

The explicit expressions of the parameters in (19) are

$$\begin{array}{l}\mathcal{X}(\text{k})=\left[\begin{array}{c}\text{X}(\text{k}+1|\text{k})\\ \text{X}(\text{k}+2|\text{k})\\ ⋮\\ \text{X}(\text{k}+{\text{H}}_{\text{p}}|\text{k})\end{array}\right],{\text{A}}_{}=\left[\begin{array}{c}\text{A}\\ {\text{A}}^2\\ ⋮\\ {\text{A}}^{{\text{H}}_{\text{p}}}\end{array}\right],{\text{B}}_{}=\left[\begin{array}{ccc}\text{B}& \cdots & 0\\ \text{AB}& \cdots & 0\\ ⋮& \ddots & ⋮\\ {\text{A}}^{{\text{H}}_{\text{p}}-1}\text{B}& \cdots & {\text{A}}^{{\text{H}}_{\text{p}}-{\text{H}}_{\text{u}}}\text{B}\end{array}\right],\\ \text{and}\text{Δ}\mathcal{U}(\text{k})=\left[\begin{array}{c}\text{ΔU}(\text{k}|\text{k})\\ \text{ΔU}(\text{k}+1|\text{k})\\ ⋮\\ \text{ΔU}(\text{k}+{\text{H}}_{\text{u}}-1|\text{k})\end{array}\right]\end{array}$$E25

The initial state ${\hat{\text{X}}}^0(\text{k}|\text{k}-1)$ at instant k can be estimated from the previous state estimation $\hat{\text{X}}(\text{k}-1)$ and the previous actuator move$\text{ΔU}(\text{k}-1)$, i.e.,

$$\begin{array}{cc}{\hat{\text{X}}}^0(\text{k}|\text{k}-1)& =\text{A}\hat{\text{X}}(\text{k}-1)+\text{BΔU}(\text{k}-1)\end{array}$$E26

The measurement information at instant k can be used to improve the estimation,

$$\begin{array}{cc}\hat{\text{X}}(\text{k})& ={\hat{\text{X}}}^0(\text{k}|\text{k}-1)+\text{L}(\end{array}\text{Y}(\text{k})-\text{C}{\hat{\text{X}}}^0(\text{k}|\text{k}-1))$$E27

where L$\in {ℝ}^{(2\cdot {\text{N}}_{\text{y}}\cdot \text{m})\times ({\text{N}}_{\text{y}}\cdot \text{m})}$ is the state observer matrix.

Replace the state $\text{X}(\text{k})$ by its estimation$\hat{\text{X}}(\text{k})$, and perform the output prediction$\mathcal{Y}(\text{k})$,

$$\begin{array}{cc}\mathcal{Y}(\text{k})& ={\text{c}}_{}{\text{A}}_{}\hat{\text{X}}(\text{k})+{\text{c}}_{}{\text{B}}_{}\text{Δ}\mathcal{U}(\text{k})\end{array}$$E28

where ${\text{c}}_{}\in {ℝ}^{({\text{N}}_{\text{y}}\cdot \text{m}\cdot {\text{H}}_{\text{p}})\times (2\cdot {\text{N}}_{\text{y}}\cdot \text{m}\cdot {\text{H}}_{\text{p}})}$ is the output prediction matrix, given by

$$\begin{array}{c}{\text{c}}_{}=\text{diag}(\text{C},\cdots \text{C})=\left[\begin{array}{cccc}\text{C}& 0& \cdots & 0\\ 0& \text{C}& ⋮& 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & \text{C}\end{array}\right]\text{. Also},\mathcal{Y}(\text{k})=\left[\begin{array}{c}\text{Y}(\text{k}+1|\text{k})\\ \text{Y}(\text{k}+2|\text{k})\\ ⋮\\ \text{Y}(\text{k}+{\text{H}}_{\text{p}}|\text{k})\end{array}\right]\end{array}\text{.}$$E29

From the expression in (24), one can define the objective function of a CD-MPC problem,

$${\min }_{\text{Δ}\mathcal{U}(\text{k})}\left|\right|\mathcal{Y}(\text{k})-{\mathcal{Y}}_{\text{tgt}}|{|}_{{\mathcal{Q}}_1}^2+||\text{Δ}\mathcal{U}(\text{k})|{|}_{{\mathcal{Q}}_2}^2+\left|\right|\mathcal{U}(\text{k})-{\mathcal{U}}_{\text{tgt}}|{|}_{{\mathcal{Q}}_3}^2+||{ℱ}_{\text{b}}\mathcal{U}(\text{k})|{|}_{{\mathcal{Q}}_4}^2$$E30

${\mathcal{Y}}_{\text{tgt}}={[{\text{Y}}_{\text{tgt}}^{\text{T}},{\text{Y}}_{\text{tgt}}^{\text{T}},\cdots ,{\text{Y}}_{\text{tgt}}^{\text{T}}]}^{\text{T}}$defines the measurement targets over the prediction horizon${\text{H}}_{\text{p}}$. Similarly, ${\mathcal{U}}_{\text{tgt}}={[{\text{U}}_{\text{tgt}}^{\text{T}},{\text{U}}_{\text{tgt}}^{\text{T}},\cdots ,{\text{U}}_{\text{tgt}}^{\text{T}}]}^{\text{T}}$defines the input actuator setpoint targets over the control horizon${\text{H}}_{\text{u}}$. $({\mathcal{Q}}_1,{\mathcal{Q}}_2,{\mathcal{Q}}_3,{\mathcal{Q}}_4)$are the diagonal weighting matrices. ${\mathcal{Q}}_1$defines the relative importance of the individual quality measurements. ${\mathcal{Q}}_2$defines the relative aggressiveness of the individual CD actuators. ${\mathcal{Q}}_3$defines the relative deviation from the targets of the individual CD actuators. ${\mathcal{Q}}_4$defines the relative picketing penalty of the individual CD actuators. The matrix${ℱ}_{\text{b}}=\text{diag}({\text{F}}_{\text{b}},\cdots ,{\text{F}}_{\text{b}}$ ) is the augmented actuator bending matrix. The detailed definition of ${\text{F}}_{\text{b}}$ will be covered in Section 4.2.2. $\left|\right|\cdot |{|}_{{ℛ}_i}^2$is the square of weighted 2-norm, i.e., $\left|\right|\cdot |{|}_{{ℛ}_i}^2$= ${(\cdot )}^T{ℛ}_i(\cdot ).$ In general, $({\mathcal{Q}}_1,{\mathcal{Q}}_2,{\mathcal{Q}}_3,{\mathcal{Q}}_4)\text{are}$used as the tuning parameters for CD-MPC.$\mathcal{U}(\text{k})$is the future input prediction. It can be expressed by

$$\begin{array}{c}\begin{array}{c}\mathcal{U}(\text{k})=\left[\begin{array}{c}\text{U}(\text{k}|\text{k})\\ \text{U}(\text{k}+1|\text{k})\\ ⋮\\ \text{U}(\text{k}+{\text{H}}_{\text{u}-1})|\text{k})\end{array}\right]\\ \end{array}\begin{array}{c}=\\ \end{array}\begin{array}{c}\underset{}{\left[\begin{array}{c}\text{I}\\ \text{I}\\ ⋮\\ \text{I}\end{array}\right]}\\ {\mathcal{S}}_{\text{I}}\end{array}\begin{array}{c}\text{U}(\text{k}-1)\\ \end{array}\begin{array}{c}+\\ \end{array}\begin{array}{c}\underset{}{\left[\begin{array}{cccc}\text{I}& 0& \cdots & 0\\ \text{I}& \text{I}& ⋮& 0\\ ⋮& \cdots & \ddots & ⋮\\ \text{I}& \text{I}& \text{I}& \text{I}\end{array}\right]}\\ {\mathcal{S}}_{\text{p}}\end{array}\end{array}\begin{array}{c}\text{Δ}\mathcal{U}(\text{k})\\ \end{array}$$E31

where $\text{I}\in {ℝ}^{({\displaystyle \sum }_{\text{j}=1}^{{\text{N}}_{\text{u}}}{\text{n}}_{\text{j}})\times ({\displaystyle \sum }_{\text{j}=1}^{{\text{N}}_{\text{u}}}{\text{n}}_{\text{j}})}$ is the identity matrix. Inserting (26) into (25) and replacing $\mathcal{U}(\text{k})$by$\text{Δ}\mathcal{U}(\text{k})$, the QP problem can be recast into

$${\min }_{\text{Δ}\mathcal{U}(\text{k})}\frac{1}{2}\text{Δ}{\mathcal{U}}^{\text{T}}(\text{k})\text{ΦΔ}\mathcal{U}(\text{k})+{\text{φ}}^{\text{T}}\text{Δ}\mathcal{U}(\text{k})$$E32

where$\text{Φ}$ is the Hessian matrix and $\text{φ}$ is the gradient matrix. Both can be derived from the prediction matrices (${\text{A}}_{},{\text{B}}_{},{\text{C}}_{})$and weighting matrices (${\mathcal{Q}}_1,{\mathcal{Q}}_2,{\mathcal{Q}}_3,{\mathcal{Q}}_4)\text{.}$Refer to (Fan 2003) for the detailed expressions of $\text{Φ}$ and $\text{φ}\text{.}$

By solving the QP problem in (27), one can derive the predicted optimal array$\text{Δ}\mathcal{U}(\text{k})$. Only the first component of$\text{Δ}\mathcal{U}(\text{k})$, i.e., $\text{ΔU}(\text{k}),$is sent to the real process and the rest are rejected. By repetition of this procedure, the optimal MV moves at any instant are derived for unconstrained CD-MPC problems.

4.2.2.Constraints

In Section 4.2.1 the CD-MPC controller is formulated as an unconstrained QP problem. In practice the new actuator setpoints given by the CD-MPC controller in (27) should always respect the actuator’s physical limits. In other words, the hard constraints on $\text{Δ}\mathcal{U}(\text{k})$ should be added into the problem in (27).

The CD actuator constraints include:

• First and second order bend limits;

• Average actuator setpoint maintenance;

• Maximum actuator setpoints;

• Minimum actuator setpoints; and

• Maximum change of actuator setpoints between consecutive CD-MPC iterations.

Of these five types of actuator constraints, most of them are very common for the typical MPC controllers, except for the bend limits which are special for papermaking CD processes. The first and second bend limits define the allowable first and second order difference between the adjacent actuator setpoints of the actuator beam. It typically applies to slice lips and induction heaters to prevent the actuator beams from being overly bent or locally over-heated. The bending matrix of the j^th actuator beam, $F_{b,j}$($\text{j}=1,\cdots ,{\text{N}}_{\text{u}})$can be defined by

$$\begin{array}{c}-\\ \end{array}\begin{array}{c}\underset{}{\left[\begin{array}{c}{\text{δ}}_{1,\text{j}}\\ {\text{δ}}_{2,\text{j}}\\ ⋮\\ {\text{δ}}_{2,\text{j}}\\ {\text{δ}}_{1,\text{j}}\end{array}\right]}\\ {\text{γ}}_{\text{b},\text{j}}\end{array}\begin{array}{c}\le \\ \end{array}\begin{array}{cc}\begin{array}{c}\underset{}{\left[\begin{array}{cccccc}-1& 1& 0& \cdots 0& \cdots & 0\\ 0.5& -1& 0.5& 0& \cdots & 0\\ 0& 0.5& -1& 0.5& \cdots & 0\\ ⋮& \cdots & \cdots & \cdots & \ddots & ⋮\\ 0& 0& \cdots & 0.5& -1& 0\\ 0& 0& \cdots & 0& 1& -1\end{array}\right]}\\ {\text{F}}_{\text{b},\text{j}}\end{array}& \begin{array}{c}\underset{}{\left[\begin{array}{c}{\text{u}}_{1,\text{j}}\\ {\text{u}}_{2,\text{j}}\\ ⋮\\ {\text{u}}_{{\text{n}}_{\text{j}}-1,\text{j}}\\ {\text{u}}_{{\text{n}}_{\text{j}},\text{j}}\end{array}\right]}\\ {\text{u}}_{\text{j}}\end{array}\end{array}\begin{array}{c}\le \\ \end{array}\begin{array}{c}\underset{}{\left[\begin{array}{c}{\text{δ}}_{1,\text{j}}\\ {\text{δ}}_{2,\text{j}}\\ ⋮\\ {\text{δ}}_{2,\text{j}}\\ {\text{δ}}_{1,\text{j}}\end{array}\right]}\\ {\text{γ}}_{\text{b},\text{j}}\end{array}$$E33

where ${\text{δ}}_{1,\text{j}}$ and ${\text{δ}}_{2,\text{j}}$ are the first order and the second order bend limit of the j^th actuator beam${\text{u}}_{\text{j}}$. ${\text{γ}}_{\text{b}}$and$F_{b,j}$ define the bend limit vector and the bend limit matrix of the j^th actuator${\text{u}}_{\text{j}}$, respectively. The bend limit matrix $F_{b,j}$ is not only part of the constraints, but also the objective function in (27). In (27),${ℱ}_{\text{b}}=\text{diag}({\text{F}}_{\text{b}},\cdots ,{\text{F}}_{\text{b}})$ and${\text{F}}_{\text{b}}=\text{diag}({\text{F}}_{\text{b},1},\cdots ,{\text{F}}_{\text{b},{\text{N}}_{\text{u}}}$ ).

The individual bend limit constraint on the j^th actuator beam ${\text{u}}_{\text{j}}$ in (28) can be extended to the overall bend limit matrix ${\text{F}}_{\text{b}}$ for the augmented actuator setpoint array U, i.e.,

$$\left[\begin{array}{c}{\text{F}}_{\text{b}}\\ -{\text{F}}_{\text{b}}\end{array}\right]\text{U}\le \left[\begin{array}{c}{\text{γ}}_{\text{b}}\\ {\text{γ}}_{\text{b}}\end{array}\right]$$E34

where ${\text{γ}}_{\text{b}}$ is the overall bend limit vector, and${\text{γ}}_{\text{b}}={[{\text{γ}}_{\text{b},1}^{\text{T}},\cdots {\text{γ}}_{\text{b},{\text{N}}_{\text{u}}}^{\text{T}}]}^{\text{T}}$.

Similar to the bend limits, other types of actuator physical constraints can be formulated as the matrix inequalities,

$$\left[\begin{array}{c}{\text{F}}_{\max }\\ -{\text{F}}_{\min }\\ {\text{F}}_{\text{avg}}\\ -{\text{F}}_{\text{avg}}\\ {\text{F}}_{\Delta \text{U}}\\ -{\text{F}}_{\Delta \text{U}}\end{array}\right]\text{U}\le \left[\begin{array}{c}{\text{γ}}_{\max }\\ {\text{γ}}_{\min }\\ {\text{γ}}_{\text{avg}}\\ {\text{γ}}_{\text{avg}}\\ {\text{γ}}_{\Delta \text{U}}\\ {\text{γ}}_{\Delta \text{U}}\end{array}\right]$$E35

where the subscripts “max”, “min”, “avg”, and “$\Delta \text{U}"$stand for the maximum, minimum, average limit, and maximum setpoint changes between two consecutive CD-MPC iterations of the augmented actuator setpoint array, U. It is straightforward to derive the expressions of${\text{F}}_{\max }$, ${\text{F}}_{\min }$, ${\text{F}}_{\text{avg}}$,${\text{F}}_{\Delta U}$. Therefore the detailed discussion is omitted.

From (29) and (30), one can see that the constraints on the augmented actuator setpoint array $\text{U}$ can be represented by a linear matrix inequality, i.e.,

$$\text{FU}\le \text{γ}$$E36

where F and $\text{γ}$ are constant coefficients used to combine the inequalities in (29) and (30) together.

(26) is inserted into (31).The constraint in (31) is then added to the objective function in (27). Finally the CD-MPC controller is formulated as a constrained QP problem,

$$\begin{array}{cc}\text{subject to},& \begin{array}{c}{\min }_{\text{Δ}\mathcal{U}(\text{k})}\frac{1}{2}\text{Δ}{\mathcal{U}}^{\text{T}}(\text{k})\text{ΦΔ}\mathcal{U}(\text{k})+{\varphi }^{\text{T}}\text{Δ}\mathcal{U}(\text{k})\\ \\ \\ ℱ({\mathcal{S}}_{\text{I}}\text{U}(\text{k}-1)+{\mathcal{S}}_{\text{p}}\text{Δ}\mathcal{U}(\text{k}))\le \text{Γ}\end{array}\end{array}$$E37

where $ℱ=\text{diag}(\text{F},\text{F},\cdots ,\text{F})$ and$\text{Γ}=\text{diag}(\text{γ},\text{γ},\cdots ,\text{γ})$. By solving the QP problem in (32), the optimal actuator move at instant kcan be achieved.

4.2.3. CD-MPC tuning

Figure19 illustrates the implementation of the CD-MPC controller. First, the process model is identified offline from input/output process data. Then the CD-MPC tuning algorithm is executed to generate optimal tuning parameters. Subsequently these tuning parameters are deployed to the CD-MPC controller. The controller generates the optimal actuator setpoints continuously based on the feedback measurements.

Figure 19.

The implementation of the CD-MPC controller

The objective for CD-MPC tuning algorithm in Figure 19 is to determine the values of${\mathcal{Q}}_1,{\mathcal{Q}}_2,{\mathcal{Q}}_3,\text{and}{\mathcal{Q}}_4$ in (25). It has been proven that${\mathcal{Q}}_1$ defines the relative importance of quality measurements,${\mathcal{Q}}_2$ defines the dynamic characteristics of the closed-loop CD-MPC system, and${\mathcal{Q}}_3$ and${\mathcal{Q}}_4$ define the spatial frequency characteristics of the closed-loop CD-MPC system.${\mathcal{Q}}_3$ is for the high spatial frequencybehaviours and ${\mathcal{Q}}_4$ is for the low spatial frequencies (Fan 2004).

Strictly speaking, the CD-MPC tuning problem requires analyzing the robust stability of a closed-loop control system with nonlinear optimization. An analytic solution to the QP problem in (32) is the prerequisite for the CD-MPC tuning algorithm. However, in practice it is very challenging; almost impossible to derive the explicit solution to (32) due to the large size of CD-MPC problems. A novel two-dimensional loop shaping approach is proposed in (Fan 2004) to overcome limitations for large scaled MPC systems. The algorithm consists of four steps:

1. Ignore the inequality constraint in (32) such that the closed-loop system given by (27) is linear.

2. Compute the closed-loop transfer function of the unconstrained CD-MPC system given by (27).

3. By performing two-dimensional loop shaping, optimize the weighting matricesto get the best trade off between the performance and robustness of the unconstrained CD-MPC system.

4. Finally, re-introduce the constraint in (32) for implementation.

Figure 20 shows the closed-loop diagram of the unconstrained CD-MPC system with unstructured model uncertainties. The derivation of the pre-filtering matrix ${\text{K}}_{\text{r}}$ and feedback controller K is standard and can be found in (Fan 2003).

Figure 20.

Closed-loop diagram of unconstrained CD-MPC system with unstructured model uncertainties

From the small gain theory (Khalil 2001), the linear closed-loop system in Figure 20 is robustly stable if the closed-loop in (32) is nominally stable and,

$$\left|\right|{\text{G}}_{\text{ud}}(\text{z}){△}_{\text{G}}(\text{z})|{|}_{\infty }<1\Leftarrow \stackrel{}{\text{σ}}({\text{G}}_{\text{ud}}({\text{e}}^{\text{jω}}))<\frac{1}{\stackrel{}{\text{σ}}({△}_{\text{G}}({\text{e}}^{\text{jω}}))},\forall \text{ω}$$E38

Here ${\text{G}}_{\text{ud}}(\text{z})$ is the control sensitivity function which defines the linear transfer function from the output disturbance D(k) to the actuator setpoint U(k),

$${\text{G}}_{\text{ud}}(\text{z})=\text{K}(\text{z}){[\text{I}-\text{G}(\text{z})\text{K}(\text{z})]}^{-1}$$E39

The sensitivity function of the system in Figure 20 defines the linear transfer function from the output disturbance D(k) to the output Y(k),

${\text{G}}_{\text{yd}}(\text{z})={[\text{I}-\text{G}(\text{z})\text{K}(\text{z})]}^{-1}$.

By properly choosing the weighting matrices ${\mathcal{Q}}_1$ to${\mathcal{Q}}_4$, both the control sensitivity function ${\text{G}}_{\text{ud}}(\text{z})$ and the sensitivity function ${\text{G}}_{\text{yd}}(\text{z})$ can be guaranteed stable, and also the small gain condition in (33) can be satisfied. The two-dimensional loop shaping approach uses ${\text{G}}_{\text{ud}}(\text{z})$ and ${\text{G}}_{\text{yd}}(\text{z})$ to analyze the behaviour of the closed-loop system in Figure 20.

It has been shown that both ${\text{G}}_{\text{ud}}(\text{z})$ and ${\text{G}}_{\text{yd}}(\text{z})$ can be approximated as rectangular circulant matrices. One important property of circulant matrixes is that the circulant matrix can be block-diagonalized by left- and right-multiplying Fourier matrices. Fourier matrices multiplication is equivalent to performing the standard discrete Fourier transformation. Therefore, the two-dimensional frequency representation of ${\text{G}}_{\text{ud}}(\text{z})$ and ${\text{G}}_{\text{yd}}(\text{z})$ can be obtained by,

$${\hat{\text{g}}}_{\text{ud}}(\text{ν},{\text{e}}^{\text{jω}})={\text{F}}_{\text{n}}{\text{G}}_{\text{ud}}({\text{e}}^{\text{jω}}){\text{F}}_{\text{m}}^{\text{H}},\text{and}{\hat{\text{g}}}_{\text{yd}}(\text{ν},{\text{e}}^{\text{jω}})={\text{F}}_{\text{m}}{\text{G}}_{\text{yd}}({\text{e}}^{\text{jω}}){\text{F}}_{\text{m}}^{\text{H}}$$E40

where $\text{ν}$ represents the spatial frequency. ${\text{F}}_{\text{m}}$and ${\text{F}}_{\text{n}}$ are m-points and n-points Fourier matrices, respectively. The detailed definitions of Fourier matrices can be found in (Fan 2004). The two-dimensional frequent representation ${\hat{\text{g}}}_{\text{ud}}(\text{ν},{\text{e}}^{\text{jω}})$ and ${\hat{\text{g}}}_{\text{yd}}(\text{ν},{\text{e}}^{\text{jω}})$ are block diagonal matrices. The singular values of ${\hat{\text{g}}}_{\text{ud}}(\text{ν},{\text{e}}^{\text{jω}})$ and ${\hat{\text{g}}}_{\text{yd}}(\text{ν},{\text{e}}^{\text{jω}})$ are directly linked to the spatial frequencies.

Instead of tune ${\hat{\text{g}}}_{\text{ud}}(\text{ν},{\text{e}}^{\text{jω}})$ and ${\hat{\text{g}}}_{\text{yd}}(\text{ν},{\text{e}}^{\text{jω}})$ in full $\text{ν and ω}$ frequency ranges, two dimensional loop shaping approach decouples the spatial tuning and dynamic tuning by firstly tuning the controller at zero spatial frequency, i.e., setting$\text{ν}=0$, and then tuning the controller at zero dynamic frequency, i.e., setting$\text{ω}=0$. The theoretical proof of this strategy can be founded in (Fan 2004).

From spatial tuning, the value of the weighting matrices${\mathcal{Q}}_3\text{and}{\mathcal{Q}}_4$ can be determined, and from the dynamic tuning, the value of ${\mathcal{Q}}_2$ are determined.${\mathcal{Q}}_1$, as mentioned above, defines the relative importance of quality measurements and its value is defined by a CD-MPC user.

In practice, the process gain matrix ${\text{P}}_{\text{ij}}$ in (16) is ill-conditioned. Similar to MD-MPC tuning, the scaling matrices have to be applied before tuning the controller. A scaling approach discussed in (Lu 1996) is used by CD-MPC to reduce the condition number of the gain matrices.

4.2.4.Fast QP solver

The technical challenge of the CD-MPC optimization is how to solve the problem in (32) efficiently and accurately. The typical scanning rate of the paper machine is 10 - 30seconds. Also considering the time cost of software implementation and data acquisition, the computation time of the problem in (32) is typically limited to 5 to 10 seconds.

Different optimization techniques have been developed to solve QP problems efficiently, such as the active set method, interior point method, QR factorization, etc. This section presents a fast QP solver, called QPSchur, which is specifically designed to solve a large scaled CD-MPC problem. QPSchur is a dual space algorithm, where an unconstrained optimal solution is found first and violated constraints are added until the solution is feasible (Bartlett 2002).

Let’s consider the Lagrangian of the constrained QP in (32)

$$\text{Λ}(\text{Δ}\mathcal{U},\text{λ})=\frac{1}{2}\text{Δ}{\mathcal{U}}^{\text{T}}(\text{k})\text{ΦΔ}\mathcal{U}(\text{k})+{\varphi }^{\text{T}}\text{Δ}\mathcal{U}(\text{k})+{\text{λ}}^{\text{T}}({\text{ξ}}^{\text{T}}\text{Δ}\mathcal{U}(\text{k})-\text{ψ})$$E41

where $\text{ξ}$ = ${\mathcal{S}}_{\text{p}}^{\text{T}}{ℱ}^{\text{T}}$ and $\text{ψ}=\text{Γ}-ℱ{\mathcal{S}}_{\text{I}}\text{U}(\text{k}-1)\text{.}$ In (36), $\text{Δ}\mathcal{U}(\text{k})$is called the primary variable and $\text{λ}\le 0$ is called as the dual variable (also known as the Lagrangian variable).

At the starting point, QPSchur ignores all the constraints in (32) and solves unconstrained QP problem. This is equivalent to set the dual variable$\text{λ}=0$. By this means, the initial optimal solution $\text{Δ}{\mathcal{U}}^{\text{*}}(\text{k})$ is determined,

$$\text{Δ}{\mathcal{U}}^{\text{*}}(\text{k})=-{\text{Φ}}^{-1}\varphi$$E42

If $\text{Δ}{\mathcal{U}}^{\text{*}}(\text{k})$ satisfies all the inequality constraints, i.e., ${\text{ξ}}^{\text{T}}\text{Δ}{\mathcal{U}}^{\text{*}}(\text{k})\le \text{ψ}$, then $\text{Δ}{\mathcal{U}}^{\text{*}}(\text{k})$ is the optimal solution, i.e.$\text{Δ}{\mathcal{U}}^{\text{o}}(\text{k})=\text{Δ}{\mathcal{U}}^{\text{*}}(\text{k})$. The first elements of ${\text{ΔU}}^{\text{o}}(\text{k})$ are sent to the real process, and the CD-MPC optimization stops the search iteration. If $\text{Δ}{\mathcal{U}}^{\text{*}}(\text{k})$ violates one or more of the inequality constraints in (32), all the violation inequalities are noted, such that

$${\text{ξ}}_{\text{sub}}^{\text{T}}\text{Δ}{\mathcal{U}}^{\text{*}}(\text{k})\ge {\text{ψ}}_{\text{sub}}$$E43

where$({\text{ξ}}_{\text{sub}}^{\text{T}}$, ${\text{ψ}}_{\text{sub}})$is the violating subset of the inequality constraints in (32), and called the active set matrix and the active set vector, respectively. The Lagrangian in (36) is redefined by using$({\text{ξ}}_{\text{sub}}^{\text{T}}$,${\text{ψ}}_{\text{sub}})$. The Karush-Kuhn-Tucker (KKT) condition of the updated Lagrangian is,

$$\left[\begin{array}{cc}\text{Φ}& {\text{ξ}}_{\text{sub}}^{}\\ {\text{ξ}}_{\text{sub}}^{\text{T}}& \text{Ξ}\end{array}\right]\left[\begin{array}{c}\text{Δ}\mathcal{U}(\text{k})\\ {\text{λ}}_{\text{sub}}\end{array}\right]=\left[\begin{array}{c}-\varphi \\ {\text{ψ}}_{\text{sub}}\end{array}\right]$$E44

Here $\text{Ξ}=0$ for the first searching iteration. Since $\text{Φ}$ is non-singular (refer to Fan 2003), the problem in (39) can be solved by using Gaussian elimination. The Schur complement of the block $\text{Ξ}$ is given by

$$\mathbb{S}=\text{Ξ}-{\text{ξ}}_{\text{sub}}^{\text{T}}{\text{Φ}}^{-1}{\text{ξ}}_{\text{sub}}$$E45

The Schur complement theorem guarantees that $\mathbb{S}$ is non-singular if the Hessian matrix $\text{Φ}$ is non-singular. From$\mathbb{S}$, (39) can be solved by

$$\begin{array}{rr}\hfill {\text{λ}}_{\text{sub}}& \hfill ={\mathbb{S}}^{-1}({\text{ψ}}_{\text{sub}}+{\text{ξ}}_{\text{sub}}^{\text{T}}{\text{Φ}}^{-1}\varphi )\\ \hfill & \hfill ={\mathbb{S}}^{-1}({\text{ψ}}_{\text{sub}}+{\text{ξ}}_{\text{sub}}^{\text{T}}\text{Δ}{\mathcal{U}}^{\text{*}}(\text{k}))\\ \hfill & \hfill \\ \hfill \text{Δ}\mathcal{U}(\text{k})& \hfill ={\text{Φ}}^{-1}(-\varphi -{\text{ξ}}_{\text{sub}}{\text{λ}}_{\text{sub}})\end{array}$$E46

The inequality constraints in (32)are re-evaluated, and the new active constraints (violated constraints) and the positive dual variables inequalities are added into the subset pair$({\text{ξ}}_{\text{sub}}^{\text{T}}$,${\text{ψ}}_{\text{sub}})$. The KKT condition of (39) is updated to derive

$$\left[\begin{array}{cc}\text{Φ}& {\widehat{\text{ξ}}}_{\text{sub}}^{}\\ {\widehat{\text{ξ}}}_{\text{sub}}^{\text{T}}& \widehat{\text{Ξ}}\end{array}\right]\left[\begin{array}{c}\text{Δ}\mathcal{U}(\text{k})\\ {\widehat{\text{λ}}}_{\text{sub}}\end{array}\right]=\left[\begin{array}{c}-\varphi \\ {\widehat{\text{ψ}}}_{\text{sub}}\end{array}\right]$$E47

where

$${\widehat{\text{ξ}}}_{\text{sub}}^{}=[{\text{ξ}}_{\text{sub}},{\text{ξ}}_{\text{new}}],\widehat{\text{Ξ}}=\left[\begin{array}{cc}\text{Ξ}& \text{ρ}\\ {\text{ρ}}^{\text{T}}& \text{χ}\end{array}\right],{\widehat{\text{λ}}}_{\text{sub}}=\left[\begin{array}{c}{\text{λ}}_{\text{sub}}\\ {\text{λ}}_{\text{new}}\end{array}\right],\text{and}{\widehat{\text{ψ}}}_{\text{sub}}=\left[\begin{array}{c}{\text{ψ}}_{\text{sub}}\\ {\text{ψ}}_{\text{new}}\end{array}\right]$$E48

In the same fashion, the Schur complement of the block $\widehat{\text{Ξ}}$ can be represented by,

$$\begin{array}{cc}\hat{\mathbb{S}}& =\widehat{\text{Ξ}}-{\widehat{\text{ξ}}}_{\text{sub}}^{\text{T}}{\text{Φ}}^{-1}{\widehat{\text{ξ}}}_{\text{sub}}\\ & =\left[\begin{array}{cc}\text{Ξ}& \text{ρ}\\ {\text{ρ}}^{\text{T}}& \text{χ}\end{array}\right]-\left[\begin{array}{c}{\text{ξ}}_{\text{sub}}^{\text{T}}\\ {\text{ξ}}_{\text{new}}^{\text{T}}\end{array}\right]{\text{Φ}}^{-1}[{\text{ξ}}_{\text{sub}},{\text{ξ}}_{\text{new}}]\\ & =\left[\begin{array}{cc}\mathbb{S}& \text{ρ}-{\text{ξ}}_{\text{sub}}^{\text{T}}{\text{Φ}}^{-1}{\text{ξ}}_{\text{new}}\\ {\text{ρ}}^{\text{T}}-{\text{ξ}}_{\text{new}}{\text{Φ}}^{-1}{\text{ξ}}_{\text{sub}}^{\text{T}}& \text{χ}-{\text{ξ}}_{\text{new}}^{\text{T}}{\text{Φ}}^{-1}{\text{ξ}}_{\text{new}}\end{array}\right]\end{array}$$E49

From (44), the new Schur complement $\hat{\mathbb{S}}$ can be easily derived from$\mathbb{S}$. The Schur complement update requires only multiplication with ${\text{Φ}}^{-1}$ that is calculated in the initial search step and stored for reuse. This feature makes the SchurQP much faster than a standard QP solver. Removing the non-active constraints (zero dual variables) of each search step is achieved easily: the columns of the Schur complement $\hat{\mathbb{S}}$ corresponding to the non-active constraints is removed before pursuing the next search iteration.

At the current search iteration, if all the inequality constraints in (32) and the sign of dual variables are satisfied, the solution to (42) will be the final optimal solution of the CD-MPC controller, i.e.,

$$\begin{array}{rr}\hfill {\widehat{\text{λ}}}_{\text{sub}}& \hfill ={\hat{\mathbb{S}}}^{-1}({\widehat{\text{ψ}}}_{\text{sub}}+{\widehat{\text{ξ}}}_{\text{sub}}^{\text{T}}\text{Δ}{\mathcal{U}}^{\text{*}}(\text{k}))\\ \hfill \text{Δ}{\mathcal{U}}^{\text{o}}(\text{k})& \hfill ={\text{Φ}}^{-1}(-\varphi -{\widehat{\text{ξ}}}_{\text{sub}}{\widehat{\text{λ}}}_{\text{sub}})\end{array}$$E50

${\text{ΔU}}^{\text{o}}(\text{k})$ (the first component of the optimal solution$\text{Δ}{\mathcal{U}}^{\text{o}}(\text{k})$ ) is sent to the real process, and a new constrained QP problem is formed at the end of the next scan.

4.3. Mill implementation results

CD-MPC has been implemented in Honeywell’s quality control system (QCS) and widely deployed on different types of paper mills including fine paper, newsprint, liner board, and tissue, etc. In this chapter, a CD-MPC application for a fine paper machine will be used as an example to demonstrate the effectiveness of the CD-MPC controller.

4.3.1. Paper machine configuration

The paper machine discussed here is a fine paper machine, equipped with three CD actuator beams and two measurement scanner frames. The CD actuators include headbox slice lip (63 zones), infrared dryer (40 zones), and induction heater (79 zones). The two scanner frames hold the paper quality gauges for dry weight, moisture, and caliper. Each measurement profile includes 250 measurement points with the measurement interval equal to 25.4 mm (CD bin width). The production range of this machine is from 26 gsm (gram per square meter) to 85 gsm. The machine speed varies from 2650 feet per minute (13.5 meter/second) to 3100 feet per minute (15.7 meter/second). The scanning rates of the two scanners are 32 and 34 seconds, respectively. In order to capture the nonlinearity of the process, three model groups are setup to represent the products of light weight paper, medium weight paper, and heavy weight paper, respectively. All three CD actuator beams and three quality measurement profiles are included into the CD-MPC controller. In this section, the medium weight scenario is used to illustrate the control performance of the CD-MPC controller.

4.3.2. Multiple actuator beams and multiple quality measurements model

Figure21 shows the two-dimensional process models from the slice lip actuators (Autoslice) to the measurements of dry weight, moisture and caliper profiles. The system identification algorithm discussed in Section 4.1.2 is used to derive these models. The plots on the left are the spatial responses, and the plots on the right are the dynamic responses. The purple profiles are the average of the real process data, and the white profiles are the estimated profiles based on identified process model. It can be seen from comparison to the model for Autoslice to caliper that the models for Autoslice to dry weight and to moisture have high model fit. In general, the bump test with a larger bump magnitude and longer bump duration will lead to a more accurate process model (better model fit). However, the open-loop bump tests degrade the quality of the finished product and excessive bump tests are always prevented. The criterion of the CD model identification is to provide a process model accurate enough for a CD-MPC controller.

From the model identification results in Figure21, we can see the strong input-output coupling properties of papermaking CD processes. The response width from slice lip to dry weight equals to 226.8mm. This is equivalent to 2.3 times the zone width of the slice lip CD actuator. Therefore, each individual zone of the slice lip affects not only its own spatial zone but also adjacent zones. As we discussed above, a CD-MPC process has two-fold process couplings: one is the coupling between different actuator beams; and the other is the coupling between the different zones of the same actuator beams. Considering these strong coupling characteristics, MPC strategy is a good candidate for CD control design.

Figure 21.

Themultiple CD actuator beams and quality measurement model display

4.3.3. Control performance of the CD-MPC controller

Table 3 summarizes the performance comparison between the CD-MPC controller and the traditional single-input-single-output (SISO) CD controller (a Dahlin controller). Although traditional CD control is still quite common in paper mills, CD-MPC is becoming more and more popular. The significant performance improvement can be observed after switching CD control into multivariable CD-MPC.

Paper Properties	Traditional CD; Control 2σ	Multivariable CD-MPC Control 2σ	Improvement (%)
Dry Weight (gsm)	0.40	0.24	40%
Moisture (%)	0.31	0.19	39%
Caliper (mil)	0.032	0.025	22%

Table 3.

Traditional CD versus CD-MPC

Figures 22–24 provide a visual performance comparison for the different quality measurements in both spatial domain and spatial frequency domain. It can be seen that the peak-to-peak values (the proxy of$2{\text{σ}}_{\text{CD}}$ indexes) are smaller when using the CD-MPC controller. Also the controllable disturbances (the disturbances with the spatial frequency less than Xc) are effectively rejected by the CD-MPC controller. Here $\text{X}3\text{db}$ represents the spatial frequency where the spatial process power drops to 50% of the maximum spatial power over the full spatial frequency band, $\text{Xc}$represents the frequency where the spatial power drops to 4% of the maximum power, and $1/2\text{Xa}$ represents the Nyquist frequency.

Figure 22.

Performance comparison of dry weight profiles

Figure 23.

Performance comparison of moisture profiles

Figure 24.

Performance comparison of caliper profiles

5. Conclusion and perspective

We have seen that MPC has a number of applications in paper machine control. MPC performs basic MD control, and allows for enhanced MD-MPC control that incorporates economic optimization, and orchestrates transitions between paper grades. MPC can also be used for CD controls, using a carefully chosen solution technique to handle the large scale nature of the problem within the required time scale.

While MD-MPC provides robust and responsive control, and also easily scales to demanding paper machine applications with larger numbers of CV’s and MV’s. The MD-MPC formulation may also be augmented with an economic objective function so that paper machine operational efficiency can be optimized (maximum production, minimum energy costs, maximizing filler to fibre ratio etc.) while all quality variables continue to be regulated.

In the future as new online sensors, such as the extensional stiffness sensor, gain acceptance additional quality variables can be adding to MD-MPC. In the case of extensional stiffness, this online strength measurement could allow economic optimization to minimize fibre use while maintaining paper strength.

The papermaking CD process is a large scaled two-dimensional system. It shows strong input-output coupling properties. MPC is a standard technique in controlling multivariable systems, and has become a standard advanced control strategy in papermaking systems. However, there are several barriers for the acceptance of CD-MPC by mill personnel: one is the novel multivariable control concept and the other is the non-trivial tuning technique. Commercial offline tools, such as IntelliMap, facilitate the acceptance of CD-MPC by providing automatic model identification and easy-to-use offline CD-MPC tuning. Such packages enable the CD-MPC users to review the predicted CD steady states before they update their CD control to CD-MPC (Fan et al. 2005). CD-MPC has been successfully deployed in over 70 paper mills and applied to practically all types of existing CD processes from fine paper, to board, to newsprints, to tissues, etc. Without doubt, CD-MPC will have a significant impact in papermaking CD control applications over the next decade.

CD-MPC offers the significant capability to include multiple CD actuator arrays and multiple CD measurement arrays into one single CD controller. The next generation CD-MPC applications are most likely to include non-standard CD measurement, such as fibre orientation, gloss, web formation, and web porosity into the existing CD-MPC framework. A successful CD-MPC application for fibre orientation control has been reported in (Chu et al. 2010a). However there still exist technical challenges of controlling non-standard paper properties by using CD-MPC; for example, the derivation of accurate parametric models and the effectiveness of CD-MPC tuners for non-standard CD measurements.

In the current CD-MPC framework, system identification and controller design are clearly separated. The efforts towards integrating system identification and controller design may bring significant benefits to CD control. Online CD model identification has drawn extensive attention in both academia and industries. A closed-loop CD alignment identification algorithm is presented in (Chu et al. 2010b). Closed loop identification of the entire CD model remains an open problem.