Method of Designing Model Predictive Control for Cross Directional Flat Sheet Manufacturing Processes to Guarantee Temporal Robust Stability and Performance

Abstract

Automated parameter tuning techniques for cross-directional model predictive control for paper-making under user-specified parametric uncertainties are developed. The CD-MPC design explicitly accounts for parametric model uncertainty while finding a value for the CD profile trajectory generation that minimizes the CD settling time of the measurement 2 sigma spread and does not exceed an overshoot limit for actuator 2 sigma spread. The inventive technique includes: (i) providing uncertainty specifications for the temporal parameters of the process model, (ii) specifying 2 sigma overshoot limits for the CD actuator profile, (iii) using robust stability theory to find a minimum bound for the profile trajectory tuning parameter, (iv) using a frequency domain technique to reduce the search range for the profile trajectory tuning parameter, and (v) performing an intelligent search for the tuning parameter that minimizes measurement 2 sigma settling time without exceeding the actuator 2 sigma overshoot limit.

Claims

1. A system which forms a material in a spatially-distributed multivariable-array cross-directional process wherein the system comprises: at least one set of actuator arrays each distributed adjacent to the material in the cross direction (CD), wherein each set of actuator arrays is controllable to vary the properties of the material; means for measuring and acquiring properties data about the properties of the material, and a multivariable model predictive controller (MPC) for providing CD control of the CD process, wherein the MPC, in response to signals that are indicative of the properties of the material, provides signals to the at least one set of actuator arrays to control properties of the material, wherein the MPC includes means for temporally tuning the MPC which comprises: (a) means for inputting nominal temporal model parameters and corresponding parametric uncertainty specifications; (b) means for specifying a worst-case 2 sigma (2) overshoot limit for a CD actuator profile; (c) means for developing a robust temporal stability condition based on the parametric uncertainty specifications; (d) means for calculating a lower bound for a profile trajectory tuning parameter that assures robust stability; (e) means for reducing a possible range of the tuning parameter by using a frequency domain technique; and (f) means for performing an intelligent search for the tuning parameter that minimizes worst-case measurement 2 settling time without exceeding the actuator worst-case 2 overshoot limit.

2. The system of claim 1 wherein the nominal temporal model parameters comprise parameters of a temporal mathematical model of the process.

3. The system of claim 2 wherein the parametric uncertainty comprises trust ranges around the nominal temporal model parameters which characterize possible mismatch between the temporal mathematical model and real process.

4. The system of claim 1 wherein the worst-case 2 overshoot limit comprises the largest overshoot of all possible 2 spreads given parametric uncertainty specifications.

5. The system of claim 1 wherein the profile trajectory tuning parameter comprises a parameter in a temporal filter that affects the aggressiveness of a reference trajectory.

6. The system of claim 1 wherein worst-case measurement 2 settling time comprises the longest settling time of all possible 2 spreads given parametric uncertainty specifications.

7. The system of claim 1 wherein the robust temporal stability condition is characterized by a requirement that on the temporal frequency response of the closed-loop system which guarantees robust stability under the parametric uncertainty.

8. The system of claim 1 wherein the lower bound of the profile trajectory tuning parameter comprises the smallest value of the tuning parameter that guarantees the robust stability.

9. The system of claim 1 wherein the possible range of the tuning parameter comprises a range of tuning parameters within which an optimal tuning parameter exists.

10. The system of claim 1 wherein the MPC includes means for temporally tuning the MPC with respect to one of the actuator arrays.

11. In a process control system having a multivariable model predictive controller (MPC) for providing control to a spatially-distributed multiple-array, sheetmaking cross-directional (CD) process having at least one manipulated actuator array and at least one controlled measurement array, a method of providing control of the multiple-array process that comprises the steps of: (a) tuning the MPC by the steps of: (i) inputting nominal temporal model parameters and corresponding parametric uncertainty specifications; (ii) specifying a worst-case 2 sigma (2) overshoot limit for a CD actuator profile; (iii) developing a robust temporal stability condition based on the parametric uncertainty specifications; (iv) calculating a lower bound for a profile trajectory tuning parameter that assures robust stability; (v) reducing a possible range of the tuning parameter by using a frequency domain technique; and (vi) performing an intelligent search for the tuning parameter that minimizes worst-case measurement 2 settling time without exceeding the actuator worst-case 2 overshoot limit; (b) inputting the tuning parameters into the MPC; and (c) controlling the multiple-array CD process with the MPC.

12. The system of claim 11 wherein the nominal temporal model parameters comprise parameters of a temporal mathematical model of the process.

13. The system of claim 12 wherein the parametric uncertainty comprises trust ranges around the nominal temporal model parameters which characterize possible mismatch between the temporal mathematical model and real process.

14. The system of claim 11 wherein the worst-case 2 overshoot limit comprises the largest overshoot of all possible 2 spreads given parametric uncertainty specifications.

15. The system of claim 11 wherein the profile trajectory tuning parameter comprises a parameter in a temporal filter that affects the aggressiveness of a reference trajectory.

16. The system of claim 11 wherein worst-case measurement 2 settling time comprises the longest settling time of all possible 2 spreads given parametric uncertainty specifications.

17. The system of claim 11 wherein the robust temporal stability condition is characterized by a requirement that on the temporal frequency response of the closed-loop system which guarantees robust stability under the parametric uncertainty.

18. The system of claim 11 wherein the lower bound of the profile trajectory tuning parameter comprises the smallest value of the tuning parameter that guarantees the robust stability.

