Optimal closed-loop input design for identification of flat-sheet process models

Abstract

Sheetmaking cross-directional (CD) control requires a sophisticated model-based controller whose operation requires an accurate model of process behavior, but due to the complexity of the process, identifying these process models is challenging. Current techniques rely on open-loop process experimentation. Using non-causal scalar transfer functions for the steady-state CD process model and controller model avoid the problem of large dimensions associated with the CD process. These non-causal transfer functions can be represented by causal transfer functions that are equivalent to the non-causal ones in the sense of the output spectrum. A closed-loop optimal input design framework is proposed based on these causal equivalent models. CD actuators have responses in both sides along the cross direction which can be viewed as a non-causal behavior. Techniques to perform the non-causal modeling are demonstrated and developed in a closed-loop optimal input design framework based on non-causal modeling of the closed-loop CD process.

Claims

1. A method of closed-loop identification of process models for the model predictive control (MPC) by a profile analyzer, of an industrial sheetmaking system having a plurality of actuators arranged in the cross-direction (CD) and a scanning sensor supported on a supporting frame wherein the profile analyzer manipulates the actuators to control the sheet parameters in a spatially distributed sheet process, the method comprises the steps of: (a) receiving from the scanning sensor, data indicative of the magnitude of a measured sheet property of the spatially-distributed sheet process; (b) selecting by the profile analyzer a process model for the spatially-distributed sheet process wherein the process model is defined by a matrix; (c) converting the matrix into a non-causal transfer function; (d) converting the non-causal transfer function into a causal model that has an equivalent spectrum; (e) using the causal model to design an optimal input spectrum for process excitation; (f) transforming a frequency domain representation of the optimal input spectrum to a time domain realization which is a sequence of actuator movements; (g) providing signals by the profile analyzer to the plurality of actuators, of the sequence of actuator movements and collecting data via the scanning sensor changes in cross-directional sheet properties due to the sequence of actuator movements to determine measured actuator response profiles; (h) analyzing by the profile analyzer, the data to extract new model parameters; and (i) using the new model parameters for selecting the process model.

2. The method of claim 1 wherein the matrix defines steady-state gains between actuator positions and spatially distributed process measurements.

3. The method of claim 1 wherein step (a) comprises of using an initial process model that is an existing process model or that is a developed process model to approximate process behavior wherein the initial process model allows for the design of an excitation sequence that is tailored to a specific process of interest.

4. The method of claim 1 wherein step (b) a non-causal spatial finite impulse response model is generated by taking parameters from a single column of a spatial gain matrix and the non-causal spatial impulse response model is factored into a causal transfer function and an identical but anti-causal transfer function.

5. The method of claim 1 wherein step (c) comprises taking the square of the causal transfer functions.

6. The method of claim 1 wherein in step (d) an input spectrum is designed to minimize the covariance of parameter estimates of the causal-equivalent model subject to constrains on input and output power.

7. The method of claim 6 wherein in step (d) to minimize the covariance of parameter estimates of the causal-equivalent model comprises making a finite dimensional parameterization of the input spectrum.

8. The method of claim 1 wherein the process model for the spatially-distributed process is with respect to one actuator array and a corresponding measurement array.

9. The method of claim 1 wherein the spatially-distributed sheetmaking process is a paper-making process.

10. A profile analyzer for providing multivariable model predictive control (MPC) to a cross-direction (CD) process having at least one manipulated actuator array and a scanning sensor supported on a supporting frame, and at least one controlled measurement array, wherein the profile analyzer manipulates the actuators to control the sheet parameters in a spatially distributed sheet process the profile analyzer comprises a processor that is configured to: receive from the scanning sensor, data indicative of the magnitude of a measured sheet property of the spatially-distributed sheet process; select a process model for the spatially-distributed sheet process wherein the process model is defined by a matrix; convert the matrix into a non-causal transfer function; convert the non-causal transfer function into a causal model that has an equivalent spectrum; use the causal model to design an optimal input spectrum for process excitation; transform a frequency domain representation of the optimal input spectrum to a time domain realization which is a sequence of actuator movements; provide, signals to the plurality of actuators, of the sequence of actuator movements and collecting data via the scanning sensor changes in cross-directional sheet properties due to the sequence of actuator movements to determine measured actuator response profiles; analyze the data to extract new model parameters; and use the new model parameters for the process model.

11. The system of claim 10 wherein the matrix defines steady-state gains between actuator positions and spatially distributed process measurements.

12. The system of claim 10 wherein the MPC employs an initial process model that is an existing process model or that is a developed process model to approximate process behavior wherein the initial process model allows for the design of an excitation sequence that is tailored to a specific process of interest.

