Method of designing model predictive control for cross directional flat sheet manufacturing processes to guarantee spatial robustness and to prevent actuator picketing

Abstract

Automated parameter tuning techniques for cross-directional model predictive control for paper-making under user-specified parametric uncertainties to reduce variability of the actuator and measurement profiles in the spatial domain is proposed. Decoupling properties of the spatial and temporal frequency components permit separate controller design and parameter tuning. CD-MPC design that explicitly accounts for parametric model uncertainty while finding MPC cost function weighing matrices that prevent actuator picketing and guarantee robust stability of the spatial CD profile. Picketing refers to periodic variation patterns in the actuator array. The inventive technique includes: (i) determining the worst case cutoff frequency of all process models, given parametric uncertainty, (ii) designing a weighing matrix to penalize high frequency actuator variability based on the process model and worst case cutoff frequency, and (iii) finding a multiplier for the spatial frequency weighted actuator variability term in the MPC cost function that assures robust spatial stability.

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 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 to the cross-directional process, wherein the MPC, in response to signals that are indicative of the properties of the properties data, provides signals to the at least one set of actuator arrays to vary properties of the material, wherein the MPC includes means for spatially tuning the MPC which comprises: (a) means for inputting nominal spatial model parameters and corresponding parametric uncertainty specifications; (b) means for calculating worst-case cut-off frequency v.sub.w for all possible process models given the parametric uncertainty specification wherein the worst-case cut-off frequency comprises the smallest cut-off frequency of all possible process models given the parametric uncertainty specifications; (c) means for designing a weighting matrix S.sub.b to penalize high frequency actuator variability based on the process model and worst-case cut-off frequency; (d) means for developing a robust spatial condition based on the parametric uncertainty specifications; (e) means for adjusting a multiplier of S.sub.b in a MPC cost function to satisfy the robust stability condition wherein S.sub.b is a weighting matrix in penalty term of a MPC cost function which controls actuator bending and picketing behavior; and (f) means for outputting the weighting matrix S.sub.b and its multiplier.

2. The system of claim 1 wherein the S.sub.b is a spatial filter whose spatial frequency response is a mirror of an actuator frequency response.

3. The system of claim 1 wherein the parametric uncertainty comprises trust ranges around the nominal spatial model parameters which characterize possible mismatch between an identified model and an actual process.

4. The system of claim 1 wherein high frequency actuator variability comprises actuator bending and picketing behavior that includes frequency components beyond the worst-case cut-off frequency.

5. The system of claim 1 wherein the nominal spatial model parameters comprise parameters of a mathematical model of the process obtained from a bump test and model identification procedure.

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

7. 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 spatial model parameters and corresponding parametric uncertainty specifications; (ii) calculating worst-case cut-off frequency v.sub.w for all possible process models given the parametric uncertainty specification; (iii) designing a weighting matrix S.sub.b to penalize high frequency actuator variability based on the process model and worst-case cut-off frequency wherein the worst-case cut-off frequency comprises the smallest cut-off frequency of all possible process models given the parametric uncertainty specifications; (iv) developing a robust spatial condition based on the parametric uncertainty specifications; (v) adjusting a multiplier of S.sub.b in a MPC cost function to satisfy the robust stability condition wherein S.sub.b is a weighting matrix in penalty term of a MPC cost function which controls actuator bending and picketing behavior; and (vi) outputting the weighting matrix S.sub.b and its multiplier; (b) inputting the tuning parameters into the MPC; and (c) controlling the multiple-array CD process with the MPC.

8. The system of claim 7 wherein the S.sub.b is a spatial filter whose spatial frequency response is a mirror of an actuator frequency response.

9. The system of claim 7 wherein the parametric uncertainty comprises trust ranges around the nominal spatial model parameters which characterize possible mismatch between an identified model and an actual process.

10. The system of claim 7 wherein high frequency actuator variability comprises actuator bending and picketing behavior that includes frequency components beyond the worst-case cut-off frequency.

11. The system of claim 7 wherein the nominal spatial model parameters comprise parameters of a mathematical model of the process obtained from a bump test and model identification procedure.

12. The system of claim 7 wherein the spatially-distributed sheetmaking process is a paper-making process.

13. 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 spatial model parameters and corresponding parametric uncertainty specifications; (b) calculating worst-case cut-off frequency v.sub.w for all possible process models given the parametric uncertainty specifications wherein the worst-case cut-off frequency comprises the smallest cut-off frequency of all possible process models given the parametric uncertainty specifications; (c) designing a weighting matrix Sb to penalize high frequency actuator variability based on the process model and worst-case cut-off frequency; (d) developing a robust spatial condition based on the parametric uncertainty specifications; (e) adjusting a multiplier of St in a MPC cost function to satisfy the robust stability condition wherein S.sub.b is a weighting matrix in penalty term of a MPC cost function which controls actuator bending and picketing behavior; (f) outputting the weighting matrix S.sub.b and its multiplier; (g) inputting tuning parameters into the MPC; and (h) controlling the cross-directional process with the MPC.