19. The system of claim 11 wherein the possible range of the tuning parameter comprises a range of tuning parameters within which an optimal tuning parameter exists.

20. A non-transitory computer readable medium embodying a computer program for automatically tuning a model predictive controller (MPC) employed to control a cross-directional process having a manipulated actuator array comprising a plurality of actuators and at least one controlled measurement array wherein the program comprises readable program code for: (a) inputting nominal temporal model parameters and corresponding parametric uncertainty specifications; (b) specifying a worst-case 2 sigma (2) overshoot limit for a CD actuator profile; (c) developing a robust temporal stability condition based on the parametric uncertainty specifications; (d) calculating a lower bound for a profile trajectory tuning parameter that assures robust stability; (e) reducing a possible range of the tuning parameter by using a frequency domain technique; and (f) performing an intelligent search for the tuning parameter that minimizes worst-case measurement 2 settling time without exceeding the actuator worst-case 2 overshoot limit;

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0033] FIGS. 1, 2, and 3 are schematic illustrations of a papermaking system;

[0034] FIG. 4 is a block diagram of a closed-loop CD-MPC system;

[0035] FIG. 5 is a block diagram of a closed-loop CD-MPC system;

[0036] FIGS. 6(a) and 6(b) illustrate robust performance visualization for the 2 spreads for the output and input, respectively;

[0037] FIG. 7 illustrates the minimum limit on the maximum peak on the sensitivity function;

[0038] FIG. 8 is a flow chart of the automated temporal tuning;

[0039] FIGS. 9(a) and 9(b) illustrate visualization technique for both worst- and best-case 2 spreads;

[0040] FIGS. 10(a) and 10(b) are the temporal tuning result for the proposed approach; and

[0041] FIGS. 11(a) to 11(b) are performance illustrations for the proposed robust temporal tuning method.

DETAILED DESCRIPTION OF PREFERRED EMBODIMENTS

[0042] As shown in FIG. 1, a system for producing continuous sheet material includes various processing stages such as headbox 10, steambox 12, a calendaring stack 14 and reel 16. The array of actuators 18 in headbox 10 controls the discharge of wet stock (or feedstock) material through a plurality of slices onto supporting web or wire 30 which rotates between rollers 22 and 24. Similarly, actuators 20 on steambox 12 can control the amount of steam that is injected at points across the moving sheet. Sheet material exiting the wire 30 passes through a dryer 34 which includes actuators 36 that can vary the cross directional temperature of the dryer. A scanning sensor 38, which is supported on supporting frame 40, continuously traverses and measures properties of the finished sheet in the cross direction. The finished sheet product 42 is then collected on reel 16. As used herein, the wet end portion of the system includes the headbox, the web, and those sections just before the dryer, and the dry end comprises the sections that are downstream from the dryer. Typically, the two edges of the wire in the cross direction are designated front and back (alternatively, referred as the high and low) with the back side being adjacent to other machinery and less accessible than the front side.

[0043] The system further includes a profile analyzer 44 that is connected, for example, to scanning sensor 38 and actuators 18, 20, 32 and 36 on the headbox 10, steam box 12, vacuum boxes 28, and dryer 34, respectively. The profile analyzer is a computer which includes a control system that operates in response to the cross-directional measurements from scanner sensor 38. In operation, scanning sensor 38 provides the analyzer 44 with signals that are indicative of the magnitude of a measured sheet property, e.g., caliper, dry basis weight, gloss or moisture, at various cross-directional measurement points. The analyzer 44 also includes software for controlling the operation of various components of the sheetmaking system, including, for example, the above described actuators. To implement to the control system of the present invention, analyzer 44 can include memory 62 and processing devices 64 to execute software/firmware instructions for performing various operations related to MPC control of an industrial process. Interface 60 allows processing devices to receive data and provide signals to actuators or controllers.

[0044] FIG. 2 depicts a slice lip control system which is mounted on a headbox 10 for controlling the extent to which a flexible slice lip member 46 extends across the discharge gap 48 at the base of the headbox 10. The slice lip member 46 extends along the headbox 10 across the entire width of the web in the cross-direction. The actuator 18 controls of the slice lip member 46, but it should be understood that the individual actuators 18 are independently operable. The spacing between the individual actuators in the actuator array may or may not be uniform. Wetstock 50 is supported on wire 30 which rotates by the action of rollers 22 and 24.

[0045] As an example shown in FIG. 3, the amount of feedstock that is discharged through the gap between the slice lip member and the surface of the web 30 of any given actuator is adjustable by controlling the individual actuator 18. The feed flow rates through the gaps ultimately affect the properties of the finished sheet material, i.e., the paper 42. Specifically, as illustrated, a plurality of actuators 18 extend in the cross direction over web 30 that is moving in the machine direction indicated by arrow 6. Actuators 18 can be manipulated to control sheet parameters in the cross direction. A scanning device 38 is located downstream from the actuators and it measures one or more the properties of the sheet. In this example, several actuators 18 are displaced as indicated by arrows 4 and the resulting changes in sheet property is detected by scanner 38 as indicated by the scanner profile 54. By averaging many scans of the sheet, the peaks of profile 54 indicated by arrows 56 can be determined. This type of operation is typically used in traditional open and closed-loop bump tests. It should be noted that besides being positioned in the headbox, actuators can be placed at one or more strategic locations in the papermaking machine including, for example, in the steamboxes, dryers, and vacuum boxes. The actuators are preferably positioned along the CD at each location.

