Patent Publication Number: US-10309059-B2

Title: Method of designing model predictive control for cross directional flat sheet manufacturing processes to guarantee temporal robust stability and performance

Description:
FIELD OF THE INVENTION 
     The present invention generally relates to monitoring and controlling continuous sheetmaking systems using model-based controllers and more specifically to techniques for designing a cross-directional (CD) model predictive controller (MPC) that provides robust temporal control performance in the presence of parametric model uncertainty. 
     BACKGROUND OF THE INVENTION 
     The control problems involved in paper-making processes can be divided into machine-directional (MD) control and cross-directional (CD) control. MD control concerns the paper properties along the direction towards which the paper moves while CD control focuses on the direction that is perpendicular to the machine direction. 
     CD control aims to reduce the variability of the paper property along the cross direction and to tune the dynamical property to meet the end users&#39; specifications. The paper property is measured by a scanner mounted downstream traversing back and forth across the paper sheet; various feedback control strategies are proposed to achieve consistency of the paper profile. CD control is a challenging control problem that may involves hundreds of process actuators and hundreds or thousands of process measurements, and process models typically have a large amount of uncertainty associated with them. There are spatial and temporal aspects to this problem. The spatial aspect relates to variability of the process measurements across the sheet while the temporal aspect relates to variability of each process measurement over time. 
     Model predictive control (MPC), a control strategy which takes control and state constraints explicitly into consideration, has seen thousands of applications in industry, and has been recently introduced into CD control in paper-making processes with the advance of computational capability as well as the development of fast QP solvers. The spatially-distributed CD process is a two-dimensional (spatial and temporal) frequency process, and the spatial response and temporal response are decoupled. Consequently, the controller tuning of CD processes can be separated into spatial tuning and temporal tuning. Spatial tuning aims to tune the weighting matrices such that the steady-state paper property across the paper sheet is consistent; temporal tuning concerns more about the satisfaction of performance indices in terms of settling times and control signal overshoots. 
     Several research results about temporal tuning have been reported. In particular, a constructive procedure to design spatially-distributed feedback controllers has been developed and applied for paper-making processes. Some stability margin and parameter tuning criteria were obtained via rectangular circulant matrices (RCMs) for the unconstrained CD-MPC, which provided a guide in the parameter tuning algorithms. Furthermore, an approximate steady-state performance prediction technique was proposed to speed up the parameter tuning procedure for the constrained CD-MPC. An automated tuning method was presented for the CD process such that the performance and robustness could be simultaneously satisfied under unstructured uncertainty. 
     As the spatial and temporal aspects of the CD control design are handled independently, robust control performance in the time domain is focused for the temporal aspect. The corresponding design objective is that the control remains temporally stable and the desired CD profile temporal performance is guaranteed in spite of uncertainty about the process model. However, there is lack of an easy-to-use technique to tune a CD controller to provide a robust temporal control performance that is characterized by intuitive indices, e.g., settling time, and to consider parametric model uncertainty, which is easy to understand by the user but is hard to handle for the design. 
     SUMMARY OF THE INVENTION 
     The present invention is based in part on the development of a CD-MPC design that 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 (from an initial high variability profile to a final low variability profile) and does not exceed an overshoot limit for actuator 2 sigma spread. The inventive technique includes: (i) providing parametric 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. 
     In one aspect, the invention is directed to a system which forms a material in a spatially-distributed multivariable-array cross-directional process wherein the system includes: 
     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 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 includes:
         (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.       

     In another aspect, the invention is directed to 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 includes 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.       

     The present invention is particularly suited for so-called “single-beam” controller tuning which is tuning the MPC with respect to a single actuator array comprising of a plurality of manipulated actuators that are arranged in the CD and a corresponding single controlled measurement array. A feature of the present invention is a tuning algorithm in which temporal tuning aims to tune the parameters such that the performance indices such as the settling time of the measurement profile and the control signal overshoot, which can be easily understood by the technician-level operators and maintenance personnel, are satisfied. This can usually be realized by adding a weighted control effort in the associated cost function. However, if a smaller weighting parameter is chosen, then a larger control signal overshoot will be observed. On the other hand, if a larger weighting parameter is chosen, then it will take a longer time for the output profile to achieve a steady-state value. With the present invention, a temporal filter is adopted to smooth the reference trajectory, and a systematic procedure is developed to tune the parameter in the temporal filter to realize a trade off between the settling time and the control signal overshoot. 
     Ubiquitously existing uncertainties make the temporal tuning more challenging. The parameters in CD processes are usually identified through bump tests and are inevitably subject to model uncertainties, which may be highly structured. Thus, it is necessary that temporal tuning algorithms be robust to the uncertainties. In existing tuning techniques, the mentioned uncertainty is normally characterized by the unstructured uncertainty in the frequency domain, which is hard to understand by technician-level operators and maintenance personnel. In the present invention, the uncertainty is considered by parametric uncertainty to improve the user-friendliness of the technique. 
     Novel aspects of the present invention include: 
     (1) worst-case temporal performance indices are defined and a technique to evaluate all the possible performance under the pre-specified parametric uncertainty is proposed. 
     (2) an automatic tuning algorithm for the parameter in the temporal filter is presented to guarantee satisfactory performance in terms of the worst-case settling time and control signal overshoot under user-specified parametric uncertainty. In the presence of parametric uncertainties, the tuning is based on the worst-case situation, which ensures that the required performance indices are satisfied. 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
         FIGS. 1, 2, and 3  are schematic illustrations of a papermaking system; 
         FIG. 4  is a block diagram of a closed-loop CD-MPC system; 
         FIG. 5  is a block diagram of a closed-loop CD-MPC system; 
         FIGS. 6( a ) and 6( b )  illustrate robust performance visualization for the 2σ spreads for the output and input, respectively; 
         FIG. 7  illustrates the minimum limit on the maximum peak on the sensitivity function; 
         FIG. 8  is a flow chart of the automated temporal tuning; 
         FIGS. 9( a ) and 9( b )  illustrate visualization technique for both worst- and best-case 2σ spreads; 
         FIGS. 10( a ) and 10( b )  are the temporal tuning result for the proposed approach; and 
         FIGS. 11( a ) to 11( b )  are performance illustrations for the proposed robust temporal tuning method. 
     
    
    
     DETAILED DESCRIPTION OF PREFERRED EMBODIMENTS 
     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. 
     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. 
       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 . 
     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. 
     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. Nos. 6,086,237 to Gorinevsky et al. and 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 
     As shown in  FIG. 4 , the closed-loop CD-MPC structure for a paper-making process consists of four parts: the real process G p (z) ( 74 ), the nominal model G(z) ( 76 ), the MPC ( 72 ) and the temporal filter F α (z) ( 70 ), which are illustrated in detail in the following. Moreover, the associated signals are defined as follows: y tgt , y sp , u(z), d(z), y 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 
     The nominal model G(z) of a CD paper-making process is characterized by 
                       G   ⁡     (   z   )       =       G   0     ⁢     h   ⁡     (   z   )           ,       h   ⁡     (   z   )       =         (     1   -   a     )     ⁢     z     -     t   d             1   -     az     -   1             ,           (   1   )               
where G 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 d  are the time constant and time delay in the discretized version.
 
     The spatial gain matrix G 0  has the parameterized structure as shown below: 
     
       
         
           
             
               
                 
                   
                     
                       G 
                       0 
                     
                     = 
                     
                       
                         [ 
                         
                           
                             g 
                             1 
                           
                           , 
                           … 
                           ⁢ 
                           
                               
                           
                           , 
                           
                             g 
                             n 
                           
                         
                         ] 
                       
                       ∈ 
                       
                         ℝ 
                         
                           m 
                           × 
                           n 
                         
                       
                     
                   
                   , 
                   
                     
 
                   
                   ⁢ 
                   
                     
                       
                         
                           
                             g 
                             k 
                           
                           = 
                             
                           ⁢ 
                           
                             f 
                             ⁡ 
                             
                               ( 
                               
                                 x 
                                 , 
                                 γ 
                                 , 
                                 η 
                                 , 
                                 ξ 
                                 , 
                                 β 
                                 , 
                                 
                                   c 
                                   k 
                                 
                               
                               ) 
                             
                           
                         
                       
                     
                     
                       
                         
                           
                             = 
                               
                             ⁢ 
                             
                               
                                 γ 
                                 2 
                               
                               ⁢ 
                               
                                 { 
                                 
                                   
                                     
                                       
                                         
                                           
                                             e 
                                             
                                               - 
                                               
                                                 
                                                   
                                                     η 
                                                     ⁡ 
                                                     
                                                       ( 
                                                       
                                                         x 
                                                         - 
                                                         
                                                           c 
                                                           k 
                                                         
                                                         + 
                                                         
                                                           β 
                                                           ⁢ 
                                                           
                                                               
                                                           
                                                           ⁢ 
                                                           ξ 
                                                         
                                                       
                                                       ) 
                                                     
                                                   
                                                   2 
                                                 
                                                 
                                                   ξ 
                                                   2 
                                                 
                                               
                                             
                                           
                                           ⁢ 
                                           
                                             cos 
                                             ⁡ 
                                             
                                               [ 
                                               
                                                 
                                                   π 
                                                   ⁡ 
                                                   
                                                     ( 
                                                     
                                                       x 
                                                       - 
                                                       
                                                         c 
                                                         k 
                                                       
                                                       + 
                                                       
                                                         β 
                                                         ⁢ 
                                                         
                                                             
                                                         
                                                         ⁢ 
                                                         ξ 
                                                       
                                                     
                                                     ) 
                                                   
                                                 
                                                 ξ 
                                               
                                               ] 
                                             
                                           
                                         
                                         + 
                                       
                                     
                                   
                                   
                                     
                                       
                                         
                                           e 
                                           
                                             - 
                                             
                                               
                                                 
                                                   η 
                                                   ⁡ 
                                                   
                                                     ( 
                                                     
                                                       x 
                                                       - 
                                                       
                                                         c 
                                                         k 
                                                       
                                                       - 
                                                       
                                                         β 
                                                         ⁢ 
                                                         
                                                             
                                                         
                                                         ⁢ 
                                                         ξ 
                                                       
                                                     
                                                     ) 
                                                   
                                                 
                                                 2 
                                               
                                               
                                                 ξ 
                                                 2 
                                               
                                             
                                           
                                         
                                         ⁢ 
                                         
                                           cos 
                                           ⁡ 
                                           
                                             [ 
                                             
                                               
                                                 π 
                                                 ⁡ 
                                                 
                                                   ( 
                                                   
                                                     x 
                                                     - 
                                                     
                                                       c 
                                                       k 
                                                     
                                                     - 
                                                     
                                                       β 
                                                       ⁢ 
                                                       
                                                           
                                                       
                                                       ⁢ 
                                                       ξ 
                                                     
                                                   
                                                   ) 
                                                 
                                               
                                               ξ 
                                             
                                             ] 
                                           
                                         
                                       
                                     
                                   
                                 
                                 } 
                               
                             
                           
                           , 
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   2 
                   ) 
                 
               
             
             
               
                 
                   
                     x 
                     = 
                     1 
                   
                   , 
                   … 
                   ⁢ 
                   
                       
                   
                   , 
                   m 
                   , 
                   
                     k 
                     = 
                     1 
                   
                   , 
                   … 
                   ⁢ 
                   
                       
                   
                   , 
                   n 
                   , 
                 
               
               
                 
                     
                 
               
             
           
         
       
     
     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 th  actuator, c k  is the alignment parameter that determines the center of the corresponding spatial response. 
     Model Uncertainty 
     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 p (z) is then described by 
                         G   p     ⁡     (   z   )       =       G   0   p     ⁢       h   p     ⁡     (   z   )           ,     
     ⁢       G   0   p     =     [       g     1   ⁢   p       ,   …   ⁢           ,     g   np       ]       ,     
     ⁢       g   kp     =     f   ⁡     (     x   ,     γ   p     ,     η   p     ,     ξ   p     ,     β   p     ,     c   k   p       )         ,     
     ⁢       γ   p     =       (     1   +     r   γ       )     ⁢   γ       ,       η   p     =       (     1   +     r   η       )     ⁢   η       ,     
     ⁢       ξ   p     =       (     1   +     r   ξ       )     ⁢   ξ       ,       β   p     =       (     1   +     r   β       )     ⁢   β       ,     
     ⁢       c   k   p     =       c   k     +   ɛ       ,     
     ⁢         h   p     ⁡     (   z   )       =         (     1   -     a   p       )     ⁢     z     -     t   dp             1   -       a   p     ⁢     z     -   1               ,     
     ⁢       a   p     =       (     1   +     r   a       )     ⁢   a       ,       t   dp     =       (     1   +     r     t   d         )     ⁢     t   d         ,     
     ⁢     k   =   1     ,   …   ⁢           ,   n   ,     x   =   1     ,   …   ⁢           ,   m   ,           (   3   )               