14. The computer readable medium claim 13 wherein the S.sub.b is a spatial filter whose spatial frequency response is a mirror of an actuator frequency response.

15. The computer readable medium of claim 13 wherein the parametric uncertainty comprises trust ranges around the nominal spatial model parameters which characterize possible mismatch between an identified model and an actual process.

16. The system of claim 1 wherein the at least one set of actuator arrays comprise a single array comprising of a plurality of manipulated actuators that are arranged in the CD and wherein the means for measuring and acquiring properties data about the properties of the material comprises a corresponding single controlled measurement array.

17. The system of claim 7 wherein the at least one manipulated actuator array comprises a single array comprising of a plurality of manipulated actuators that are arranged in the CD and wherein the at least one controlled measurement array comprises a corresponding single controlled measurement array.

18. The computer readable medium of claim 13 wherein the manipulated actuator array comprises a single array comprising of a plurality of manipulated actuators that are arranged in the CD and wherein the at least one controlled measurement array comprises a corresponding single controlled measurement array.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) FIGS. 1, 2, and 3 are schematic illustrations of a papermaking system;

(2) FIG. 4 is a block diagram of a closed-loop CD-MPC system;

(3) FIG. 5 is a block diagram of a closed-loop CD-MPC system;

(4) FIG. 6 is the spatial frequency response of the existing S.sub.b;

(5) FIG. 7 illustrates the definition of the cut-off frequency v.sub.c;

(6) FIG. 8 illustrates the spatial frequency response of the new S.sub.b designed by Algorithm 2;

(7) FIG. 9 illustrates the performance comparison of different S.sub.b (existing one and the proposed one) based on the sensitivity function T.sub.yd (z) at steady state (z=1);