[0046] It is understood that the inventive technique is sufficiently flexible as to be applicable for online implementation with any large-scale industrial multiple actuator array and multiple product quality measurements cross-directional process that is controlled by a single-input-single-output (SISO) controller or by a multivariable model predictive controller (MPC) such as in papermaking. Suitable paper machine processes where paper is continuously manufactured from wet stock are further described, for instance, in U.S. Pat. No. 6,807,510 to Backstrom et al., and U.S. Pat. No. 8,224,476 to Chu et al., and U.S. 2015/0268645 to Shi et al., which are incorporated herein by reference. The present invention employs so-called bump tests by disturbing selected actuators and detecting their responses. In particular, operating parameter on the sheetmaking system, such as a papermaking machine, is altered and changes of certain dependent variables resulting therefrom are measured. Bump tests techniques are described in U.S. Pat. No. 6,086,237 to Gorinevsky et al. and U.S. Pat. No. 7,459,060 to Stewart, which are incorporated herein by reference. While the invention will be described with respect to a paper-making machine, it is understood that the invention is applicable to industrial plastic sheetmaking, rubber sheetmaking, and sheet metal operations.

I. CD-MPC Structure

[0047] As shown in FIG. 4, the closed-loop CD-MPC structure for a paper-making process consists of four parts: the real process G.sub.p(z) (74), the nominal model G(z) (76), the MPC (72) and the temporal filter F.sub.(z) (70), which are illustrated in detail in the following. Moreover, the associated signals are defined as follows: y.sub.tgt, y.sub.sp, u(z), d(z), y.sub.p(z), and y(z) are the output target, the reference trajectory, the actuator profile, the disturbance profile, the measurement profile, and the predicted output profile, respectively.

Nominal Model

[0048] The nominal model G(z) of a CD paper-making process is characterized by

[00001]$\begin{array}{cc}G\left(z\right)=G_0.Math.h\left(z\right),h\left(z\right)=\frac{\left(1-a\right).Math.z^{-t_d}}{1-{\mathrm{az}}^{-1}},& \left(1\right)\end{array}$

where G.sub.0 is a constant matrix that characterizes the spatial response/gain of the CD process; h(z) is the temporal transfer function of the process, in which a and t.sub.d are the time constant and time delay in the discretized version.

[0049] The spatial gain matrix G.sub.0 has the parameterized structure as shown below:

[00002]$\begin{array}{cc}{G}_0=\left[g_1,\mathrm{.Math.}.Math.,g_n\right]{}^{mn},.Math.\begin{array}{c}{g}_k=.Math.f\left(x,,,,,c_k\right)\\ =.Math.\frac{}{2}.Math.\left\{\begin{array}{c}{e}^{-\frac{{\left(x-c_k+.Math..Math.\right)}^2}{{}^2}}.Math.\cos \left[\frac{\left(x-c_k+.Math..Math.\right)}{}\right]+\\ e^{-\frac{{\left(x-c_k-.Math..Math.\right)}^2}{{}^2}}.Math.\cos \left[\frac{\left(x-c_k-.Math..Math.\right)}{}\right]\end{array}\right\},\end{array}& \left(2\right)\\ x=1,\mathrm{.Math.}.Math.,m,k=1,\mathrm{.Math.}.Math.,n,& \end{array}$

[0050] where , , , and are the process gain, attenuation, width, and divergence, respectively. They are utilized to characterize the spatial response of each specific actuator. For the k.sup.th actuator, c.sub.k is the alignment parameter that determines the center of the corresponding spatial response.

Model Uncertainty

[0051] Since model-plant mismatch is unavoidable in process operation and identification, model uncertainties are considered in this disclosure. Based on traditional definitions in robust control to represent model mismatch, it is assumed that the real process model belongs to a set of possible models, characterized by an unstructured or parametric perturbation on the nominal model in (1). As parametric uncertainty is easy to understand and specify by the end users in the pulp and paper industry, the real process model G.sub.p(z) is then described by

[00003]$\begin{array}{cc}{G}_p\left(z\right)=G_0^p.Math.h_p\left(z\right),.Math.G_0^p=\left[g_{1.Math.p},\mathrm{.Math.}.Math.,g_{\mathrm{np}}\right],.Math.g_{\mathrm{kp}}=f\left(x,{}_p,{}_p,{}_p,{}_p,c_k^p\right),.Math.{}_p=\left(1+r_{}\right).Math.,{}_p=\left(1+r_{}\right).Math.,.Math.{}_p=\left(1+r_{}\right).Math.,{}_p=\left(1+r_{}\right).Math.,.Math.c_k^p=c_k+\mathrm{.Math.},.Math.h_p\left(z\right)=\frac{\left(1-a_p\right).Math.z^{-t_{\mathrm{dp}}}}{1-a_p.Math.z^{-1}},.Math.a_p=\left(1+r_a\right).Math.a,t_{\mathrm{dp}}=\left(1+r_{t_d}\right).Math.t_d,.Math.k=1,\mathrm{.Math.}.Math.,n,x=1,\mathrm{.Math.}.Math.,m,& \left(3\right)\end{array}$

where r.sub.[r.sub.y, r.sub.], r.sub.[r.sub., r.sub.], r.sub.[r.sub., r.sub.], r.sub.[r.sub., r.sub.], r.sub.[r.sub., r.sub.], r.sub.td[r.sub.td, r.sup.td], and [(m/n), (m/n),] are utilized to characterize the parametric uncertainties. These trust ranges are also represented as .sub.p[, ], .sub.p[, ], .sub.p[, ], .sub.p[, ], .sub.p[, ], and t.sub.dp[t.sub.d, t.sub.d] for brevity; note that as is a global perturbation on all c.sub.k.sup.p, k=1, . . . , n, [, ] is utilized to represent the trust region of the alignment parameter. Therefore, a set of perturbed models can be characterized by the uncertain model parameters .sub.p, .sub.p, .sub.p, .sub.p, , .sub.p, t.sub.dp.