where r γ ∈[ r   γ ,  r   γ ], r η ∈[ r   η ,  r   η ], r ξ ∈[ r   ξ ,  r   ξ ], r β ∈[ r   β ,  r   β ], r α ∈[ r   α ,  r   α ], r td ∈[ r   td ,  r   td ], and ε∈[− ε (m/n),  ε (m/n),] are utilized to characterize the parametric uncertainties. These trust ranges are also represented as γ p ∈[ γ ,  γ ], η p ∈[ η ,  η ], ξ p ∈[ ξ ,  ξ ], β p ∈[ β ,  β ], α p ∈[ α ,  α ], and t dp ∈[ t   d ,  t   d ] for brevity; note that as ε is a global perturbation on all c k   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 γ p , η p , β p , ξ p , ε, α p , t dp .
 
CD Model Predictive Controller
 
     For industrial CD-MPC controllers that are applied in paper mills, the following optimization problem is solved: 
                       min     Δ   ⁢           ⁢     u   ⁡     (   k   )           ⁢     {         ∑     i   =   1       H   p       ⁢         (       y   ⁡     (     k   +   i     )       -       y   sp     ⁡     (     k   +   i     )         )     T     ⁢       Q   1     ⁡     (       y   ⁡     (     k   +   i     )       -       y   sp     ⁡     (     k   +   i     )         )           +       ∑     i   =   0         H   u     -   1       ⁢     [       Δ   ⁢           ⁢       u   ⁡     (     k   +   i     )       T     ⁢     Q   2     ⁢   Δ   ⁢           ⁢     u   ⁡     (     k   +   i     )         +         (       u   ⁡     (     k   +   i     )       -       u   sp     ⁡     (     k   +   i     )         )     T     ⁢       Q   3     ⁡     (       u   ⁡     (     k   +   i     )       -       u   sp     ⁡     (     k   +   i     )         )         +         u   ⁡     (     k   +   i     )       T     ⁢     Q   4     ⁢     u   ⁡     (     k   +   i     )           ]         }       ,           (   4   )               