(8) FIGS. 10(a) to 10(e) are graphs that illustrate the relationship between the spatial parameters and v.sub.c (note that .sub.p, .sub.p, .sub.p, .sub.p, increased from 50% to 150% of the nominal values, and is increased from m/n to m/n;

(9) FIGS. 11(a) and 11(b) illustrate the sensitivity function T.sub.yd (z) and control sensitivity T.sub.ud (z), respectively, in the spatial frequency domain with (z=1);

(10) FIGS. 12(a) and 12(b) illustrate the sensitivity function T.sub.yd (z) and control sensitivity T.sub.ud (z), respectively, in the spatial frequency domain with (z=1);

(11) FIGS. 13(a) to 13(e) illustrate the relationship between the spatial parameters and the spatial gain | (v.sub.j, 1)|;

(12) FIG. 14 is the performance illustration of the maximum singular values calculated based on the optimized problem in (31);

(13) FIG. 15 is the flow chart of the automated spatial tuning;

(14) FIG. 16 illustrates spatial tuning result for the proposed approach (note that the upper line is for the |1/ (v, 1)|; and the lower line is for the |t.sub.ud (v, 1)|; and

(15) FIGS. 17(a) to 17(b) are performance illustrations for the proposed robust tuning method.

DETAILED DESCRIPTION OF PREFERRED EMBODIMENTS

(16) 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.

(17) 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.

(18) 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.

(19) 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.

(20) 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. 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.

(21) I. CD-MPC Structure

(22) 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.a(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.

(23) Nominal Model

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

(25) $\begin{array}{cc}G\left(z\right)=G_{0}h\left(z\right),h\left(z\right)=\frac{\left(1-a\right)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 and t.sub.d are the time constant and time delay in the discretized version.

(26) The spatial gain matrix G.sub.0 has the parameterized structure as shown below:

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

(28) Since model-plant mismatch is unavoidable in process operation and identification, model uncertainties are considered in this invention. 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

(29) $\begin{array}{cc}{G}_p\left(z\right)=G_0^p{h}_p\left(z\right),G_0^p=\left[g_{1p},\mathrm{.Math.},g_{\mathrm{np}}\right],g_{\mathrm{kp}}=f\left(x,{}_p,{}_p,{}_p,{}_p,c_k^p\right),{}_p=\left(1+r_{}\right),{}_p=\left(1+r_{}\right),{}_p=\left(1+r_{}\right),{}_p=\left(1+r_{}\right),c_k^p=c_k+\mathrm{.Math.},h_p\left(z\right)=\frac{\left(1-a_p\right)z^{-t_{\mathrm{dp}}}}{1-a_p{z}^{-1}},a_p=\left(1+r_a\right)a,t_{\mathrm{dp}}=\left(1+r_{t_d}\right)t_d,k=1,\mathrm{.Math.},n,x=1,\mathrm{.Math.},m,& \left(3\right)\end{array}$
where 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., r.sub.], r.sub.td[r.sub.td, r.sub.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

(30) For industrial CD-MPC controllers that are applied in paper mills, the following optimization problem is solved:

(31) $\begin{array}{cc}\underset{u\left(k\right)}{\min }\left\{\underset{i=1}{\overset{H_p}{.Math.}}{\left(y\left(k+i\right)-y_{\mathrm{sp}}\left(k+i\right)\right)}^T{Q}_1\left(y\left(k+i\right)-y_{\mathrm{sp}}\left(k+i\right)\right)+\underset{i=0}{\overset{H_u-1}{.Math.}}\left[{u\left(k+i\right)}^T{Q}_{2}u\left(k+i\right)+{\left(u\left(k+i\right)-u_{\mathrm{sp}}\left(k+i\right)\right)}^T{Q}_3\left(u\left(k+i\right)-u_{\mathrm{sp}}\left(k+i\right)\right)+{u\left(k+i\right)}^T{Q}_{4}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 Hu is the control horizon; y(k)R.sup.m and y.sub.sp(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:

(32) $\begin{array}{cc}{Q}_4=q_4{S}_b^T{S}_b,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

(33) 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. For illustrative purposes, a known temporal filter is utilized to generate the reference trajectory Y.sub.sp(k) based on
Y.sub.sp(k)=F.sub.(y.sub.trg(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

(34) $\begin{array}{cc}{f}_{}\left(z\right)=\frac{\left(1-a_r\right)z^{-t_d}}{1-a_r{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

(35) 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.

(36) 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.

(37) 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.fwdarw.(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.

(38) 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.

(39) 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.

(40) 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.

(41) 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 a are the tuning parameters in this analysis. As the spatial tuning and temporal tuning can be considered independently, we only consider spatial tuning in the work, and corresponding design objective is to tune the parameters so that: (1) the closed-loop system is robustly spatial stable and (2) the variability of the steady-state actuator and measurement profiles is minimized.

(42) 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.

(43) II. Two-Dimensional Frequency Analysis

(44) 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 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).sup.1={A(1, 1), . . . , A(k, k), A(k+1+mn, k+1), . . . , A(m, n)}, where k=n/2 if in is even, or k=(n+1)/2 if n is odd; and

(45) $\begin{array}{cc}\stackrel{}{g}\left(v_j,z\right)=\frac{\stackrel{}{g}\left(v_j\right)\left(1-a\right)z^{-t_d}}{1-{\mathrm{az}}^{-1}},\mathrm{for}j=1,\mathrm{.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.
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:
.sub.yd(Z)=F.sub.mT.sub.yd(z)F.sub.m.sup.H,
.sub.ud(Z)=F.sub.nT.sub.ud(z)F.sub.m.sup.H,
{.sub.yd(v.sub.1,z), . . . ,.sub.yd(v.sub.m,z)}=diag(.sub.yd(z)),
{.sub.ud(v.sub.1,z), . . . ,.sub.ud(v.sub.n,z)}=DIAG(.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 .sub.ud (v z) and .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 spatial tuning is explored.
III. Spatial Tuning for CD-MPC

(46) In this section, we develop the spatial tuning algorithm. In spatial tuning, the objective is to tune matrices Q.sub.3=q.sub.3I.sub.n, Q.sub.4=q.sub.4S.sub.b.sup.TS.sub.b for robust stability and for improving consistency of the actuator and measurement profiles. First a new S.sub.b matrix is designed to improve the spatial performance and then a tuning algorithm is proposed to automatically adjust the weights q.sub.3 and q.sub.4.

(47) The New S.sub.b to Shape the Spectra of the Actuator Profile

(48) In this subsection, we first perform a spatial frequency analysis on the existing S.sub.b and then propose a new S.sub.b to improve the performance. Note that the spatial frequency representation of a signal can be calculated by the Discrete Fourier Transform (DFT), i.e., given a signal yC.sup.m, its spatial frequency representation can be obtained by
{hacek over (y)}=F.sub.m.Math.y,(17)
where F.sub.m is the Fourier matrix. To further analyze the signal in the spatial frequency domain, we need the following useful lemma.
Lemma 1. Given spatial signals yC.sup.m and uC.sup.n, and assuming NC.sup.mn is an RCM, if y=N.Math.u, then
{hacek over (y)}=.Math.{hacek over (u)},(18)
where =F.sub.mNF.sub.n.sup.H.
Note that if N is a constant RCM, in Lemma 1 denotes its spatial frequency domain representation.

(49) In the cost function of CD-MPC, the cost term associated with Q.sub.4 is
u(k+i).sup.TQ.sub.4u(k+i)=(S.sub.bu(k+i)).sup.Tq.sub.4(S.sub.bu(k+i)),
where the matrix S.sub.b can be approximated as an RCM; its spatial frequency domain representation can be calculated as
{b(v.sub.1), . . . ,b(v.sub.n)}=diag(F.sub.nS.sub.bF.sub.n.sup.H).(19)

(50) FIG. 6 illustrates the spatial frequency response (gain) of the existing S.sub.b (namely, |b(v.sub.j)|), in which we can see that S.sub.b looks like a high pass filter in the spatial frequency domain. Therefore, when Q.sub.4 penalizes the actuator picketing, based on Lemma 1, it actually puts more weights on the high frequency components of the actuator profile, which works like a low pass filter to suppress the high frequency components in the actuator profile u. The reason that an actuator profile with high spatial frequency (namely, actuator picketing) is not desired is that CD processes normally have very small gains at high spatial frequencies (also known as ill-conditioning), and an actuator profile with high spatial frequency could lead to instability of CD processes. However, the spatial frequency response of different CD processes behave differently, and thus, it is beneficial to design a specific S.sub.b for a specific CD process.

(51) To achieve better robustness and performance, the new S.sub.b here is designed to be a high pass spatial filter whose stop band can be adjusted according to different CD processes. To achieve this objective, the definition of cut-off frequency, which is used to characterize the spatial frequency response gain of a CD process, is introduced as follows.

(52) Definition 1.

(53) Assume the maximum gain of the spatial frequency response of a CD process is g.sub.max. The cut-off frequency vc is the spatial frequency at which the following holds
g.sub.v.sub.c=r.Math.g.sub.max,r(0,0.5].(20)

(54) FIG. 7 illustrates the above definition using one of the typical single beam CD processes (note that we use r=0.1 in this analysis).

(55) As the process is almost uncontrollable above the cut-off frequency, the actuator profile with the spatial frequency components that is higher than v.sub.c is not desirable. Therefore, the stop band of the new S.sub.b can be selected based on v.sub.c of a given CD process, which allows the new Sb to put more weights on the high spatial frequency components (>v.sub.c) in the input profile while putting small or zero weights on the spatial frequency that is smaller than v.sub.c. The first approach is to design the new S.sub.b as a high pass Finite Impulsive Response (FIR) filter with the stop band that equals v.sub.c of a given process model. This can be realized by the following algorithm. Algorithm 1 (Design S.sub.b as a high pass FIR filter with stop band that equal v.sub.c) 1: Input G.sub.0; 2: Calculate v.sub.c based on G.sub.0; 3: Design a low pass filter based on the sinc function and the window function with the stop band that equals v.sub.c; 4: Conduct the inverse Fourier transform to obtain the filter in the spatial domain; 5: Construct the desired high pass filter S.sub.b via spectral inversion of the filter obtained in line 4; 6: End

(56) Executing Algorithm 1 yields the desired S.sub.b. However, there exist two drawbacks for this method: (1) the stop band of S.sub.b (which we denote as v.sub.b hereafter) may not equal v.sub.c exactly, because the length of the filter is limited by the size of the penalty matrix; (2) even if v.sub.b=v.sub.c can be achieved accurately, the transition band of the filter S.sub.b still needs to be optimized based on the spatial frequency response of G(z).

(57) Considering the aforementioned drawbacks, it is desirable to design S.sub.b with a spectrum that is the mirror of the system's spatial frequency response with respect to the cut-off frequency. To achieve this target, we need to introduce the real valued Fourier matrix, which is defined below.

(58) Definition 2.

(59) The real valued Fourier matrix F.sup.rR.sup.nn defined as

(60) $\begin{array}{cc}{F}_n^r\left(j,k\right)=\{\begin{array}{cc}\sqrt{\frac{1}{n}},& j=1,\\ \sqrt{\frac{2}{n}}.Math.\sin \left[\left(k-1\right)v_j\right],& j=2,\mathrm{.Math.},q,\\ \sqrt{\frac{2}{n}}.Math.\cos \left[\left(k-1\right)v_j\right],& j=q+1,\mathrm{.Math.},n\end{array},& \left(21\right)\end{array}$
where q=(n+1)/2 if n is odd, or q=n/2 if n is even; v.sub.j=2(j1)/n is the spatial frequency.

(61) The basic idea to generate the desired S.sub.b with real valued Fourier matrix is that we first obtain the mirrored frequency response of a given process model by numerical methods, and then perform the inverse Fourier transform using the real valued Fourier matrix, resulting in a real valued Ss with the expected spatial frequency response.

(62) Although the real valued Fourier matrix helps build a real valued S.sub.b, we still need to ensure that the obtained S.sub.b has the circulant property to retain the two-dimensional feature of the closed-loop system. Thus, we need to establish conditions on the pre-specified spectrum content to guarantee the circulant property of S.sub.b.

(63) Lemma 2.

(64) Define k.sub.m(j), j=1, . . . , n, as the pre-specified spatial frequency response gain of S.sub.b. If k.sub.m(j)R, j, and k.sub.m(j)=k.sub.m(nj+2) (j=2, . . . , (n+1)/2, if n is odd; or j=2, . . . , n/2, if n is even), then the S.sub.b obtained via the inverse real valued Fourier transform is a symmetric real valued circulant matrix.

(65) Proof.

(66) As S.sub.b is a square RCM, the spatial frequency response gains of S.sub.b equal its singular values. Therefore, as the real conjugate even property is guaranteed by k.sub.m(j)R, j, and k.sub.m(j)=k.sub.m(nj+2) (j=2, . . . , (n+1)/2, if n is odd; or j=2, . . . , n/2, if n is even), the obtained S.sub.b is a real valued symmetric circulant matrix.

(67) The method using the real valued Fourier matrix approach is then summarized in the following algorithm.

(68) Algorithm 2 (Design S.sub.b based on the mirror of the spectrum of the system model).

(69) 1: Input G.sub.0; 2: Calculate vc based on G.sub.0; 3: Obtain the desired mirrored spectrum by numerical methods; 4: Process the obtained spectrum content based on Lemma 2; 5: Conduct the real valued inverse Fourier transform to obtain the weighting matrix S.sub.b in the spatial domain; 6: End

(70) FIG. 8 shows the spatial frequency response of the S.sub.b designed using Algorithm 2 for a typical single beam CD process. The closed-loop performance of the new S.sub.b is illustrated in the next section.

(71) Performance Analysis of the New S.sub.b

(72) As the new S.sub.b is designed based on the frequency property of G(z), it is expected to provide better performance compared with that of the old S.sub.b. In this section, the performance improvement is illustrated based on the sensitivity analysis.

(73) To analyze the performance of the new S.sub.b, we compare the sensitivity function T.sub.yd (z) of the new S.sub.b with that of the old S.sub.b considering similar robustness (the peak value of T.sub.ud is the same). The details are shown in FIG. 9. Note that the result is only shown at the steady state (z=1) for clear display. We can see that, on the one hand, the sensitivity function with the new S.sub.b has smaller gains at low frequencies, indicating that the disturbance at low frequencies (smaller than the cut-off frequency v.sub.c) can be better handled. On the other hand, for the frequency components that are higher than the cut-off v.sub.c, the sensitivity function with the new S.sub.b has larger gains, which prevents the actuator profile from chasing the uncontrollable high frequency components of the disturbance.

(74) It is worth mentioning that if there is no model mismatch, the new S.sub.b can be designed with v.sub.b=v.sub.c, and should achieve performance improvement based on the above analysis. However, in the presence of parametric uncertainty, it is better to design v.sub.b to be robust over the model mismatch, specifically, to design v.sub.b according to the worst spatial model. To achieve this target, the prediction of the worst spatial model, given user-specified parametric uncertainties, is introduced in the next section.

(75) Worst Spatial Model Prediction Under Parametric Uncertainty

(76) In this section, a technique to find the worst spatial model under user-specified parametric uncertainties is proposed, which can help us design a robust S.sub.b that takes the pre-specified model mismatch into account. Besides, it can also be utilized to develop a robust performance visualization method to help end users predict the possible control performance. First, we propose the definition for the worst spatial model as follows.

(77) Definition 3.

(78) Given the parametric uncertainty specifications
.sub.p[,],.sub.p[,],.sub.p[,,],.sub.p[,], and [,],
the worst spatial model refers to the spatial model that has the smallest cut-off frequency among all the possible models due to parametric uncertainty.

(79) Based on Definition 3, the cut-off frequency of the worst spatial model under parametric uncertainty can be obtained by

(80) $\begin{array}{cc}\underset{{}_p,{}_p,{}_p,{}_p,\mathrm{.Math.}}{\min }{v}_c\left(G_p\left(z\right)\right)s.t.{}_p\left[\underset{\_}{},\stackrel{\_}{}\right],{}_p\left[\underset{\_}{},\stackrel{\_}{}\right],{}_p\left[\underset{\_}{},\stackrel{\_}{}\right],{}_p\left[\underset{\_}{},\stackrel{\_}{}\right],\mathrm{.Math.}\left[\underset{\_}{\mathrm{.Math.}},\stackrel{\_}{\mathrm{.Math.}}\right].& \left(22\right)\end{array}$
Lemma 3.

(81) The cut-oft frequency vc for a given process model G.sub.p(z) is independent of the process gain parameter .sub.p.

(82) Proof. As G.sub.p(z) can be approximated as an RCM, v.sub.c can be obtained by f(x, .sub.p, .sub.p, .sub.p, .sub.p, c.sub.k.sup.p) for any k. Therefore, as .sub.p is a multiplication term in f(x, .sub.p, .sub.p, .sub.p, c.sub.k.sup.p), it does not affect the normalized shape of the spectrum of the system model G.sub.p(z) based on the property of the Fourier transform, and thus is independent of v.sub.c.

(83) The approach utilized here is to hold all the spatial parameters except one to analyze how the change of this parameter affects the cut-off frequency v.sub.c, and apply similar analysis to the parameters .sub.p, .sub.p, .sub.p, .sub.p, one by one. The results of the aforementioned method are shown in FIG. 10 for one CD process. It is observed that the cut-off frequency v.sub.c is changing monotonically with respect to .sub.p, .sub.p, .sub.p, and is not affected by .sub.p. As v.sub.c(G.sub.p(z)) is independent of .sub.p, and is changing monotonically with respect to each of .sub.p, .sub.p, .sub.p, , the smallest cut-off frequency can be obtained by solving the following optimization problem:

(84) 0$\begin{array}{cc}\underset{{}_p,{}_p,{}_p,\mathrm{.Math.}}{\min }{v}_c\left(G_p\left(z\right)\right)s.t.{}_p\left\{\underset{\_}{},\stackrel{\_}{}\right\},{}_p\left\{\underset{\_}{},\stackrel{\_}{}\right\},{}_p\left\{\underset{\_}{},\stackrel{\_}{}\right\},\mathrm{.Math.}\left[\underset{\_}{\mathrm{.Math.}},\stackrel{\_}{\mathrm{.Math.}}\right].& \left(23\right)\end{array}$

(85) The obtained worst cut-off frequency is denoted as v.sub.w, and will be utilized in the robust spatial tuning of CD-MPC.

(86) The Automated Spatial Tuning

(87) The spatial tuning of robust CD-MPC contains two parts: (1) the selection of the stop band v.sub.b of the matrix S.sub.b; (2) the selection of q3 and q4. Note that following the idea in existing CD-MPC tuning methods, the weight q3 is selected to be the same as q4 to simplify the robust tuning problem (this new weight is denoted as q.sub.f).

(88) In order to achieve the tuning objective, the analysis based on the sensitivity functions is utilized to provide some tuning guidelines, and then the automated tuning approach is proposed. FIGS. 11(a) and 11(b) illustrate the effect of changing v.sub.b in the sensitivity functions, and it is observed that a larger v.sub.b gives better robustness, while a smaller v.sub.b gives better performance. If there is no model mismatch, we can select v.sub.b=v.sub.c based on the spatial frequency response of a given process model, the performance of which is also illustrated previously in FIGS. 6 to 9. As we can predict the worst spatial model that has the smallest cut-off frequency (v.sub.w) based on the approach in the previous subsection, the design of v.sub.b=v.sub.w directly takes care of the given parametric uncertainties, so that a good robustness property can be obtained.

(89) Similar analysis can also be applied to the parameter q.sub.f, and the results are shown in FIGS. 12(a) and 12(b). We can observe that a larger q.sub.f gives better robustness, while a smaller q.sub.f gives better performance. As v.sub.b=v.sub.w is selected, the tuning objective of q.sub.f is to achieve robust stability in the spatial domain with z=1. As robust stability and performance are characterized based on sensitivity functions, it is useful to derive the sensitivity functions in the spatial frequency domain with z=1 (steady state =0), which are obtained as follows:
T.sub.yd(z=1)=(I+G.sub.0(Q.sub.3+Q.sub.4).sup.1G.sub.m.sup.TQ.sub.1).sup.1,(24)
T.sub.ud(z=1)=(Q.sub.3+Q.sub.4+G.sub.m.sup.TQ.sub.1G.sub.0).sup.1G.sub.m.sup.TQ.sub.1,(25)
where G.sub.m=.sub.k=1.sup.Hp-td.sub.l=1.sup.k.sup.l-1.

(90) The tuning objective is then to find the smallest q.sub.f so that the robust stability condition in the spatial frequency domain is satisfied, because it can provide the best performance while still guaranteeing robust stability according to the tuning guidelines in FIGS. 12(a) and 12(b). Based on the robust stability condition in (11) and the two-dimensional frequency analysis, robust stability at z=1 can be achieved if the following condition is satisfied:
(T.sub.ud(1)(1))<1.(26)

(91) However, in the presence of parametric uncertainty, the singular values of (1) are hard to compute directly, and therefore the analysis on the uncertain term is carried out first. From (9), the uncertain term (1) can be represented as
(1)=G.sub.p(1)G(1).(27)

(92) As uncertainties are specified on the model parameters, G.sub.p(1) is also an RCM based on the method that generates the spatial gain matrix in (2). Then, (1) is also an RCM according to the property that the RCM is closed on the summation operation, and therefore it can also be transformed into two-dimensional frequencies
{hacek over ()}(1)=F.sub.m(1)F.sub.n.sup.H.
{(v.sub.1,1), . . . ,(v.sub.n,1)}=diag({hacek over ()}(1)).(28)

(93) Then, since the singular values equal the gains along the spatial frequencies, equation (26) can then be represented in the two-dimensional frequency domain as

(94) $\begin{array}{cc}\underset{j}{\max }\left(.Math.{\stackrel{}{t}}_{\mathrm{ud}}\left(v_j,1\right)\left(v_j,1\right).Math.\right)<1->.Math.{\stackrel{}{t}}_{\mathrm{ud}}\left(v_j,1\right).Math.<\frac{1}{.Math.\left(v_j,1\right).Math.},j,& \left(29\right)\end{array}$
where .sub.ud (v.sub.j, 1) is the sensitivity function T.sub.ud(z) with z=1.

(95) Existing tuning approaches for large-scale systems normally consider unstructured uncertainty, in which (1) is a constant for all spatial frequency. In the proposed approach, however, we consider parametric uncertainty to improve user-friendliness, and therefore the key problem here is to calculate the singular values |(v.sub.j, 1)|, j, given the user-specified parametric uncertainties in (3), which can be mathematically formulated as

(96) $\begin{array}{cc}\underset{{}_p,{}_p,{}_p,{}_p,\mathrm{.Math.}}{\max }.Math.\left(v_j,1\right).Math.,\mathrm{for}j=1,\mathrm{.Math.},ns.t.{}_p\left[\underset{\_}{},\stackrel{\_}{}\right],{}_p\left[\underset{\_}{},\stackrel{\_}{}\right],{}_p\left[\underset{\_}{},\stackrel{\_}{}\right],{}_p\left[\underset{\_}{},\stackrel{\_}{}\right],\mathrm{.Math.}\left[\underset{\_}{\mathrm{.Math.}},\stackrel{\_}{\mathrm{.Math.}}\right].& \left(30\right)\end{array}$

(97) It is worth noting that analysis cannot be applied here because of the large dimensionality of CD processes (n can be as high as 300, and m as high as 2000) and the nonlinearity of f(x, .sub.p, .sub.p, .sub.p, .sub.p, c.sub.k.sup.p).

(98) Due to high nonlinearity of f(x, .sub.p, .sub.p, .sub.p, .sub.p, c.sub.k.sup.p), based on which G.sub.p(z) and G(z) are generated, explicit representations of the maximum singular values of at v.sub.j, j=1, . . . , n, are hard to achieve. Therefore the following approach is utilized to obtain an approximate solution to the problem in (30). Given the temporal frequency z=1, as the singular values of can be obtained by its spatial frequency gains, we can analyze how the spatial gain changes along the spatial frequency with respect to each of the spatial parameters. We increase each of the parameters from its lower bound (50% of the nominal value) to the upper bound (150% times the nominal value) in (3) with all the other parameters fixed, and analyze the change of the spatial gains. The results of this analysis are shown in FIGS. 13(a) to 13(e) for a typical process, which indicate that the spatial gain at each frequency is decreasing monotonically when a specified parameter is increasing from its lower bound to its nominal value (light lines) while increasing when the parameter increases within the other half of the region (darker lines).

(99) It is observed that in FIG. 13, |(v.sub.j, 1)| for j=1, . . . , n, is charging monotonically with respect to each of .sub.p, .sub.p, .sub.p, .sub.p, , so that the extreme values should occur at the boundaries of the uncertainty regions. Thus, the optimized problem in (30) is equivalent to:

(100) $\begin{array}{cc}\underset{{}_p,{}_p,{}_p,{}_p,\mathrm{.Math.}}{\max }.Math.\left(v_j,1\right).Math.,\mathrm{for}j=1,\mathrm{.Math.},ns.t.{}_p\left\{\underset{\_}{},\stackrel{\_}{}\right\},{}_p\left\{\underset{\_}{},\stackrel{\_}{}\right\},{}_p\left\{\underset{\_}{},\stackrel{\_}{}\right\},{}_p\left\{\underset{\_}{},\stackrel{\_}{}\right\},\mathrm{.Math.}\left\{\underset{\_}{\mathrm{.Math.}},\stackrel{\_}{\mathrm{.Math.}}\right\}.& \left(31\right)\end{array}$
Note that in the optimization problem in (31), only the extreme case spatial parameters are considered, which simplifies the problem at hand.

(101) Although the monotonicity property cannot be rigorously proved, it is intuitive in that the maximum singular value, which characterizes the system behavior, usually happens at the largest amount of model mismatch. This is also verified by extensive simulations, in which we compare the maximum singular values of from the proposed approach with that of a large number (>500) of 's generated randomly within the pre-specified uncertainty region for different processes. We only show the results for one of the typical processes, which is a single beam CD process that the dry weight is controlled by the primary autoslice, and the pre-specified parametric uncertainties are

(102) $\begin{array}{cc}{}_p\left[0.5,1.5\right],{}_p\left[0.8,1.2\right],{}_p\left[0.8,1.2\right],{}_p\left[0.5,1.5\right],\mathrm{.Math.}\left[-0.6\frac{m}{n},0.6\frac{m}{n}\right].& \left(32\right)\end{array}$

(103) The results are shown in FIG. 14, in which the curves are the singular values of randomly generated systems along the spatial frequency while the lower broken bound is the maximum singular values calculated by solving the optimization problem in (31). We can see that the broken line bound is a tight bound for all possible singular values due to model mismatch.

(104) Based on the tuning guideline of qf and the robust stability condition in (29), the desired qf can be found by utilizing a bisection search for the smallest qf that satisfies (29), which can be normally completed in about 2 second on a desktop with Intel i5 core and 6G memory.

(105) As illustrated in FIG. 15, the automated robust spatial tuning of the CD-MPC includes six steps:

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

(107) (2) Calculate the worst-case cut-off frequency v.sub.w for all the possible models given the parametric uncertainty. The worst-case cut-off frequency refers to the smallest cut-off frequency of all possible models given the parametric uncertainty specifications (the definition of cut-off frequency is shown in Definition 1). Based on this definition, the cut-off frequencies of a large number of models based on the parametric uncertainty need to be calculated to find the worst-case cut-off, (illustrated in equation (22)). Based on Lemma 3 and the monotonicity property (shown in FIG. 10), the worst-case cut-off can be calculated only based on several extreme models (illustrated in equation (23)), which significantly reduces computational time; this guarantees the design/tuning procedure of the digital controller can be accomplished within the required amount of time.

(108) (3) Design a weighting matrix S.sub.b to penalize high frequency actuator variability based on the process model and worst-case cut-off frequency. The S.sub.b is a weighting matrix in the fourth penalty term (Q4) of the MPC cost function as shown in equation (6), and it controls the actuator bending/picketing behavior. The high frequency actuator variability refers to the actuator bending/picketing behaviour that includes frequency components beyond the worst-case cut-off frequency. To prevent the high frequency actuator variability without unnecessarily limit the actuator movement, it is desired to design S.sub.b as a spatial filter whose spatial frequency response is (or is close to) the mirror of the actuator frequency response (shown in FIG. 8) so that the penalty level at each frequency is determined based on the spatial frequency property of the process. Lemma 2 helps to achieve this idea and also guarantees the generated S.sub.b can be applied in the digital controller.

(109) (4) Develop a robust spatial stability condition based on the parametric uncertainty specifications. Robust spatial condition refers to the requirement on the spatial frequency response of the closed-loop system, which guarantees the robust stability under the parametric uncertainty. Given the robust spatial condition, the frequency responses of a large amount of uncertain models generated based on the parametric uncertainty need to be tested (illustrated in equation (30)), which requires a lot of computational time; based on the monotonicity property between the model parameters and the corresponding frequency response (shown in FIG. 13) discovered in the invention, the condition can be checked based on only several critical frequency responses (illustrated in equation (31)); this saves the computational time and guarantees the design/tuning procedure of the digital controller can be accomplished within the required amount of time.

(110) (5) Adjust the multiplier of S.sub.b in the MPC cost function to satisfy the robust stability condition. The multiplier of S.sub.b refers to the scalar weight of the S.sub.b matrix defined in equation (6).

(111) (6) Output the weighting matrix S.sub.b and its multiplier.

(112) IV. Simulation Results with Real-Time Simulator

(113) We apply the proposed new S.sub.b and the automated robust 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.

(114) In this test, the parametric uncertainty specifications are given as

(115) $\begin{array}{cc}{}_p\left[0.75,1.25\right],{}_p\left[0.75,1.25\right],{}_p\left[0.75,1.25\right],{}_p\left[0.75,1.25\right],\mathrm{.Math.}\left[-0.6\frac{m}{n},0.6\frac{m}{n}\right].& \left(33\right)\end{array}$

(116) By utilizing the technique in (23), the worst cut-off frequency v.sub.w of the process model for the given parametric uncertainties in (33) can be obtained, and then the spectra penalization weighting matrix S.sub.b is designed based on the stop band v.sub.b=v.sub.w. The spatial tuning results in a qf that equals 0.026 (see FIG. 16, which shows that the stability condition in (29) is satisfied, and the temporal parameters are tuned by the existing method under 50% uncertainty on the temporal parameters, i.e., a.sub.p belongs to [0.5a,1.5a]). 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 (33). The performance is shown in FIGS. 17(a) and 17(b). FIG. 17(a) shows the steady-state measurement profile, while FIG. 17(b) is for the corresponding actuator profile. In the figures, it is evident that the profiles using the proposed method are much better than that of the existing method as the input/output profile is less oscillatory. More specifically, the 2 (two times of the standard deviation) of the steady-state actuator profile has been improved by about 66% while that of the measurement profile has been improved by about 44%.

(117) 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.