CD Model Predictive Controller

[0052] For industrial CD-MPC controllers that are applied in paper mills, the following optimization problem is solved:

[00004]$\begin{array}{cc}\underset{.Math..Math.u\left(k\right)}{\min }.Math.\left\{\underset{i=1}{\overset{H_p}{.Math.}}.Math.{\left(y\left(k+i\right)-y_{\mathrm{sp}}\left(k+i\right)\right)}^T.Math.Q_1\left(y\left(k+i\right)-y_{\mathrm{sp}}\left(k+i\right)\right)+\underset{i=0}{\overset{H_u-1}{.Math.}}.Math.\left[.Math..Math.{u\left(k+i\right)}^T.Math.Q_2.Math..Math..Math.u\left(k+i\right)+{\left(u\left(k+i\right)-u_{\mathrm{sp}}\left(k+i\right)\right)}^T.Math.Q_3\left(u\left(k+i\right)-u_{\mathrm{sp}}\left(k+i\right)\right)+{u\left(k+i\right)}^T.Math.Q_4.Math.u\left(k+i\right)\right]\right\},& \left(4\right)\end{array}$

subject to the system dynamics defined in (1) and the constraints as follows:

u(k)bu(k1),(5)

where H.sub.p is the prediction horizon, and H, is the control horizon; y(k)R.sup.m and y.sub.sp.sup.(k)R.sup.m are the predicted output profile and the corresponding reference signal; u(k)R.sup.n and u.sub.sp(k)R.sup.n are the actuator profile and its reference; u(k)(=u(k)u(k1)) is the change in the actuator profile; Q.sub.1 to Q.sub.3 are diagonal weighting matrices; Q.sub.4 is the weighting matrix on the actuator bending/picketing in the following form:

[00005]$\begin{array}{cc}{Q}_4=q_4.Math.S_b^T.Math.S_b,.Math.S_b=\left[\begin{array}{ccccccc}-1& 1& 0& \mathrm{.Math.}& \mathrm{.Math.}& \mathrm{.Math.}& 0\\ 1& -2& 1& & & & \mathrm{.Math.}\\ 0& 1& -2& & & & \mathrm{.Math.}\\ \mathrm{.Math.}& & & & & & \mathrm{.Math.}\\ \mathrm{.Math.}& & & & -2& 1& 0\\ \mathrm{.Math.}& & & & 1& -2& 1\\ 0& \mathrm{.Math.}& \mathrm{.Math.}& \mathrm{.Math.}& 0& 1& -1\end{array}\right],& \left(6\right)\end{array}$

where q.sub.4 is a scalar weight and S.sub.bR.sup.nn is the bending moment matrix. Note that for the actuator profile, the first and second order derivatives are incorporated in the matrix S.sub.b, and thus the bending behavior is penalized in the cost function of CD-MPC. , and b are the constraint matrices (vectors) derived based on the physical limitations of the process.

Temporal Filter

[0053] The traditional output reference trajectory is constructed as a step change, which requires the predicted output profile to track the output target immediately after the dead time of the process. To help the tuning, a temporal filter is utilized to generate the reference trajectory Y.sub.sp(k) based on

Y.sub.sp(k)=F.sub.(y.sub.tgt(k)d.sub.y(k)),(7)

where y.sub.tgt (k) is the output target, and d.sub.y(k)=y.sub.p(k)y(k) is the disturbance estimated based on the process output y.sub.p(k) and predicted output y(k). F.sub. is the time domain implementation of f.sub. (z) based on y.sub.sp(z)=f.sub. (z)I.sub.m(y.sub.tgt (z)d.sub.y(z)) and f.sub. (z) is the temporal filter

[00006]$\begin{array}{cc}{f}_{}\left(z\right)=\frac{\left(1-a_r\right).Math.z^{-t_d}}{1-a_r.Math.z^{-1}},& \left(8\right)\end{array}$

where .sub.r=e.sup.T/; T is the sampling time, and is the continuous-time time constant of the temporal transfer function of the process; I.sub.m represents an m-by-m identity matrix. Note that based on this filter, the aggressiveness of the control signal can be adjusted by the parameter with Q.sub.2 set to a small-valued scalar matrix.

Closed-Loop Transfer Functions

[0054] In order to analyze the properties of the CD-MPC shown in (4), we calculate the equivalent transfer matrices with the unconstrained MPC, and then evaluate robust stability and performance of the corresponding closed-loop system. To aid the analysis, the closed-loop system can be rearranged in FIG. 5. K.sub.r (z) (80) and K.sub. (z) (88) are derived based on the explicit solution of the unconstrained MPC. Note that as the unconstrained MPC is exploited, the resultant closed-loop system is linear.

[0055] In FIG. 5, G.sub.p(z) (82) is the real process that is different from the nominal model G(z), (86) and can be represented in the additive uncertainty form

G.sub.p(z)=G(z)+(z),(9)

in which (z) (84) denotes the model uncertainty.

[0056] Robust stability is analyzed using the small gain theorem. Specifically, given the parametric uncertainties defined in (3), the closed-loop system in FIG. 5 is robustly stable for all G.sub.p(z) if it is nominally stable and