subject to the system dynamics defined in (1) and the constraints as follows:
 
ΩΔ u ( k )≤ b−Γu ( k− 1),  (5)
 
where H p  is the prediction horizon, and H u  is the control horizon; y(k) ∈R m  and y sp (k)∈R m  are the predicted output profile and the corresponding reference signal; u(k) ∈R n  and u sp (k) ∈R n  are the actuator profile and its reference; Δu(k)(=u(k)−u(k−1)) is the change in the actuator profile; Q 1  to Q 3  are diagonal weighting matrices; Q 4  is the weighting matrix on the actuator bending/picketing in the following form:
 
                       Q   4     =       q   4     ⁢     S   b   T     ⁢     S   b         ,     
     ⁢       S   b     =     [           -   1         1       0       …       …       …       0           1         -   2         1       ⋱       ⋱       ⋱       ⋮           0       1         -   2         ⋱       ⋱       ⋱       ⋮           ⋮       ⋱       ⋱       ⋱       ⋱       ⋱       ⋮           ⋮       ⋱       ⋱       ⋱         -   2         1       0           ⋮       ⋱       ⋱       ⋱       1         -   2         1           0       …       …       …       0       1         -   1           ]       ,           (   6   )               
where q 4  is a scalar weight and S b  ∈R n×n  is the “bending moment matrix”. Note that for the actuator profile, the first and second order derivatives are incorporated in the matrix S 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
 
     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 sp (k) based on
 
 Y   sp ( k )= F   α ( y   tgt ( k )− d   y ( k )),  (7)
 