13. The system of claim 10 wherein the processor is configured to convert the matrix into a non-transfer function whereby a non-causal spatial finite impulse response model is generated by taking parameters from a single column of a spatial gain matrix and the non-causal spatial impulse response model is factored into a causal transfer function and an identical but anti-causal transfer function.

14. The system of claim 10 wherein the processor is configured to convert the non-causal transfer function into a causal model that has an equivalent spectrum by taking the square of the causal transfer functions.

15. The system of claim 10 wherein the processor is configured to use the causal model to design an optimal input spectrum for process excitation whereby an input spectrum is designed to minimize the covariance of parameter estimates of the causal-equivalent model subject to constrains on input and output power.

16. The system of claim 15 wherein to minimize the covariance of parameter estimates of the causal-equivalent model comprises making a finite dimensional parameterization of the input spectrum.

17. The system of claim 10 wherein the process model is with respect to one manipulated actuator array and a corresponding measurement array.

18. A non-transitory computer readable medium embodying a computer program for identifying alternative process models for a profile analyzer having model predictive control (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, and a scanning sensor supported on a supporting frame, wherein the profile analyzer manipulates the actuators to control the sheet parameters in a spatially distributed sheet process, wherein the program comprises readable program code for: receiving from the scanning sensor, data indicative of the magnitude of a measured sheet property of the spatially-distributed sheet process; selecting by the profile analyzer a process model for the spatially-distributed sheet process wherein the process model is defined by a matrix; converting the matrix into a non-causal transfer function; converting the non-causal transfer function into a causal model that has an equivalent spectrum; using the causal model to design an optimal input spectrum for process excitation; transforming a frequency domain representation of the optimal input spectrum to a time domain realization which is a sequence of actuator movements; providing signals by the profile analyzer to the plurality of actuators, of the sequence of actuator movements and collecting data via the scanning sensor changes in cross-directional sheet properties due to the sequence of actuator movements to determine measured actuator response profiles; analyzing by the profile analyzer, the data to extract new model parameters; and using the new model parameters for the process model.

19. The non-transitory computer readable medium of claim 18 wherein the process models are with respect to one actuator array and a corresponding measurement array.

20. The non-transitory computer readable medium of claim 18 wherein the cross-directional process is a paper-making process.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

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

(2) FIG. 4 illustrates a closed-loop optimal input design configuration;

(3) FIGS. 5A and 5B depict the impulse response of a single actuator (solid line) and the impulse response of the estimated non-causal transfer function (dash dotted line); and

(4) FIGS. 6A and 6B depict the impulse responses of the estimated process model in the closed-loop under the optimally designed input (FIG. 6A) and the bumped input (FIG. 6B) in 100 Monte-Carlo simulations.

DETAILED DESCRIPTION OF PREFERRED EMBODIMENTS

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

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

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

(8) 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. In contrast, the inventive reverse bump test does not directly send perturbations to the actuator profile. 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.

(9) 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. In so-called bump tests operating parameters on the sheetmaking system, such as a papermaking machine, are altered and changes of certain dependent variables resulting therefrom are measured. Bump tests techniques are described, for open-loop CD process model estimation, in U.S. Pat. No. 6,086,237 to Gorinevsky et al. and, for closed-loop CD process alignment identification, in U.S. Pat. No. 7,459,060 to Stewart, which are incorporated herein by reference. While the invention will be described with respect to a papermaking machine, it is understood that the invention is applicable to other spatially-distributed processes such as plastic sheetmaking, rubber sheetmaking, and sheet metal operation.

(10) The invention will be illustrated with a closed-loop input design for the CD process. Since most of time the closed-loop CD process is operating at steady-state, the preferred process will be on the steady-state closed-loop CD process. The main challenges are how to deal with the large input-output dimensions of the process and how to incorporate the controller in order to perform the closed-loop optimal input design. With the present invention, a non-causal model for the CD process is developed to avoid the high dimensional issue associated with the conventional Multi-input-multi-output (MIMO) CD model. To eliminate the resultant difficulty for the input design, we then propose an approach to obtain a causal model with an equivalent output spectrum as the non-causal model. It is shown that the maximum likelihood estimate and the parameter covariance matrix of the causal-equivalent model will converge to those of the non-causal model asymptotically with probability one. In this sense, the optimal excitation signal can be designed directly based on the causal model.

(11) I. CD Processes and the Steady-State Model

(12) A. Open-loop CD Process Model

(13) In the following CD process model that includes controller 70 and plant 72 shown in FIG. 4,
y(t)=g(z.sup.1)Gu(t)+v(t),(1)
y(t) .sup.m represents the measured controlled variable (CV) profile, and m is the number of the measurement boxes along the cross direction. u(t) .sup.m is the manipulated variable (MV) profile and here we assume a square CD model. v(t) .sup.m is the disturbance acting on the output of the process. z.sup.1 is the unit backward-shift operator. g(z.sup.1) is used to describe the dynamics associated with each actuator and is assumed to be a scalar transfer function. In other words, we assume that all the actuators share the same dynamics, which is a common practice to simplify the modeling of the process. G .sup.mm is the steady-state gain matrix and each column of G matrix is indeed the sampled impulse response of a single actuator along the CD at steady-state. The most distinguished features of G are its ill-conditionality, Toeplitz structure and sparsity. These characteristics greatly reduce the complexity of the controller design and model identification of the CD process. For clarity, we pose the following assumption on the structure of G matrix.

(14) Assumption 1. The actuators of the CD process have the same and symmetric impulse response shape along the spatial direction at steady-state, except the center of the response shape of each actuator is different. The columns of G are indeed sampled version of these responses.

(15) For the dynamic model, we assume g(z.sup.1) to have the following form (the subscript t represents temporal),

(16) $\begin{array}{cc}g\left(z^{-1}\right)=z^{-d}\frac{B_t\left(z^{-1}\right)}{A_t\left(z^{-1}\right)}& \left(2\right)\end{array}$
where d is the time-delay. B.sub.t(z.sup.1) and A.sub.t(z.sup.1) are polynomials and mostly g(z.sup.1) is the discretization of a first-order plus time-delay model with unit gain. Similarly, the output disturbance v(t) is assumed to be filtered white noise, both temporally and spatially,

(17) $\begin{array}{cc}v\left(t\right)=\frac{D_t\left(z^{-1}\right)}{C_t\left(z^{-1}\right)}e\left(t\right)& \left(3\right)\end{array}$
where C.sub.t(z.sup.1) and D.sub.t(z.sup.1) are monic and scalar polynomials describing the filter in the temporal direction. Again, here we assume the disturbances affecting all the output channels have the same dynamic model. The constant matrix .sup.mm denotes the spatial correlation of the noise. e(t) is the white noise with zero mean and covariance [e(t)e.sup.T(t.sub.0)]=.sub.e(tt.sub.0), where is the expectation operator, .sub.e is the co-variance matrix and is the Dirac delta function.

(18) For the purpose of spatial optimal excitation signal design, the following steady-state CD process model is of interest,
y.sub.ss=G.sub.ssu.sub.ss+v.sub.ss,(4)
where y.sub.ss .sup.m is the steady-state measured CV profile, and u.sub.ss .sup.m is the steady-state MV profile. G.sub.ss=G is the steady-state process gain. v.sub.ss .sup.m is the steady-state output disturbance. For convenience, we suppose that the spatial filter is also Toeplitz-structured and sparse as G.sub.ss.
B. Closed-loop Steady-state CD Process Model

(19) It is well-known that the closed-loop dither signal design will involve the explicit expression of the controller in the formulation of the objective function. Most modern industrial CD controllers are MPC and if any of the constraints is active, the controller will become highly complicated and even nonlinear.

(20) The presence of a nonlinear controller in the control loop will add to the complexity of closed-loop optimal input design. Thus to simplify this procedure, we introduce the following assumption.

(21) Assumption 2. Throughout this analysis, the MPC is assumed to operate in the linear mode and no constraints are active.

(22) From Assumption 2, the specific expression of the model predictive controller at steady-state for the CD process is known to be,
K.sub.ss=Q.sub.3.sup.1.sub.KG.sub.ssQ.sub.1,(5)
where Q.sub.1 is the weight matrix in the MPC objective function penalizing the deviation of the CV profile from its set-points. Q.sub.3 is the corresponding weight matrix to penalize the offset of the steady-state manipulated variable (MV) from its target. .sub.K is a constant determined from the dynamic model (2) of the actuators. In practice, for convenience, the weighting matrices Q.sub.1 and Q.sub.3 are in general selected to be diagonal. From (5), it is observed that the controller K.sub.ss will possess a similar structure as the gain G.sub.ss. Awareness of this important point will greatly facilitate the proceedings of the derivations in the sequel.

(23) Combining (4) and (5), from FIG. 4, we can easily arrive at the closed-loop models of the CD process,
Y.sub.ss=(I+G.sub.ssK.sub.ss).sup.1G.sub.ssK.sub.ssr+(I+G.sub.ssK.sub.ss).sup.1v.sub.ss,(6)
u.sub.ss=(I+K.sub.ssG.sub.ss).sup.1r(I+K.sub.ssG.sub.ss).sup.1K.sub.ssv.sub.ss,(7)
where r .sup.m is the spatial excitation signal to be designed.
C. Spatial Optimal Input Design for the CD Process

(24) When it comes to the spatial optimal input, the parameters of interest for the CD process model will be those in the gain matrix G.sub.ss (or more specifically, the parameters in a column of G.sub.ss). Notice that the optimal input design directly based on the closed-loop model (6)-(7) is nontrival due to the large input-output dimensions as well as the large number of parameters in G.sub.ss. To avoid this issue, we propose to use a scalar transfer function along the spatial coordinate to represent the spatial response of the actuators. In this sense, the original optimal input design aimed for the MIMO CD model can be re-formulated into that for a scalar spatial model, which significantly reduces the associated complexity. However, the price to pay is that the scalar spatial transfer function will have to be non-causal as any bumped actuator will generate responses on two sides (see FIG. 5), analogous to the past and future for the conventional time coordinate. The following section will demonstrate the specific algorithms of converting the CD steady-state process model G.sub.ss and controller K.sub.ss into the scalar no causal then further into the causal transfer function form.

(25) II. Causal Scalar Transfer Function Representation of the CD Process

(26) In this section, we show the procedures to develop causal-equivalent models for both the CD closed-loop model and open-loop model. Let us focus on the steady-state process model G.sub.ss since the steady-state controller matrix K.sub.ss will follow the same vein as long as it possesses a similar structure as G.sub.ss.

(27) A. Non-causal Transfer Function Form of the Closed-loop CD Process

(28) From the aforementioned structure of the G.sub.ss as well as Assumption 1, one is readily able to extract a scalar non-causal FIR model from any single column of G.sub.ss to represent the spatial impulse response of the actuator,
g(,.sup.1)=g.sub.n.sup.n+. . . +g.sub.0+. . . +g.sub.n.sup.n,(8)
where is the spatially forward-shift operator. The positive and negative powers of denote the anti-causal and causal shift. The g.sub.i, i=n, . . . , n, are the spatial impulse response coefficients of each single actuator and in general the symmetry of the impulse response is enforced, i.e., g.sub.i=g.sub.i. As in most cases, the noncausal FIR model (8) will have a high order (i.e., n normally is large), a parsimonious non-causal transfer function is necessary to simplify this model. Before we show that, the following assumption is posed.

(29) Assumption 3. The CD MIMO steady-state models (e.g., G.sub.ss and K.sub.ss are Toeplitz-structured, and the corresponding spatial impulse response sequence satisfies the Wiener-Hopf factorization conditions: real, symmetric, and taking G.sub.ss as an example,
g.sub.n.sup.n+. . . +g.sub.0+. . . +g.sub.n.sup.n=M()M(.sup.1),,(9)
where =e.sup.j. Here M() has the following expression,
M()=m.sub.0+m.sub.1.sup.1+. . . +m.sub.n.sup.n,(10)
where m.sub.i, i=1, . . . n, are the coefficients.

(30) An immediate observation is that the frequency response of the left-hand side of (9) is non-negative and real for any frequency, which places certain restrictions on the scope of the possible spatial impulse response shapes that we may investigate. However, industrial experience reveals that most actual actuator response shapes are able to satisfy this condition. The relationship between G.sub.ss and K.sub.ss from (5) affirms that if G.sub.ss satisfies (5) then so does K.sub.ss.

(31) After obtaining the causal FIR model M(.sup.1), the next step would be to find a parsimonious transfer function model (e.g. output error model) to represent M(.sup.1). This process can be accomplished from the system identification toolbox in Matlab and the original noncausal g(, .sup.1) is re-written as follows,

(32) $\begin{array}{cc}\stackrel{\_}{g}\left(,{}^{-1}\right)=\frac{B\left(\right)B\left({}^{-1}\right)}{A\left(\right)A\left({}^{-1}\right)}& \left(11\right)\\ B\left(\right)=b_0+b_1{}^{-1}+\mathrm{.Math.}+b_{n_b}{}^{-n_b},& \left(12\right)\\ A\left(\right)=1+a_1{}^{-1}+\mathrm{.Math.}+a_{n_a}{}^{-n_a},& \left(13\right)\end{array}$
where n.sub.a and n.sub.b are the orders of B(.sup.1) and A(.sup.1), respectively. In a similar fashion, the noncausal transfer function form of the controller is achieved to be,

(33) $\begin{array}{cc}\stackrel{\_}{k}\left(,{}^{-1}\right)=\frac{F\left(\right)F\left({}^{-1}\right)}{E\left(\right)E\left({}^{-1}\right)}& \left(14\right)\\ F\left({}^{-1}\right)=f_0+f_1+\mathrm{.Math.}+f_{n_f}{}_f^{-n_f},& \left(15\right)\\ E\left({}^{-1}\right)=1+e_1{}^{-1}+\mathrm{.Math.}+e_{n_e}{}_e^{-n_e},& \left(16\right)\end{array}$
where n.sub.e and n.sub.f are the orders of E(.sup.1) and F(.sup.1), respectively. From (11)-(16), the original high-dimensional MIMO steady-state closed-loop model (6)-(7) can be replaced by scalar but non-causal transfer functions,

(34) $\begin{array}{cc}{y}_{\mathrm{ss}}\left(x\right)=\frac{\stackrel{\_}{g}}{1+\stackrel{\_}{g}\stackrel{\_}{k}}r\left(x\right)+\frac{1}{1+\stackrel{\_}{g}\stackrel{\_}{k}}{v}_{\mathrm{ss}}\left(x\right),& \left(17\right)\\ u_{\mathrm{ss}}\left(x\right)=\frac{\stackrel{\_}{g}}{1+\stackrel{\_}{g}\stackrel{\_}{k}}r\left(x\right)-\frac{\stackrel{\_}{k}}{1+\stackrel{\_}{g}\stackrel{\_}{k}}{v}_{\mathrm{ss}}\left(x\right),& \left(18\right)\end{array}$
where x stands for the spatial coordinate. Note that the input and output sensitivity functions have the same non-causal transfer function representation as shown in the above equations.
B. Causal Equivalent Closed-loop Models

(35) Up to now, the closed-loop scalar non-causal model of the CD process (17)-(18) is still not a convenient form for further processing such as the optimal input design. In this subsection, we will develop methods to find causal equivalent models for the non-causal transfer functions such as g(, .sup.1). First, the following Lemma is necessary.

(36) Lemma 1. Suppose that g.sub.1(, .sup.1) and g.sub.2(, .sup.1) satisfy the Wiener-Hopf factorization conditions. Then the sum g.sub.1(, .sup.1)+g.sub.2(, .sup.1) can also satisfy these conditions.

(37) Proof. Since g.sub.1(, .sup.1) and g.sub.2(, .sup.1) satisfy (9), we have,
g.sub.1(e.sup.j,e.sup.j)0,,
g.sub.2(e.sup.j,e.sup.j)0,,
Thus it follows that,
g.sub.1(e.sup.j,e.sup.j)+g.sub.2(e.sup.j,e.sup.j)0,.(19)
Besides, the coefficient sequence of (19) is real and symmetric. Thus one is always able to find an M() such that (9) is satisfied. This ends the proof.

(38) Defining

(39) $\stackrel{\_}{S}=\frac{1}{1+\stackrel{\_}{g}\stackrel{\_}{k}}$
from (17)-(18), we have,

(40) $\begin{array}{cc}\stackrel{\_}{S}=\frac{A\left(\right)A\left({}^{-1}\right)E\left(\right)E\left({}^{-1}\right)}{A\left(\right)A\left({}^{-1}\right)E\left(\right)E\left({}^{-1}\right)+B\left(\right)B\left({}^{-1}\right)F\left(\right)F\left({}^{-1}\right)}& \left(20\right)\end{array}$
From Lemma 1, it follows that the denominator of (20) can be factorized to be the product of a causal FIR filter and its anti-causal form. Therefore, the closed-loop transfer functions (17)-(18) are simplified to be,
y.sub.ss(x)=S.sub.1(,.sup.1)r(x)+S.sub.2(,.sup.1)v.sub.ss(x),(21)
u.sub.ss(x)=S.sub.2(,.sup.1)r(x)S.sub.3(,.sup.1)v.sub.ss(x),(22)
where S.sub.i(,.sup.1), i=1,2,3, has a similar structure as (11) and (14). Further notice that can also be represented by a noncausal transfer function as is assumed in the previous sections. In other words, the spatial noise v.sub.ss has the following expression,

(41) $\begin{array}{cc}{v}_{\mathrm{ss}}\left(x\right)=\stackrel{\_}{h}\left(,{}^{-1}\right)e\left(x\right)=\frac{D\left(\right)D\left({}^{-1}\right)}{C\left(\right)C\left({}^{-1}\right)},& \left(23\right)\end{array}$
where {e(x)} is a spatial white noise sequence. To find a causal-equivalent transfer function for (21)-(23), we establish the following theorem.

(42) Theorem 1. Consider a stochastic process with the output sequence {y(x), x=1, . . . , m} (in the sequel, the subscript is omitted and the argument x is used to indicate the steady-state and output sequence) generated according to the following non-causal Box-Jenkins model

(43) $\begin{array}{cc}y\left(x\right)=\frac{M\left(\right)M\left({}^{-1}\right)}{N\left(\right)N\left({}^{-1}\right)}r\left(x\right)+\frac{R\left(\right)R\left({}^{-1}\right)}{S\left(\right)S\left({}^{-1}\right)}e\left(x\right),& \left(24\right)\end{array}$
where {e(x), x=1, . . . , m} is a Gaussian white noise sequence. The polynomials with arguments .sup.1 and are the causal and anti-causal parts, respectively. Assume that all the polynomials have no zeros on the unit circle and are minimum phase. Then there exist causal polynomials {tilde over (M)}.sub.y(.sup.1), .sub.y(.sup.1), {tilde over (R)}.sub.y(.sup.1), S.sub.y(.sup.1) and a white noise sequence {.sub.y(x)} as well as a stochastic sequence {{tilde over (y)}(x)} which has the same spectra as {y(x)} such that,

(44) 0$\begin{array}{cc}\tilde{y}\left(x\right)=\frac{{\tilde{M}}_y\left({}^{-1}\right)}{{\tilde{N}}_y\left({}^{-1}\right)}r\left(x\right)+\frac{{\tilde{R}}_y\left({}^{-1}\right)}{{\tilde{S}}_y\left({}^{-1}\right)}{\tilde{e}}_y\left(x\right).& \left(25\right)\end{array}$

(45) Proof. If multiplying both sides of (24) by using N()N(.sup.1)S()S(.sup.1), we can obtain,
N()N(.sup.1)S()S(.sup.1)y(x)=M()M(.sup.1)S()S(.sup.1)r(x)+N()N(.sup.1)R()R(.sup.1)e(x),(26)
Define the roots of the anti-causal polynomials M(.sup.1), N(.sup.1), R(.sup.1), S(.sup.1) to be, respectively, .sub.i, .sub.i, .sub.i and .sub.i. Let

(46) ${}_M={}_i\frac{{}^{-1}-{}_i}{-{}_i},{}_N={}_i\frac{{}^{-1}-{}_i}{-{}_i},{}_R={}_i\frac{{}^{-1}-{}_i}{-{}_i},{}_S={}_i\frac{{}^{-1}-{}_i}{-{}_i}.$
Notice that N().sub.N=N.sup.2(.sup.1) and the same also holds for M(), R() and S(). Multiplying both sides of (26) by .sub.m.sub.s, after some manipulations, one can obtain,

(47) $\begin{array}{cc}{N}^2\left({}^{-1}\right)S^2\left({}^{-1}\right)\tilde{y}\left(x\right)=M^2\left({}^{-1}\right)S^2\left({}^{-1}\right)r\left(x\right)+R^2\left({}^{-1}\right){\tilde{e}}_y\left(x\right),& \left(27\right)\\ \mathrm{where}& \\ \tilde{y}\left(x\right)=\frac{{}_M}{{}_N}y\left(x\right),{\tilde{e}}_y\left(x\right)=\frac{{}_M{}_S}{{}_N{}_S}e\left(x\right)& \left(28\right)\end{array}$
Since .sub.M, .sub.N, .sub.R, and .sub.s are all-pass filters, {{tilde over (e)}.sub.y(x)} is a white noise sequence with the same spectra as {e(x)} but may correspond to different realizations. Besides, {{tilde over (y)}(x)} has the same spectra as {y.sub.ss(x)}. Therefore, (25) is verified by pairing {tilde over (M)}(.sup.1)=M.sup.2(.sup.1) and so on with (27), which ends this proof.

(48) Remark 1. From Theorem 1, one may interpret that the equivalence between {{tilde over (y)}(x)} and {y(x)} is in terms of the spectra, although the realizations might be different. However, this equivalence greatly facilitates the maximum likelihood estimation for the original noncausal model by reducing it into a causal-equivalent form. The rationale of performing in this way is based on the conclusion that the log-likelihood function for the non-causal model and the causal model converge to be the same with probability one as the sample number tends to infinity, which can also be extended to the non-causal Box-Jenkins model in (24).

(49) Similarly, the input signal u.sub.ss(x) in (22) can also be represented through causal filters,

(50) $\begin{array}{cc}\tilde{u}\left(x\right)=\frac{{\tilde{M}}_u\left({}^{-1}\right)}{{\tilde{N}}_u\left({}^{-1}\right)}r\left(x\right)+\frac{{\tilde{R}}_u\left({}^{-1}\right)}{{\tilde{S}}_u\left({}^{-1}\right)}{\tilde{e}}_u\left(x\right),& \left(29\right)\end{array}$
where {(x)} and {u.sub.ss(x)} have the same spectra. The equations (25) and (29) will be necessary for the optimal input design in the sequel.
C. Covariance Matrix Equivalence of the Causal and Non-causal Model Parameter Estimates

(51) It is well known that if the white noise is Gaussian distributed, the prediction error method with properly chosen criterion will coincide with the maximum likelihood estimation. It has been shown that, for the open-loop data, the log-likelihood function of the non-causal ARX model and that of the corresponding causal ARX model will converge to the same value as the sample number tends to infinity. In this subsection, we will demonstrate a similar statement for the closed-loop data.

(52) Theorem 2. Let us consider the following non-causal process model ( is the parameter in a compact set ),
y(x)=g(,.sup.1,)u(x)+h(,.sup.1,)e(x),(30)
where g was defined in (11)-(13) and h was defined in (23). e (x) is Gaussian white noise. Suppose that the data is generated in the closed-loop with controller model (14)-(16) and that all the relevant transfer functions are stable. Denote .sub.y(y) as the log-likelihood function of the non-causal model (30) and .sub.{tilde over (y)}({tilde over (y)}) as the log-likelihood function of the causal-equivalent model of (30) obtained similarly as (25). We have, as m.fwdarw.,

(53) ${\sup }_{}.Math.{}_y\left(y\right)-{}_{\tilde{y}}\left(\tilde{y}\right).Math..fwdarw.0,w.p\mathrm{.1},\underset{}{\sup }.Math.\frac{d{}_y\left(y\right)}{d}-\frac{d{}_{\tilde{y}}\left(\tilde{y}\right)}{d}.Math..fwdarw.0,w.p\mathrm{.1}$
Remark 2. Theorem 2 implies that both the log-likelihood function and its derivative with respect the parameter obtained from the original noncausal model and the causal model are identical asymptotically. Therefore, we can conclude that the parameter covariance matrices from the two schemes coincide, and hence we may perform the optimal input design based on the causal model.
III. Closed-Loop Optimal Input Design

(54) The optimal input design in the closed-loop for the steady-state CD process model was investigated. The emphasis will be placed on the non-causal CD process model due to the drawbacks associated with the MIMO CD closed-loop model (6)-(7). The actual implementation of the input design will be on the causal equivalent CD closed-loop model due to Theorem 2.

(55) Note that in practice, the noise model parameters are of less interest and thus we split the parameter to be =[.sup.T.sup.T].sup.T, where is the process model parameter vector and is the noise model parameter vector. For the input design, the focus will be on minimizing the covariance of only. Based on this motivation, due to Theorem 2, the parameter covariance matrix of , P.sub., is expressed as,

(56) $\left(31\right)$$P_{}{\frac{1}{m}\left[\frac{1}{2{}_0}{}_{-}^{}\frac{1}{{.Math.\tilde{h}\left(e^j,{}_0\right).Math.}^2}\frac{\tilde{g}\left(e^j,{}_0\right)}{}{}_{\tilde{u}}\left(\right)\frac{{\tilde{g}}^T\left(e^{-j},{}_0\right)}{}d\right]}^{-1},$
where .sub.0 is the variance of the noise {tilde over (e)}.sub.y(x). {tilde over (g)} and {tilde over (h)} are the causal equivalent forms of g and h, respectively. The input spectrum .sub.(), according to (29), is related to the excitation spectrum .sub.r()via,

(57) $\begin{array}{cc}{}_{\tilde{u}}\left(\right)={.Math.\frac{{\tilde{M}}_u\left(e^{-j}\right)}{{\tilde{N}}_u\left(e^{-j}\right)}.Math.}^2{}_r\left(\right)+{.Math.\frac{{\tilde{R}}_u\left(e^{-j}\right)}{{\tilde{S}}_u\left(e^{-j}\right)}.Math.}^2{}_0.& \left(32\right)\end{array}$

(58) The closed-loop optimal input design can be formulated as minimizing a function of the parameter covariance P.sub. subject to a set of constraints, e.g., input and output power constraints,

(59) $\begin{array}{cc}\underset{{}_r\left(\right)}{\min }{f}_0\left(P_{}\left({}_r\left(\right)\right)\right)& \left(33\right)\\ s.t.\frac{1}{2}{}_{-}^{}{}_u\left(\right)d{c}_u,& \left(34\right)\\ \frac{1}{2}{}_{-}^{}{}_y\left(\right)d{c}_y,& \left(35\right)\end{array}$

(60) where c.sub.u and c.sub.y are the limits on the input signal power and output signal power. The constraints (34)-(35) can be written in terms of the design variable .sub.r() by (32) and (25), respectively. As this optimization problem is still infinite-dimensional (since .sub.r() is a continuous function of ), a technique known as the finite dimensional parameterization can be employed to reduce it into finite-dimensional case. Specifically, .sub.r() can be parameterized by the definition of a spectrum,

(61) $\begin{array}{cc}{}_r\left(\right)=\underset{k=-m_c}{\overset{m_c}{.Math.}}{c}_k{e}^{-jk}0,,& \left(36\right)\end{array}$
where c.sub.k, k=m.sub.c, . . . , m.sub.c, are the parameters, and m.sub.c is the selected number of parameters. With (36), the original optimization problem can be cast into one with finite number of parameters. It is worth pointing out that the non-negativity of the parameterized spectrum (36) at any frequency has to be satisfied while searching for the optimal c.sub.k. This requirement is fulfilled by the KYP lemma by constructing a controllable and observable state-space realization for the spectrum. After these treatments, we will obtain a neat and solvable convex optimization (choose f.sub.0(.Math.) to be convex) problem which can be readily solved by the off-the-shelf solvers such as the CVX toolbox.
Remark 3. Note that the aforementioned optimal input design only considers the power constraints on the input and output (34)-(35). However, in practice, the hard constraints on the CVs and MVs make more sense and this is still an open problem for the frequency-domain optimal input design as posed above. Besides, specific to the CD process, the second-order bending constraints preventing the picketing on the actuators are also important.
IV. Case Study

(62) A simulation example was used to validate the proposed CD process model identification and closed-loop optimal input design methods. In particular, the effect of the optimally designed input on the identification was compared with that of the bumped excitation that is currently employed in the industry as described in Chu et al. U.S. Pat. No. 8,224,476 to Chu et al., which is incorporated herein.

(63) In practice, the spatial response shape of a single actuator is assumed to satisfy the following nonlinear equation,

(64) $\begin{array}{cc}{g}_s\left(x\right)=\frac{y}{2}\left\{e^{\frac{{a\left(x+\right)}^2}{{}^2}}\cos \frac{\left(x+\right)}{}+e^{-\frac{{a\left(x-\right)}^2}{{}^2}}\cos \frac{\left(x+\right)}{}\right\}& \left(37\right)\end{array}$
where , , , represent the gain, width, divergence and attenuation, respectively. x is the spatial coordinate. In this example, these parameters are specified with values, respectively, =0.3802, =268.6414 mm, =0.10, =3.5. The response shape under impulse signal of amplitude 5 is illustrated as the solid curve in the plot of FIG. 5A. For convenience, we assume that the CD process has 222 actuators and measurement bins. The controller is chosen to be CD-MPC with prediction horizon 25 samples and control horizon 1 sample (sampling interval is 12 seconds). The weighting matrices in the cost function is selected to be Q.sub.1=0.4I, and Q.sub.3=0.1667I. The .sub.K in (5) is computed to be 12.3212. From the above discussion of the development of the causal-equivalent models, one is able to obtain a non-causal scalar model for the CD process and the controller, respectively. The impulse response curves of these non-causal models are shown in FIGS. 5A and 5B in the dash-dotted curves. Noting that here for simplicity, we select n.sub.b=n.sub.f=1, n.sub.a=n.sub.e=2 and higher orders will improve the fitting performance but will also increase the computation of the optimal input design. The noise variance is chosen to be 0.1 with noise model =l (output-error structure). Notice that the small process gain (the causal-equivalent model has even smaller gain) requires large excitation signal for a good signal-to-noise ratio.

(65) To make the comparison between the optimal excitation signals with the bumped excitation meaningful, we set a hard constraint 10 on the amplitude of the excitation signals. For the optimally designed input, if any part of its amplitude violates this constraint, we will set that part to be saturated on this bound. For the bumped signal, the amplitude of the bumps alternate between 10 and 10. To achieve this goal, we carefully choose c.sub.u=4 and c.sub.y=0.2. The plot in FIG. 6A shows the impulse responses of the estimated models in 100 Monte-Carlo simulations. The plot in FIG. 6B shows the corresponding results from another 100 simulations under bumped excitation signal. From FIGS. 6A and 6B, the impulse responses of the estimated models from the bumped signal have large variance and the estimated gains are not precise. However, the estimates from the optimal input demonstrate smaller variance and the gains are more accurate. Therefore, the designed optimal input outperforms the bumped excitation signals in this example in terms of the closed-loop identification performance.

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