T.sub.ud(z)(z).sub.<1(T.sub.ud(e.sup.i)(e.sup.i))<1,,(10)

T.sub.ud(z)=K.sub.(z)[IG(z)K.sub.(z)].sup.1,(11)

where (A) denotes the maximum singular value of a matrix A; (z) is the model uncertainty in (9) and T.sub.ud (z)C.sup.nm is the sensitivity function from the disturbance profile d(z) to the input profile u(z) at the nominal case.

[0057] Since the performance of CD control is characterized by its capability to suppress the disturbance, the closed-loop transfer function from the disturbance profile d(z) to the output profile y(z) is used to evaluate the performance

T.sub.yd(z)=[IG(z)K.sub.(z)].sup.1,(12)

where T.sub.yd(z)C.sup.mm.

[0058] T.sub.ud (z) and T.sub.yd (z) are employed as the key transfer functions in this paper for evaluating the closed-loop system performance.

[0059] Given the CD-MPC structure shown in FIG. 4, the closed-loop system behavior is affected by the following parameters: the penalty matrices Q.sub.1 to Q.sub.4, the temporal filter parameter , the prediction horizon H.sub.p, and the control horizon H.sub.u.

[0060] Note that in industrial practices, the prediction horizon H.sub.p is normally selected to be 4 times the summation of the time constant and delay; the control horizon H.sub.u=1 is usually utilized in large-scale MPC (can also be increased based on practical situations); Q.sub.1 is often fixed in robust tuning. Consequently, only Q.sub.3, Q.sub.4, and are the tuning parameters in this analysis. As the spatial tuning and temporal tuning can be considered independently, we only consider temporal tuning under temporal uncertainty in this analysis, and the corresponding design objective is to tune the parameters so that: (1) the closed-loop system is robustly temporal stable; (2) the performance indices, namely, settling times, control signal overshoots, meet the end users' specifications.

[0061] The traditional idea is to conduct the tuning based on the sensitivity functions T.sub.ydC.sup.mm and T.sub.udC.sup.nm which, however, have very large dimensions (n can be as high as 300, and m as high as 2000). Such large dimensionality makes the analysis and tuning difficult, and therefore a two-dimensional frequency analysis is introduced in the next section.

II. Two-Dimensional Frequency Analysis

[0062] This section describes the two-dimensional frequency analysis of CD processes and illustrates the decoupling property of the two-dimensional frequencies. The spatially invariant property is one of the most important features of CD processes; this allows us to approximate the plant model G(z) as an rectangular circulant matrix (RCM), and further represent it in the two-dimensional (spatial and temporal) frequency domain. Thus, the multi-variable transfer function G(z) can be simplified into a group of single-input, single-output transfer functions. Denote {hacek over (G)}(z) as the two-dimensional frequency domain representation of G(z), and it can be calculated by

{hacek over (G)}(z)=F.sub.mG(z)F.sub.n.sup.H,(13)

where F.sub.m and F.sub.n are the complex Fourier matrices with dimension m and n, respectively. Without loss of generality, we assume m>n in this analysis. Then, all the frequency information of G(z) is contained in the following transfer functions:

{{hacek over (g)}(v.sub.1,z), . . . ,{hacek over (g)}(v.sub.n,z)}=DIAG(F.sub.mG(z)F.sub.n.sup.H),(14)

where DIAG(A) denotes the operation of getting the following elements of a rectangular matrix AC.sup.mn: DIAG(A).sub.1={A(1, 1), . . . , A(k, k), A(k+1+mn, k+1), . . . , A(m, n)}, where k=n/2 if n is even, or k=(n+1)/2 if n is odd; and

[00007]$\begin{array}{cc}\stackrel{}{g}\left(v_j,z\right)=\frac{\stackrel{}{g}\left(v_j\right).Math.\left(1-a\right).Math.z^{-t_d}}{1-{\mathrm{az}}^{-1}},\mathrm{for}.Math..Math.j=1,\mathrm{.Math.}.Math.,n.& \left(15\right)\end{array}$

Note that based on (15), the spatial frequency response gain of G(z) at v.sub.j is |{hacek over (g)}(v.sub.j)|, where v.sub.j, j=1, . . . , n, are the spatial frequencies with engineering units.

[0063] Given the special structure of CD processes, it can be further shown that the corresponding sensitivity functions T.sub.yd (z), T.sub.ud (z) are also RCMs. Thus, they can also be analyzed in the two-dimensional frequency domain as follows:

{hacek over (T)}.sub.yd(z)=F.sub.mT.sub.yd(z)F.sub.m.sup.H,

{hacek over (T)}.sub.ud(z)=F.sub.nT.sub.ud(z)F.sub.m.sup.H,

{{hacek over (t)}.sup.yd(v.sub.1,z), . . . ,{hacek over (t)}.sub.yd(v.sub.m,z)}=diag({hacek over (T)}.sub.yd(z)),

{{hacek over (t)}.sub.ud(v.sub.1,z) . . . ,{hacek over (t)}.sub.ud(v.sub.n,z)}=DIAG({hacek over (T)}.sub.ud(z)),(16)

where diag(A) represents the diagonal elements of a square matrix A. Therefore controller tuning can be implemented based on the sensitivity functions {hacek over (t)}.sub.ud (v, z) and {hacek over (t)}.sub.yd (v, z) in the two-dimensional frequency domain. Following the idea of CD-MPC tuning for unstructured uncertainty, the tuning task can be separated into two parts: spatial tuning at z=1 (z=e.sup.i, =0) via Q.sub.3 and Q.sub.4, and temporal tuning at v=0 via the parameter (which has the same function as Q.sub.2). In the following, the temporal tuning is explored.