where y tgt  (k) is the output target, and d y (k)=y p (k)−y(k) is the disturbance estimated based on the process output y p (k) and predicted output y(k). F α  is the time domain implementation of f α  (z) based on y sp (z)=f α  (z)I m (y tgt  (z)−d y (z)) and f α  (z) is the temporal filter
 
                         f   α     ⁡     (   z   )       =         (     1   -     a   r       )     ⁢     z     -     t   d             1   -       a   r     ⁢     z     -   1               ,           (   8   )               
where α r =e −ΔT/ασ ; ΔT is the sampling time, and σ is the continuous-time time constant of the temporal transfer function of the process; I 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 2  set to a small-valued scalar matrix.
 
Closed-Loop Transfer Functions
 
     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 r  (z) ( 80 ) and K α  (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. 
     In  FIG. 5 , G 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   p ( z )= G ( z )+Δ( z ),  (9)
 
in which Δ(z) ( 84 ) denotes the model uncertainty.
 
     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 p (z) if it is nominally stable and
 
∥ T   ud ( z )Δ( z )∥ ∞ &lt;1→ σ ( T   ud ( e   iω )Δ( e   iω ))&lt;1,∀ω,  (10)
 
 T   ud ( z )= K   α ( z )[ I−G ( z ) K   α ( z )] −1 ,  (11)
 
where  σ  (A) denotes the maximum singular value of a matrix A; Δ(z) is the model uncertainty in (9) and T ud  (z) ∈C n×m  is the sensitivity function from the disturbance profile d(z) to the input profile u(z) at the nominal case.
 
     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   yd ( z )=[ I−G ( z ) K   α ( z )] −1 ,  (12)
 
where T yd (z) ∈C m×m .
 
     T ud  (z) and T yd  (z) are employed as the key transfer functions in this paper for evaluating the closed-loop system performance. 
     Given the CD-MPC structure shown in  FIG. 4 , the closed-loop system behavior is affected by the following parameters: the penalty matrices Q 1  to Q 4 , the temporal filter parameter α, the prediction horizon H p , and the control horizon H u . 
     Note that in industrial practices, the prediction horizon H p  is normally selected to be 4 times the summation of the time constant and delay; the control horizon H u =1 is usually utilized in large-scale MPC (can also be increased based on practical situations); Q 1  is often fixed in robust tuning. Consequently, only Q 3 , Q 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&#39; specifications. 
     The traditional idea is to conduct the tuning based on the sensitivity functions T yd  ∈ C m×m  and T ud  ∈C n×m  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 
     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   m   G ( z ) F   n   H ,  (13)
 
where F m  and F n  are the complex Fourier matrices with dimension m and n, respectively. Without loss of generality, we assume m&gt;n in this analysis. Then, all the frequency information of G(z) is contained in the following transfer functions:
 
{{hacek over ( g )}( v   1   ,z ), . . . ,{hacek over ( g )}( v   n   ,z )}= DIAG ( F   m   G ( z ) F   n   H ),  (14)
 
where DIAG(A) denotes the operation of getting the following elements of a rectangular matrix A ∈C m×n : DIAG(A) 1 ={A(1, 1), . . . , A(k, k), A(k+1+m−n, k+1), . . . , A(m, n)}, where k=n/2 if n is even, or k=(n+1)/2 if n is odd; and
 
                         g   ˘     ⁡     (       v   j     ,   z     )       =           g   ˘     ⁡     (     v   j     )       ⁢     (     1   -   a     )     ⁢     z     -     t   d             1   -     az     -   1             ,       for   ⁢           ⁢   j     =   1     ,   …   ⁢           ,     n   .             (   15   )               
Note that based on (15), the spatial frequency response gain of G(z) at v j  is |{hacek over (g)}(v j )|, where v 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 yd  (z), T ud  (z) are also RCMs. Thus, they can also be analyzed in the two-dimensional frequency domain as follows:
 
 Ť   yd ( z )= F   m   T   yd ( z ) F   m   H ,
 
 Ť   ud ( z )= F   n   T   ud ( z ) F   m   H ,
 
{ ť   yd ( v   1   ,z ), . . . , ť   yd ( v   m   ,z )}=diag( Ť   yd ( z )),
 
{ ť   ud ( v   1   ,z ) . . . , ť   ud ( v   n   ,z )}= DIAG ( Ť   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 ť ud  (v, z) and ť 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 iω , ω=0) via Q 3  and Q 4 , and temporal tuning at v=0 via the parameter α (which has the same function as Q 2 ). In the following, the temporal tuning is explored.
 
III. Temporal Tuning for CD-MPC
 
     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. 
     Based on the small gain theory and the two dimensional analysis, the robust stability condition for temporal tuning can be represented as 
                                  t   ˘     ud     ⁡     (     0   ,     e     i   ⁢           ⁢   ω         )       ⁢     δ   ⁡     (     0   ,     e     i   ⁢           ⁢   ω         )              &lt;   1     ,     
     ⁢       ∀   ω     -&gt;                t   ˘     ud     ⁡     (     0   ,     e     i   ⁢           ⁢   ω         )            &lt;     1          δ   ⁡     (     0   ,     e     i   ⁢           ⁢   ω         )                  ,     ∀   ω     ,           (   17   )               
where ť ud  (0, e iω ) and δ (0, e iω ) are T ud  (z) and Δ(z) at spatial frequency v=0. Note that as G p (z) is an RCM because of the parametric uncertainty, Δ(z)=G 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 iω ) denotes its two dimensional form).
 
     As the extreme system behavior normally happens at the extreme model parameters, the maximum value of |δ(0, e iω )| can be obtained by solving 
     
       
         
           
             
               
                 
                   
                     
                       
                         max 
                         
                           
                             a 
                             p 
                           
                           , 
                           
                             t 
                             dp 
                           
                         
                       
                       ⁢ 
                       
                          
                         
                           δ 
                           ⁡ 
                           
                             ( 
                             
                               0 
                               , 
                               
                                 e 
                                 
                                   i 
                                   ⁢ 
                                   
                                       
                                   
                                   ⁢ 
                                   ω 
                                 
                               
                             
                             ) 
                           
                         
                          
                       
                     
                     , 
                     
                       ∀ 
                       ω 
                     
                   
                   ⁢ 
                   
                     
 
                   
                   ⁢ 
                   
                     
                       
                         s 
                         . 
                         t 
                         . 
                         
                             
                         
                         ⁢ 
                         
                           a 
                           p 
                         
                       
                       ∈ 
                       
                         { 
                         
                           
                             a 
                             _ 
                           
                           , 
                           
                             a 
                             _ 
                           
                         
                         } 
                       
                     
                     , 
                     
                       
                         t 
                         dp 
                       
                       ∈ 
                       
                         
                           { 
                           
                             
                               
                                 t 
                                 _ 
                               
                               d 
                             
                             , 
                             
                               
                                 t 
                                 _ 
                               
                               d 
                             
                           
                           } 
                         
                         . 
                       
                     
                   
                 
               
               
                 
                   ( 
                   18 
                   ) 
                 
               
             
           
         
       
     
     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 |ť ud (0, e 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. 
     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. 
     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. 
     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. 
     As uncertainties are considered on the dynamic parameters α and t 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: 
                         min       a   p     ,     t   dp         ⁢     2   ⁢     σ   ⁡     (       y   i     ⁡     (       a   p     ,     t   dp       )       )           ,       for   ⁢           ⁢   i     =   1     ,   …   ⁢           ,   l     ⁢     
     ⁢         s   .   t   .           ⁢     a   p       ∈     [       a   _     ,     a   _       ]       ,       t   dp     ∈     [         t   _     d     ,       t   _     d       ]       ,             (   19   )                     min       a   p     ,     t   dp         ⁢           ⁢     2   ⁢     σ   ⁡     (       y   i     ⁡     (       a   p     ,     t   dp       )       )           ,       for   ⁢           ⁢   i     =   1     ,   …   ⁢           ,   l     ⁢     
     ⁢         s   .   t   .           ⁢     a   p       ∈     [       a   _     ,     a   _       ]       ,       t   dp     ∈     [         t   _     d     ,       t   _     d       ]       ,             (   20   )               
where l is the entire scan time for the process. Note that as the spatial gain matrix G 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 2σs).
 
     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 α p  ∈{ α ,  α } and t dp  ∈{ t   d ,t   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:
 
α p ∈[0.5 α, 1.5α],t dp ∈[0.5 td ,1.5 t   d ].  (21)
 
     The tuning objective of the temporal part is then formulated as 
                       min   α     ⁢       T   s     ⁡     (     2   ⁢     σ   ⁡     (   y   )         )         ⁢     
     ⁢         s   .   t   .           ⁢     OS   ⁡     (     2   ⁢     σ   ⁡     (   u   )         )         ≤     OS   *       ,             (   22   )               
where T 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 yd  (z) at v=0 is first selected, and then we find an a to meet the bound while still guaranteeing robust stability.
 
     Define κ l  as the minimum limit of the maximum peak of the sensitivity function (see  FIG. 7 ), the objective here is to adjust α so that
 
| ť   yd (0, e   iω )|&lt;κ l ,∀ω,  (23)
 
is satisfied. Then, the frequency domain tuning can be realized by solving:
 
                     α   f     =       max       k   =   0     ,   …   ⁢           ,     k     ma   ⁢           ⁢   x           ⁢     (     α     ω   k     m     )               (   24   )                   α     ω   k     m     =     min   ⁡     (     α     ω   k       )         ⁢     
     ⁢         s   .   t   .           ⁢     α     ω   k         &gt;     α   *       ,                t   ˘       y   ⁢           ⁢   d       ⁡     (     0   ,     e     i   ⁢           ⁢   ω         )            &lt;     κ   l       ,             (   25   )               
where ω k , k=0, . . . , k max , are the temporal frequencies distributed from 0 to the Nyquist frequency. κ 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&gt;OS*, and thus the time domain tuning is implemented to fine tune the parameter.
 
     In the time domain tuning, the parameter α is fine tuned to find the optimal α (denoted as α 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 α, the bisection search can be utilized to find the α t  based on the closed-loop simulation with the following procedure. If the OS with α f  is larger than the specification, the search is implemented in the region α∈[α f , α u ], where α u  is the upper bound for α with κ l =1.2, and α f  is obtained in the frequency tuning. Otherwise, the search is implemented in the region α ∈ [α i , α f ], where α l  is the lower bound for a with κ l =1.8. To sum up, the temporal tuning part is implemented using the following algorithm: 
     Algorithm 1 Temporal Tuning 
     
         
         
           
             1: Input OS*; 
             2: Input the parametric uncertainty specifications [ α ,  α ] and [ t   d , t   d   ]; 
             3: Calculate the lower bound for robust stability α*; 
             4: Implement the frequency tuning to achieve α f ; 
             5: if OS(α f )&gt;OS* then 
             6: Search in [α f , α u ] for α t ; 
             7: else 
             8: Search in [α l , α f ] for α t ; 
             9: end if 
             10: End 
           
         
       
    
     Implementing the tuning procedure in Algorithm 1, the automated temporal tuning costs about 6 seconds on a computer with i5 CORE and 6G RAM. 
     Remark 1. To save more computation time, we can also stop the algorithm if OS(α f )≤OS*. This can avoid the tuning in the time domain once the requirement is satisfied in the frequency tuning. 
     Remark 2. 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. 
     As illustrated in  FIG. 8 , automatic temporal tuning of the CD-MPC includes of the following 7 steps: 
     (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. 
     (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. 
     (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. 
     (4) Calculate a lower bound for the profile trajectory tuning parameter α 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. 
     (5) Use a frequency domain technique to further reduce the search range for α. The “search range for α” refers to the range of the tuning parameter within which the optimal tuning parameter may exist. 
     (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 2σ 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. 
     (7) Output the tuning parameter α. 
     IV. Simulation Results with Real-Time Simulator 
     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. 
     In this test, the parametric uncertainty specifications are given as 
     
       
         
           
             
               
                 
                   
                     
                       γ 
                       p 
                     
                     ∈ 
                     
                       [ 
                       
                         
                           0.75 
                           ⁢ 
                           γ 
                         
                         , 
                         
                           1.25 
                           ⁢ 
                           γ 
                         
                       
                       ] 
                     
                   
                   , 
                   
                     
                       η 
                       p 
                     
                     ∈ 
                     
                       [ 
                       
                         
                           0.75 
                           ⁢ 
                           η 
                         
                         , 
                         
                           1.25 
                           ⁢ 
                           η 
                         
                       
                       ] 
                     
                   
                   , 
                   
                     
 
                   
                   ⁢ 
                   
                     
                       β 
                       p 
                     
                     ∈ 
                     
                       [ 
                       
                         
                           0.75 
                           ⁢ 
                           β 
                         
                         , 
                         
                           1.25 
                           ⁢ 
                           β 
                         
                       
                       ] 
                     
                   
                   , 
                   
                     
                       ξ 
                       p 
                     
                     ∈ 
                     
                       [ 
                       
                         
                           0.75 
                           ⁢ 
                           ξ 
                         
                         , 
                         
                           1.125 
                           ⁢ 
                           ξ 
                         
                       
                       ] 
                     
                   
                   , 
                   
                     
 
                   
                   ⁢ 
                   
                     
                       a 
                       p 
                     
                     ∈ 
                     
                       [ 
                       
                         
                           0.5 
                           ⁢ 
                           a 
                         
                         , 
                         
                           1.5 
                           ⁢ 
                           a 
                         
                       
                       ] 
                     
                   
                   , 
                   
                     
                       t 
                       
                         d 
                         p 
                       
                     
                     ∈ 
                     
                       [ 
                       
                         
                           0.5 
                           ⁢ 
                           
                             t 
                             d 
                           
                         
                         , 
                         
                           1.5 
                           ⁢ 
                           
                             t 
                             d 
                           
                         
                       
                       ] 
                     
                   
                   , 
                   
                     
 
                   
                   ⁢ 
                   
                     ɛ 
                     ∈ 
                     
                       
                         [ 
                         
                           
                             
                               - 
                               0.6 
                             
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             
                               m 
                               n 
                             
                           
                           , 
                           
                             0.6 
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             
                               m 
                               n 
                             
                           
                         
                         ] 
                       
                       . 
                     
                   
                 
               
               
                 
                   ( 
                   26 
                   ) 
                 
               
             
           
         
       
     
     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%. 
     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))&lt;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 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. 
     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.