III. Temporal Tuning for CD-MPC

[0064] In this section, we develop the temporal tuning algorithm. In the temporal tuning part, we first calculate a lower bound of the tuning parameter to guarantee robust stability at the spatial frequency v=0. Then, a tuning technique is proposed to find the parameter to achieve the user-specified temporal performance.

[0065] Based on the small gain theory and the two dimensional analysis, the robust stability condition for temporal tuning can be represented as

[00008]$\begin{array}{cc}.Math.{\stackrel{}{t}}_{\mathrm{ud}}\left(0,e^{i.Math..Math.}\right).Math.\left(0,e^{i.Math..Math.}\right).Math.<1,.Math.->.Math.{\stackrel{}{t}}_{\mathrm{ud}}\left(0,e^{i.Math..Math.}\right).Math.<\frac{1}{.Math.\left(0,e^{i.Math..Math.}\right).Math.},,& \left(17\right)\end{array}$

where {hacek over (t)}.sub.ud (0, e.sup.i) and (0, e.sup.i) are T.sub.ud (z) and (z) at spatial frequency v=0. Note that as G.sub.p(z) is an RCM because of the parametric uncertainty, (z)=G.sub.p (z)G(z) is also an RCM based on the property that the RCM is closed on the summation operation. Therefore, (z) can be transformed into two dimensional frequency domain ((v, e.sup.i) denotes its two dimensional form).

[0066] As the extreme system behavior normally happens at the extreme model parameters, the maximum value of |(0, e.sup.i)| can be obtained by solving

[00009]$\begin{array}{cc}\underset{a_p,t_{\mathrm{dp}}}{\max }.Math..Math.\left(0,e^{i.Math..Math.}\right).Math.,.Math..Math.s.t..Math.a_p\left\{\underset{\_}{a},\stackrel{\_}{a}\right\},t_{\mathrm{dp}}\left\{{\underset{\_}{t}}_d,{\stackrel{\_}{t}}_d\right\}.& \left(18\right)\end{array}$

[0067] Then, a bisection search can be utilized for the smallest alpha () that satisfies (17), as increasing a can provide better robustness (i.e., leads to |{hacek over (t)}.sub.ud(0, e.sup.i)|). Then, the obtained is a lower bound for a that guarantees robust stability, and is denoted as * hereafter. Note that only a finite number of temporal frequencies distributed from 0 to the Nyquist frequency may be considered in solving (18) due to the implementation issue.

[0068] In the temporal frequency domain, once robust stability is obtained, we move on to performance tuning. As existing robust tuning approaches normally consider frequency domain performance evaluation, which may not be intuitive to end users in the pulp and paper industry, the time domain performance indices, namely, measurement profile settling times and control signal overshoots, are utilized in this section. The traditional time domain indices are characterized based on the step response of each output measurement, which, however, is not applicable here due to the fact that CD processes, even for single beam cases, have hundreds of output measurements (which we call the measurement profile). Therefore, it is more intuitive and efficient to transform m dimensional output measurements into just one dimension and then to adopt the time domain performance indices. In this analysis, we select the spread of the two times of the standard deviation (2) of the input/output profile as the one dimensional representation which is widely used in the pulp and paper industry and is familiar to end users. Before proceeding, the temporal tuning indices are defined as follows.

[0069] Definition 1 (Overshoot of the 2 spread). The overshoot of a stable 2 spread is its maximum value minus the final value divided by the final value.

[0070] Definition 2 (Settling time of the 2 spread). The settling time of a stable 2 spread is the time required for the spread to reach and stay at its final value.

[0071] As uncertainties are considered on the dynamic parameters a and t.sub.d, which results in a set of process models, an efficient visualization technique needs to be designed to evaluate the possible 2 spreads. This visualization for the output 2 spread can be obtained by solving the following optimization problems:

[00010]$\begin{array}{cc}\underset{a_p,t_{\mathrm{dp}}}{\min }.Math.2.Math.\left(y_i\left(a_p,t_{\mathrm{dp}}\right)\right),\mathrm{for}.Math..Math.i=1,\mathrm{.Math.}.Math.,l.Math..Math.s.t..Math.a_p\left[\underset{\_}{a},\stackrel{\_}{a}\right],t_{\mathrm{dp}}\left[{\underset{\_}{t}}_d,{\stackrel{\_}{t}}_d\right],& \left(19\right)\\ \underset{a_p,t_{\mathrm{dp}}}{\min }.Math..Math.2.Math.\left(y_i\left(a_p,t_{\mathrm{dp}}\right)\right),\mathrm{for}.Math..Math.i=1,\mathrm{.Math.}.Math.,l.Math..Math.s.t..Math.a_p\left[\underset{\_}{a},\stackrel{\_}{a}\right],t_{\mathrm{dp}}\left[{\underset{\_}{t}}_d,{\stackrel{\_}{t}}_d\right],& \left(20\right)\end{array}$

where l is the entire scan time for the process. Note that as the spatial gain matrix G.sub.0 only affects the steady-state value of the 2, it is selected to be the worst spatial model based on a existing method, which produces the worst spatial performance (on the basis of the steady-state values of the 2s).

[0072] Based on the FOPDT structure of h(z), the extreme closed-loop system behavior happens when all the extreme parameters are taken simultaneously, and therefore the above optimization problems can be solved only at a.sub.p{a, } and t.sub.dp{t.sub.d,t.sub.d} to obtain the envelope of the 2 spreads. Similarly, this approach can also be applied on the input profile to achieve visualization. The performance of the visualization method based on optimization problems in (19) and (20) is verified through many industrial examples. To save space, only one typical result is shown in FIG. 6. Note that the model used here is a process in which weight is controlled by the primary autoslice, with the following uncertainties:

a.sub.p[0.5a,1.5a],t.sub.dp[0.5.sub.td,1.5t.sub.d].(21)

[0073] The tuning objective of the temporal part is then formulated as

[00011]$\begin{array}{cc}\underset{}{\min }.Math.T_s\left(2.Math.\left(y\right)\right).Math..Math.s.t..Math.\mathrm{OS}\left(2.Math.\left(u\right)\right){\mathrm{OS}}^{*},& \left(22\right)\end{array}$

where T.sub.s(2 (y)) and OS(2 (u)) are the worst-case (longest) settling time of the output 2 and the worst-case (largest) over-shoot of the input 2 under pre-specified parametric uncertainties, respectively. OS* is the requirement on the worst-case overshoot. In the above tuning objective, the users can specify the maximum allowable overshoot on the input 2, and then solve the optimization problem in (22), yielding the tuning parameter that provides the smallest worst-case settling time of (22) efficiently, we separate the temporal tuning into two parts: (1) tuning in the frequency domain; (2) tuning in the time domain. In the frequency domain tuning, an upper limit on the sensitivity function T.sub.yd(z) at v=0 is first selected, and then we find an a to meet the bound while still guaranteeing robust stability.

[0074] Define .sub.l as the minimum limit of the maximum peak of the sensitivity function (see FIG. 7), the objective here is to adjust so that

|{hacek over (t)}.sub.yd(0,e.sup.i)|<.sub.l,,(23)

is satisfied. Then, the frequency domain tuning can be realized by solving:

[00012]$\begin{array}{cc}{}_f=\underset{k=0,\mathrm{.Math.}.Math.,k_{\mathrm{ma}.Math..Math.x}}{\max }.Math.\left({}_{{}_k}^m\right)& \left(24\right)\\ {}_{{}_k}^m=\min \left({}_{{}_k}\right).Math..Math.s.t..Math.{}_{{}_k}>{}^{*},.Math.{\stackrel{}{t}}_{y.Math..Math.d}\left(0,e^{i.Math..Math.}\right).Math.<{}_l,& \left(25\right)\end{array}$

where .sub.k, k=0, . . . , k.sub.max, are the temporal frequencies distributed from 0 to the Nyquist frequency. .sub.l can be selected according to OS* and the longest time constant of a process model for the given parametric uncertainties, and is normally between 1.2 and 1.8. It is worth noting that the frequency domain tuning may result in an with OS>OS*, and thus the time domain tuning is implemented to fine tune the parameter.

[0075] In the time domain tuning, the parameter is fine tuned to find the optimal a (denoted as .sup.t) that provides the smallest worst-case settling time while still satisfying the requirement on OS. Since OS(2(u)) is a monotonically decreasing function of a, the bisection search can be utilized to find the .sup.t based on the closed-loop simulation with the following procedure. If the OS with .sub.f is larger than the specification, the search is implemented in the region [.sub.f, .sub.u], where .sub.u is the upper bound for with .sub.l=1.2, and .sub.f is obtained in the frequency tuning. Otherwise, the search is implemented in the region [.sub.i, .sub.f], where .sub.l is the lower bound for a with .sub.l=1.8. To sum up, the temporal tuning part is implemented using the following algorithm:

Algorithm 1 Temporal Tuning

[0076] 1: Input OS*;

[0077] 2: Input the parametric uncertainty specifications [, ] and [t.sub.d, t.sub.d];

[0078] 3: Calculate the lower bound for robust stability *;

[0079] 4: Implement the frequency tuning to achieve .sub.f;

[0080] 5: if OS(.sub.f)>OS* then

[0081] 6: Search in [.sub.f, .sub.u] for .sup.t;

[0082] 7: else

[0083] 8: Search in [.sub.l, .sub.f] for .sup.t;

[0084] 9: end if

[0085] 10: End

[0086] Implementing the tuning procedure in Algorithm 1, the automated temporal tuning costs about 6 seconds on a computer with i5 CORE and 6G RAM.

[0087] Remark 1.

[0088] To save more computation time, we can also stop the algorithm if OS(.sub.f)OS*. This can avoid the tuning in the time domain once the requirement is satisfied in the frequency tuning.

[0089] Remark 2.

[0090] The visualization technique can also be extended to characterize the envelope of the 2 spreads for the best spatial performance under the given parametric uncertainties, based on the best spatial model (with the best spatial frequency property and can be achieved similarly as the worst spatial model). Combining the best and worst envelope spreads results in a more detailed visualization result, which is shown in FIG. 9. Here the lower broken line represents the best 2 spread (with smallest steady state value and settling time), while the upper broken line is for the worst 2 spread (with largest steady state value and settling time); the area in between indicates all possible 2 spreads for the given parametric uncertainties; and the solid line represents the nominal case.

[0091] As illustrated in FIG. 8, automatic temporal tuning of the CD-MPC includes of the following 7 steps:

[0092] (1) Input the nominal temporal model parameters and the corresponding parametric uncertainty specifications. The nominal temporal model parameters refers to the parameters of the temporal mathematical model of the process obtained from the model identification procedure. The parametric uncertainty refers to the trust ranges around the nominal temporal model parameters, which characterize the possible mismatch between the identified model and the real process.

[0093] (2) Specify a 2 sigma (2) overshoot limit for the CD actuator profile. The worst-case 2 overshoot refers to the largest overshoot of all possible 2 spreads given the parametric uncertainty specifications. The definition of 2 overshoot is set forth in Definition 1.

[0094] (3) Develop a robust temporal stability condition based on the parametric uncertainty specifications. The robust temporal stability condition refers to the requirement on the temporal frequency response of the closed-loop system, which guarantees the robust stability under the parametric uncertainty. Given the robust temporal condition, the frequency responses of a large number of uncertain models based on the parametric uncertainty need to be tested. Based on the analysis in the present invention, the condition can be checked based on only 4 critical frequency responses (as illustrated in equation (8)). This saves computational time and guarantees the design/tuning procedure of the digital controller can be accomplished within the required amount of time.

[0095] (4) Calculate a lower bound for the profile trajectory tuning parameter a that guarantees the robust stability. The profile trajectory tuning parameter refers to a parameter () in the temporal filter (as defined by equation (8)) that can affect the aggressiveness of the reference trajectory (as defined by equation (7)). The a lower bound for (the profile trajectory tuning parameter) refers to the smallest value of the tuning parameter that guarantees the robust stability condition.

[0096] (5) Use a frequency domain technique to further reduce the search range for a. The search range for refers to the range of the tuning parameter within which the optimal tuning parameter may exist.

[0097] (6) Perform an intelligent search for that minimizes measurement 2 settling time without exceeding the actuator 2 overshoot limit. The worst-case measurement 2 settling time refers to the longest settling time of all possible 2a spreads given the parametric uncertainty specifications (2 setting time is described in Definition 2). Based on the definition of worst-case overshoot/settling time, the 2 overshoot/settling time of a large number of 2 spreads given the parametric uncertainty need to be calculated to find the worst-case 2 overshoot/settling time, which requires a lot of computational time (as illustrated in equations (19)-(20)). Based on analysis in the invention, the worst-case 2 overshoot/settling time can be calculated only based on 4 extreme models, which significantly reduces the computational time; this guarantees the design/tuning procedure of the digital controller can be accomplished within the required amount of time. The intelligent search refers to a two-step frequency- and time-domain combined search method that can find the optimal tuning parameter within a small amount of time. Given the performance requirement, the intelligent search first finds a potential tuning parameter based on equations (24)-(25) with the temporal frequency analysis, and then fine tunes the parameter for optimality with line 5-10 in Algorithm 1. The proposed two-step search in the invention guarantees the optimality and the computational efficiency of the intelligent search, which makes the tuning/design procedure of the digital controller can be accomplished within the required amount of time.

[0098] (7) Output the tuning parameter .

IV. Simulation Results with Real-Time Simulator

[0099] We apply the proposed automated robust temporal tuning algorithms to a process model obtained from a paper mill in Canada. This is a single beam process in which the dry weight of the output profile is controlled by the primary autoslice.

[0100] In this test, the parametric uncertainty specifications are given as

[00013]$\begin{array}{cc}{}_p\left[0.75.Math.,1.25.Math.\right],{}_p\left[0.75.Math.,1.25.Math.\right],.Math.{}_p\left[0.75.Math.,1.25.Math.\right],{}_p\left[0.75.Math.,1.125.Math.\right],.Math.a_p\left[0.5.Math.a,1.5.Math.a\right],t_{d_p}\left[0.5.Math.t_d,1.5.Math.t_d\right],.Math.\mathrm{.Math.}\left[-0.6.Math..Math.\frac{m}{n},0.6.Math..Math.\frac{m}{n}\right].& \left(26\right)\end{array}$

[0101] Then, the requirement of the robust temporal tuning of CD-MPC is to guarantee that the time-domain performance indices are satisfied, i.e., a fast settling time with the worst-case overshoot of the actuator 2 less than OS*=20%.

[0102] The overall tuning costs 6 seconds on a computer with Intel i5 core and 6G memory. The temporal tuning provides =1.6821 (see FIG. 10, which shows that the requirement OS(2(u))<OS* is satisfied), and the spatial parameters are chosen based on the existing robust spatial tuning method. We then apply the obtained tuning results to the Honeywell real time MPC+simulator environment and compare the obtained results with that of the existing controller. Note that the real process model G.sub.p is selected satisfying (26). The performance is shown in FIGS. 11(a)-11(b). FIGS. 11(a) and 11(b) illustrate the 2 spreads for the measurement profile and actuator profile, in which we can see that the proposed tuning approach can achieve much better temporal performance in both the measurement and actuator profiles. More specifically, it is observed that the measurement 2 spread even becomes larger than its initial value for the existing controller under the pre-specified parametric uncertainty when it approaches the steady-state, which illustrates the robustness of the proposed method.

[0103] In some embodiments, various functions described above are implemented or supported by a computer program that is formed from computer readable program code and that is embodied in a computer readable medium. The phrase computer readable program code includes any type of computer code, including source code, object code, and executable code. The phrase computer readable medium includes any type of medium capable of being accessed by a computer, such as read only memory (ROM), random access memory (RAM), a hard disk drive, a compact disc (CD), a digital video disc (DVD), or any other type of memory. A non-transitory computer readable medium excludes wired, wireless, optical, or other communication links that transport transitory electrical or other signals. A non-transitory computer readable medium includes media where data can be permanently stored and media where data can be stored and later overwritten, such as a rewritable optical disc or an erasable memory device.