Low-order ODE approximations and receding horizon control of surface roughness during thin-film growth

Chemical Engineering Science 63 (2008) 1246 – 1260 www.elsevier.com/locate/ces

Low-order ODE approximations and receding horizon control of surface roughness during thin-film growth Amit Varshney, Antonios Armaou ∗ Department of Chemical Engineering, The Pennsylvania State University, University Park, PA 16802, USA Available online 2 August 2007

Abstract The problem of feedback controller synthesis with objective to control the microstructure during thin-film growth is considered. The problem of the non-availability of closed form dynamic models for the evolution of the microstructure is circumvented by deriving low-order statespace models that approximate the underlying kinetic Monte Carlo simulations. Initially, a finite set of “coarse” observables is identified from spatial correlation functions to represent the coarse microscopic state and capture the dominant characteristics of the microstructure during the deposition process. Subsequently, a state-space model is identified, employing proper orthogonal decomposition and Carleman linearization, that describes the evolution of the coarse observables. The state-space model is subsequently employed to design receding horizon controllers that regulate the surface roughness of the thin-film at a specified set-point during the growth process by manipulating the substrate temperature. The above approach is applied to: (i) a deposition process modeled using solid-on-solid model on a one-dimensional lattice; and (ii) an anisotropic deposition process on a two-dimensional lattice. Closed-loop simulations at various growth rates and in the presence of disturbances are performed to demonstrate the effectiveness of the proposed controller design scheme. 䉷 2007 Elsevier Ltd. All rights reserved. Keywords: Thin-film growth; Model-reduction; Receding horizon control; Surface roughness; System identification; Multiscale modeling

1. Introduction Thin-film deposition is an important process in semiconductor manufacturing. For example, an integrated circuit consists of several layers of chemically altered and patterned thin-films in order to achieve desired electrical characteristics (Edgar et al., 2000; Badgwell et al., 1995). The performance of such devices depends not only on the sharpness of the patterns and the interfaces between the layers of these films, but also heavily on the resulting microstructure of the film material. Driven by a highly competitive global market and reduced profit margins, recent focus has been towards process design, operation and control with tighter processing and device performance specifications, which translate into restrictions on pattern definition and microstructure deviation from the desired ones. These trends have motivated significant research efforts on model-based optimization and feedback control of thin-film deposition processes with emphasis to control thin-film spatial ∗ Corresponding author. Tel.: +1 8148655316; fax: +1 8148657846.

E-mail address: [email protected] (A. Armaou). 0009-2509/$ - see front matter 䉷 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.ces.2007.07.058

nonuniformity and its microstructure. The necessary process models are generally derived by combining continuum conservation laws, which are in the form of distributed nonlinear partial differential equations (PDEs), to describe reactor scale (macroscopic) phenomena with atomistic evolution rules, such as kinetic Monte Carlo (kMC) or molecular dynamics, to describe feature scale (microscopic) phenomena. Previous studies that address the issue of control of (macroscopic) spatial uniformity of thin-films include (Theodoropoulou et al., 1998; Park et al., 1999; Emami-Naeini et al., 2003) for rapid thermal processing (Armaou et al., 2001; Armaou and Christofides, 1999) for plasma-assisted etching and plasma-assisted deposition and (Banks et al., 2002; Varshney and Armaou, 2005; Varshney and Armaou, 2006b) for chemical vapor deposition. To synthesize practically implementable controllers these studies rely on the reduction of the underlying macroscopic process model since the detailed models are computationally expensive. Such an order reduction is possible because the underlying macroscopic process models are highly dissipative, i.e., their dominant dynamics are low dimensional (Christofides and Daoutidis, 1997; Christofides, 2001; Armaou and Christofides, 2002).

A. Varshney, A. Armaou / Chemical Engineering Science 63 (2008) 1246 – 1260

The issue of control of film properties such as micro-structure and composition is also significant since they define the film electromechanical properties (Akiyama et al., 2002). The possibility of obtaining real-time measurements of microscopic properties using scanning tunneling microscopy (Voigtlander, 2001) and spectroscopic ellipsometry (Zapien et al., 2001) has motivated research in real-time control of the process. Unlike continuum models for reactor scale phenomena, atomistic models (such as kMC) that describe the evolution of microscopic properties of the film are unavailable in closed form which limits their applicability towards model-based controller synthesis. Results that address the issue of control of the microstructure of thin-films include (Lou and Christofides, 2003) where kMC simulations were employed as the underlying process model and (Ni and Christofides, 2005; Lou and Christofides, 2005) which employ linear stochastic PDEs as the underlying growth model to enable the use of proportional-integral and predictive controllers in lieu of measurements. In Gallivan and Murray (2004) a methodology for reduced order modeling was developed through statistical reduction of the master equation. However, application of the above approaches for more complex deposition processes is not a straightforward task, for instance, a linear SPDE may not be able to accurately describe a complex deposition process. An extension of the SPDE-based approach to control surface roughness using a nonlinear feedback control scheme has been recently published (Lou and Christofides, 2006). However, all of the previous approaches do not explicitly consider state and input constraints which are critical to semiconductor industry owing to increasingly higher product quality specifications and productivity requirements combined with economical considerations. Motivated by the above, an alternative approach is proposed in the current work. Based on the methodology presented in Varshney and Armaou (2006a), we identify a finite set of coarse observables describing the dominant traits of the microstructure of the thin-film during deposition. These coarse variables are derived from the low-order spatial correlation functions of the film surface and characterize the surface uniquely in the sense that the (unknown) coarse variable dynamics are observable, i.e., surfaces which are statistically identical with respect to these variables exhibit identical coarse dynamical behavior. Consequently, these variables are conceived as the coarse realization of the microscopic state. Subsequently, we employ the above coarse observables to develop low-order approximate (discrete time) state-space models, which can be employed as an alternative to kMC simulations. To this end, initially we employ proper orthogonal decomposition (POD) to obtain an optimal representation of the coarse data. Subsequently, we approximate the unidentified system dynamics using polynomial expressions and reformulate the approximate polynomial system as a bilinear system through Carleman linearization. The parameters of the bilinear system are first estimated offline through the subspace-based algorithm N4SID (Favoreel et al., 1999; Favoreel, 1999) and further refined through solution of a series of nonlinear dynamic optimization programs which minimize the prediction error between the model and the kMC simulations in least-squares sense. The developed nonlinear model is

1247

subsequently utilized to design a receding horizon controller to control the surface roughness during the thin-film growth process. An advantage of receding horizon control is that state and input constraints can be readily incorporated into the control scheme. The reader is referred to Mayne et al. (2000) and Morari and Lee (1999) for excellent reviews on receding horizon control. 2. Mathematical preliminaries We seek to construct discrete-time nonlinear dynamic models for thin-film growth with the following state-space representation: z˙ = f (z) +

m 

gj (z)uj ,

z(0) = z0 ,

t = 0, 1, . . . ,

j =1

ym = Cz,

yci = hi (z),

i = 1, . . . , m,

(1)

where z(t) ∈ Rn is the vector of state variables, u(t) ∈ Rm is the vector of manipulated variables, ym (t) ∈ Rl is the measured variables vector, and yci (t) is the ith control output. f (z(t)): Rn → Rn and gj (z(t)): Rn → Rn are nonlinear vector functions of state which are sufficiently smooth with respect to their arguments. We assume that the functions exist but are unavailable. C is the measurement sensor shape function, and hi is a nonlinear function of state representing ith control objective. Without loss of generality, the origin is assumed to be a steady-state of the system. The Kronecker product between matrices A ∈ RN×M and B ∈ RL×K is defined as the matrix C ∈ RLN×KM , ⎡ A B A B ··· A B ⎤ 1,1

1,2

1,M

⎢ A2,1 B ⎢ C=A⊗B ≡⎢ ⎢ .. ⎣ .

A2,2 B .. .

···

AN,1 B

AN,2 B

· · · AN,M B

..

.

A2,M B .. .

⎥ ⎥ ⎥. ⎥ ⎦

Also, the kth order Kronecker product is recursively defined as A[k] = A ⊗ A[k−1] , ∀k 1, A[0] = 1. The symbol  is used to denote the Khatri–Rao product, which is a column-wise Kronecker product for two matrices with an equal number of columns such that given two matrices A = [a1 , a2 , . . . , aN ] ∈ RI ×N and B =[b1 , b2 , . . . , bN ] ∈ RJ ×N , the Khatri–Rao product is equal to a matrix C ∈ RI J ×N defined as C = A  B = [a1 ⊗ b1 , a2 ⊗ b2 , . . . , aN ⊗ bN ].

(2)

The orthogonal projection of the row space of matrix A ∈ RN×M into the row space of matrix B ∈ RL×M is defined as A/B = AB † B. The Moore–Penrose pseudoinverse of a matrix A ∈ RN×M is a matrix A† ∈ RM×N such that AA† A = A. A function : R  0 ×Z  0 → R  0 is said to belong to class-KL if it is nondecreasing in its first argument, nonincreasing in its second argument, and lims→0+ (s, k) = limk→∞ (s, k) = 0 (Teel, 2004). The sets R  0 and Z  0 denote the sets of nonnegative real numbers and nonnegative integers, respectively. A function from R  0 to R  0 is said to belong to class-G if it

1248

A. Varshney, A. Armaou / Chemical Engineering Science 63 (2008) 1246 – 1260

is continuous, zero at zero and nondecreasing. Also we define |x|S := inf z∈S |x − z|, where S ⊂ Rn is compact. 3. System identification procedure The microscopic state of the thin-film during homoepitaxial growth using kMC simulations is a vector comprising of lattice heights at various lattice locations. Without loss of generality, we consider a square lattice L := {(lx , ly ): lx , ly ∈ [1, 2, . . . , N]} of size NT = N × N . Xit , Xt ∈ Z  0 , i ∈ L denote the height at lattice location i at any time instant t; the t ]T defines the microscopic state vector Xt = [X1t , X2t , . . . , XN T at any time instant t. However, it is not practical to employ X as state in Eq. (1) owing to its high dimensionality (e.g., for a simple growth model lattices of at least 50 × 50 size are usually recommended). To address this problem, we identify a relatively small set of “coarse” variables that represent X in a statistical sense, and can be conveniently employed as the state in Eq. (1). Subsequently, we aim to obtain the estimates of the unknown nonlinear functions f, g1 , . . . , gm that describe the dynamics of coarse state. We note that the above identification problem is truly a black-box system identification problem, since no assumptions can be made about the functions a priori. The identification is performed offline based on the data obtained from kMC simulations. We first identify a set of coarse variables based on the method presented in Varshney and Armaou (2006a). Next we employ POD to analyze the multidimensional coarse data. This essentially provides an orthonormal basis for representing the available data in a certain least-squares optimal sense. Galerkin projection onto the subspace spanned by the dominant POD modes results in a compact (often of reduced dimension) representation of the original data. Subsequently, we employ Carleman linearization to transform the unknown nonlinear system into an extended bilinear system whose parameters are estimated using subspace method N4SID followed by least-squares minimization such that the error between the kMC data and model prediction is minimized. N4SID is an noniterative algorithm, thus circumvents the convergence problem common to most prediction error methods. Furthermore, it is fast and numerically stable since it is based on linear algebra techniques. The proposed system identification approach comprises of the following steps: 1. Identify the coarse microscopic state from appropriate correlation functions that statistically characterize the deposition surface. 2. Reduce the state-space dimension by obtaining the minimal basis for coarse state using POD. 3. Obtain a bilinear representation of the unknown nonlinear state-space model using Carleman linearization. 4. Compute the first estimates of the system matrices of the bilinear model using subspace identification algorithm N4SID. 5. Refine the initial estimates of the system matrices by nonlinear least-squares minimization. The steps are elaborated on in the following subsections.

3.1. Identification of coarse observables The sequence Xit , i ∈ L, which defines the microstructure at time t, can also be considered as a stochastic process which is completely characterized by the joint distribution function (JDF) (Rozman and Utz, 2002; Sobczyk and Kirkner, 2001; Torquato, 2001) given by t FXt ,...,Xt (x1 , . . . , xNT ) = P (X1t x1 , . . . , XN xNT ). T 1

NT

(3)

In practice, however, it is rarely possible, and seldom required, to know the entire JDF and sufficient information about the stochastic process can be obtained from its low-order statistics. It can be argued that the random deposition processes result in isotropic microstructures; this in turn implies that the corresponding stochastic processes, Xt , are strict sense stationary, i.e., the statistical properties are invariant under translation of the origin. Mathematically, it is expressed as FXt ,...,Xt (x1 , . . . , xNT ) = FXt 1

NT

t 1+c ,...,XNT +c

(x1 , . . . , xNT ),

(4)

where c is a positive constant. In order to characterize the structure of the deposition surface, we begin by obtaining an approximation of the first-order statistics of X t in the form of the probability density function (PDF). The unknown PDF can be approximated using series expansion around a known distribution (x), as stated in the following theorem (Hildebrandt, 1931): Theorem 3.1 (Charlier’s Theorem for Series B). Any real valued function F (x) which vanishes for x = +∞ and x = −∞ may be formally expanded in terms of another function (x) and its successive differences in the form F (x) = lim (B0 (x) + B1 (x) n→∞

+ B2 2 (x)+ · · · +Bn n (x)), n (x) = n−1 (x)−n−1 (x − h),

0 (x) = (x)

(5)

provided that (x) possesses the following properties: 3.1.1. (x) and its differences are defined for all real values of x; 3.1.2. (x) and its differences vanish for x = ±∞; 3.1.3. lim x m n (x) = 0 for all real values of m and n; ±∞ +∞ 3.1.4. −∞ (x)  = 0. In this work we employ the Poisson distribution as our expansion basis, (x), given as ⎧ −m x ⎨e m if x 0, (x) = x! ⎩ 0 otherwise. Note that the Poisson distribution satisfies assumptions 2.1.1–2.1.4. The number of terms in the expansion are chosen according to the accuracy required, and in general, will depend upon the skewness of the distribution. The coefficients

A. Varshney, A. Armaou / Chemical Engineering Science 63 (2008) 1246 – 1260

1249

B0 , B1 , . . . , Bn can be computed using the method of moments (Aroian, 1937; Boas, 1949). The method is based on the following identities:   [B0 (x) + · · · + Bn n (x)] = F (x) = 1,   x[B0 (x) + · · · + Bn n (x)] = xF (x) = 1 ,

level. Hence, we require the coarse variables to be chosen such that they capture the dominant statistics of the microstructure and the rest of the unaccounted statistics to quickly become functionals of the coarse variables. A systematic approach to compute coarse variables is presented in Varshney and Armaou (2006a).

.. . 

Remark 1. It should be noted that the correlation functions such as (r), saturate for r > r0 due to the finite length of lateral interactions. For cases when r0 is small, a finite number of correlation values corresponding to the discrete positions r < r0 are sufficient to accurately describe the microstructure. If long range interactions are present, parameterization of the correlation functions may be considered. For example, for random surfaces, the ACF can be approximated as (Fenner, 2004)

x n [B0 (x) + · · · + Bn n (x)] =



x n F (x) = n ,

(6)

where 1 , . . . , n denote the first-order moments of the process. We note that Eq. (6) describes a system of n+1 linear equations from which the unknown coefficients B0 , B1 , . . . , Bn can be easily computed. In the absence of correlations between different lattice locations, the PDF is sufficient for the characterization of deposition surfaces. However, when lateral interactions are also present, information regarding the neighborhood of a lattice site is also required (i.e., lattice heights of neighboring lattice sites) for the characterization of deposition surfaces. We employ various two-point spatial correlation functions of the surface to characterize this information, since it is difficult to compute the corresponding JDF. Examples include the two-point autocorrelation function (ACF) defined as (i, j, t) = (Xit − Xit )(Xjt − Xjt ) , = Xit Xjt − Xit Xjt

i, j ∈ L

(7) (8)

and the height difference correlation function (HDCF) defined as, (i, j, t) = |Xit − Xjt | ,

i, j ∈ L,

(9)

where · represents the expectation. For statistically isotropic microstructures, (i, j, t) and (i, j, t) depend only on the distance r = |i − j | between two lattice points and the above expressions can be written as (i, j, t) = (r, t), (i, j, t) = (r, t), r = |i − j |.

(10)

For random surfaces the behavior of ACF can also be described in terms of the fractal dimension and Hurst coefficient of the surface (Fenner, 2004; Gneiting and Schlather, 2004) One can also define m-point correlation functions as Xit1 Xit2 . . . Xitm and utilize them to characterize the stochastic process Xi , if information contained in the two-point ACF is insufficient. However, in practice, m-point correlation functions with m 3 are cumbersome to work with since they result in m − 1 dimensional hyper-surfaces. One can, nevertheless, employ a number of two-point correlation functions to capture increasingly greater amount of spectral information. A discussion on various two-point correlation functions can be found in Torquato (2002). The selection of spatial correlations that serve as coarse variables is an ill-posed problem. This is due to the mismatch between degrees of freedom at microscopic level and the coarse

(r) = (0) exp(−(r/)2H ), where H and  are the Hurst coefficient and the correlation length, respectively. 3.2. POD to obtain optimal data representation The obtained coarse microscopic state is not necessarily minimal (due to overlapping information contained with the two-point correlation functions) which motivates the computation of dominant modes in the data set obtained from correlation functions. POD (also known as Karhunen–Loève decomposition or principle component analysis) provides a method for computing the optimal subspace to represent a given data set (Rathinam and Petzold, 2003). In addition to being optimal, POD has the property that it uses a modal decomposition that is completely data dependent and does not rely on the knowledge of process that generates the data. The method has been employed extensively in the context of both finite and infinite dimensional systems for data analysis, model reduction, optimization and control (Armaou and Christofides, 2002; Shvartsman and Kevrekidis, 1998; Banks et al., 2002) and is included for completeness. Let, z(t), z ∈ Rn , t ∈ [0, T ], denote a trajectory of process data (coarse variables, in the current case). Let, E(z, N, [0, T ]) := {z1 (t), z2 (t), . . . , zN (t)} denote a collection (or ensemble) of trajectories. We define the inner product between two ensembles, E(z1 , N, [0, T ]) and E(z2 , N, [0, T ]) as, (E(z1 ), E(z2 )) =

N  i=1

T 0

z1i (t)T z2i (t) dt.

(11)

Given an ensemble of trajectories, obtained from simulations or experiments, POD seeks to identify a subspace S ⊆ Rn , such that following total squared distance is minimized, E(z) − S E(z) =

N  i=1

T

(zi (t)

0

− S zi (t))T (zi (t) − S zi (t)) dt,

(12)

1250

A. Varshney, A. Armaou / Chemical Engineering Science 63 (2008) 1246 – 1260

where  ·  denotes the norm induced by Eq. (11) and S zi (t) denotes the orthogonal projection of zi (t) onto subspace S. The above minimization problem can be reformulated as an eigenvalue problem for the following covariance matrix: R=

N T  i=1

zi (t)zi (t)T dt.

(13)

we obtain, f¯(¯z) =

 ∞  1 ¯  jf[i]   i!

z¯ [i] ,

z¯ =0

i=1

 ∞   1  g¯ j (¯z) = g¯ j (0) + jg¯j [i]   i!

(17)

z¯ =0

i=1

0

z¯ [i] ,

where jf¯[i] |z¯ =0 ∈ Rk×k and jg¯j [i] |z¯ =0 ∈ Rk×k are the ith partial derivatives of f¯(¯z) and g¯ j (¯z), respectively, evaluated at z¯ = 0. For notational simplicity we denote Ai = (1/ i!)jf¯[i] |z¯ =0 , Bj i =(1/ i!)jg¯j [i] |z¯ =0 and Bj 0 = g¯j (0) for the rest of the article. With z¯ [i] , we denote ith Kronecker product of vector z¯ and (·)! represents the standard factorial of an integer. Thus, Eq. (15) can be equivalently written in the following form: i

The minimum value of E(z) − S E(z) over all k(n) is given by ni=k+1 i , where 1  2  · · ·  n are the eigenvalues of R. The orthogonal projection onto the subspace S is then given by S ≡ [1 , 2 , . . . , k ]T ,

(14)

where i is the eigenvector of R corresponding to eigenvalue i . Note that by construction R is symmetric and positive semidefinite which implies that i > 0, ∀i. Moreover, the linear subspace S is an invariant subspace (Holmes et al., 1996). If of the covariance matrix are such that nthe eigenvalues n /  , where is a small parameter, for some i=k+1 i i=1 i k < n then a lower dimensional representation of the original data can be obtained from x = S E(z). 3.3. Carleman linearization

m 

g¯ j (¯z)uj , z¯ (0) = S z0 ,

(15)

j =1

where f¯: Rk → Rk and g¯ j : Rk → Rk are defined as f¯(¯z) = S f (TS z¯ ), g¯ j (¯z) = S gj (TS z¯ ).

(16)

If k < n, then Eq. (15) represents a reduced order model which approximates the original system given by Eq. (1). Referring back to system defined by Eq. (15), applying McLaurin series expansion to vector fields f¯(¯z) and g¯ j (¯z),

m 

g¯ j (¯z)uj =

j =1

∞ 

Ai z¯ [i] +

∞ m  

Bj i z¯ [i] uj . (18)

j =1 i=0

i=1

We focus on finite order polynomial approximations of f¯(¯z) and g¯ j (¯z) obtained by truncating the corresponding infinite series upto order pf and pg , respectively. Without any loss of generality we can assume pf = pg + 1 = p, in which case Eq. (18) take the following form: z˙¯ ≈

Following the identification of the minimal coarse state, in this section we discuss the methodology adopted in this work for the identification of the unknown nonlinear dynamic system describing the evolution of coarse variables. We initially approximate the unidentified system dynamics, Eq. (1), using polynomial expressions and reformulate the approximate polynomial system as a bilinear system through Carleman linearization. At this point, it is important to note that even though kMC simulations result in noisy process data, we are still interested in the expected process behavior which is (approximately) obtained from averaging the kMC data generated from several independent runs. Let S ⊆ Rn be the k dimensional approximating subspace obtained from POD and S be the orthonormal projection onto S. If z¯ (t) denotes the projection of z(t) onto S, i.e., z¯ = S z, then from Eq. (1), z˙¯ = f¯(¯z) +

z˙¯ = f¯(¯z)+

i

p 

Ai z¯ [i] +

m p−1  

Bj i z¯ [i] uj .

(19)

j =1 i=0

i=1

To linearize the system of equation (19), we compute the dynamic behavior of z¯ [k] as follows: p−i p−i+1 m    d¯z[i] [l+i−1] = Ai,l z¯ + Bj i,l z¯ [l+i−1] uj , dt

(20)

j =1 l=0

l=1

q i−1+q and Bj i,l is defined simwhere Ai,l = i−1 q=0 In ⊗ Al ⊗ In ilarly. In denotes an unitary matrix of dimension n. T T Defining x = [¯zT z¯ [2] . . . z¯ [p] ]T , the system of equation (20) can be written in the following bilinear form: x˙ = Ax +

m 

[Bj xuj + Bj 0 uj ],

(21)

j =1

where A, Bj and Bj 0 are matrices of ⎤ ⎡ A1,1 A1,2 · · · A1,p ⎢ 0 A2,1 · · · A2,p−1 ⎥ ⎥ ⎢ ⎥ A=⎢ .. .. ⎥ , ⎢ .. .. ⎣ . . . . ⎦ 0 0 · · · Ap,1 ⎡B j 1,1 Bj 1,2 · · · Bj 1,p−1 ⎢ Bj 2,0 Bj 2,1 · · · Bj 2,p−2 ⎢ ⎢ 0 Bj 3,0 · · · Bj 3,p−3 Bj = ⎢ ⎢ ⎢ . .. .. .. ⎣ .. . . . 0

0

···

Bjp,0

the form:

0⎤ 0⎥ ⎥ ⎥ 0⎥ , ⎥ .. ⎥ .⎦ 0

A. Varshney, A. Armaou / Chemical Engineering Science 63 (2008) 1246 – 1260

⎡B

j 1,0

⎤

1251

white inputs (Favoreel et al., 1999). A modification for bilinear systems excited by general inputs has also appeared (Favoreel, 1999; Verdult and Verhaegen, 2001). In the current work, the N4SID algorithm of Favoreel (1999) is employed to obtain the initial estimates of matrices A, B and N. The algorithm for white inputs is briefly outlined in the following:

⎢ 0 ⎥ ⎥ ⎢ ⎥ ⎢ ⎥. ⎢ 0 Bj 0 = ⎢ ⎥ ⎢ .. ⎥ ⎣ . ⎦ 0 The above operation is also known as Carleman linearization (Kowalski and Steeb, 1991; Svoronos et al., 1994; Armaou, 2005) and allows us to transform a nonlinear system into an augmented linear system albeit at some loss of accuracy. Defining N = [B1 , B2 , . . . , Bm ] and B = [B10 , B20 , . . . , Bm0 ] the system of equation (21) can be written in the following compact form: x˙ = Ax + Nu ⊗ x + Bu.

(22)

Time-discretization of the system of equation (22) with piecewise constant control action, u(t) = u(kT ), kT t < (k + 1)T , where T is the sample time, leads to x((k + 1)T ) = Ax(kT ) +

∞ 

Ni u[i] ⊗ x(kT )

i=1

+

∞ 

Bi u[i] (kT ),

(23)

i=1

where the matrices A, Ni , and Bi are defined as A = eAT ,

T

eA(T −1 ) N Ni = 0

T

eA(1 −2 ) N . . .

0

T

eA(i−1 −i )

0

T

Bi = 0

eA(T −1 ) N

T

eA(1 −2 ) N . . .

0

x(k + 1) = Ax(k) + N u⊗ (k) ⊗ x(k) + Bu⊗ (k), y(k) = Cx(k). 1. Compute the projections:   Yi−1|0 , Oi = Yi|2i−1 Ui−1|0   Yi|0 , Oi+1 = Yi+1|2i−1 Ui|0

(27)

where Yp|q =[Yp Yp−1|q Up Yp−1|q ]T , Yq|p =[Yq Yq+1|p Uq  Yq+1|p ]T , Yp|p = Yp , Yp = [y(p)y(p + 1) . . . y(p + j )] for p > q and j  p, q. Up is defined similarly to Yp . The parameter i should be chosen such that p−1 di−1 = i−1 l > n, where l is the number of p=1 (m + 1) outputs. 2. Compute the singular value decomposition of Oi :  T  V1 S1 0 Oi = (U1 U2 ) 0 S2 V2T = i Xˆ i , 1/2

T

(26)

(28)

where S1 ∈ Rn×n , and i and Xˆ i are defined as

× NeAi di di−1 . . . d1 ,

0. The objective is to identify the matrices A, B, N, C of the observed system

i = U1 S1 , eA(i−1 −i )

0

× B di di−1 . . . d1 .

(24)

Without loss of generality we assume T =1. Truncating the infiT T nite series up to order q, and defining u⊗ = [uT u[1] . . . u[q] ]T , N = [N1 , N2 , . . . , Nq ] and B = [B1 , B2 , . . . , Bq ], we obtain the following bilinear discrete time system: x(k + 1) = Ax(k) + N u⊗ (k) ⊗ x(k) + Bu⊗ (k).

(25)

3.4. Bilinear system identification using N4SID and least-squares The aim of our system identification efforts is now recast as the one to estimate the unknown matrices A, B and N such that the error between kMC data and model prediction is minimized. For linear systems, prediction error methods (Ljung, 1987) and subspace identification algorithms such as N4SID (Favoreel et al., 2000) are commonly employed for system identification. Recently, N4SID algorithm was also extended for identification of bilinear state-space systems excited by

1/2 Xˆ i = S1 V1 .

(29)

3. Calculate the estimated state vector sequence at time step i + 1 as Xˆ i+1 = †i−1 Oi+1 ,

i−1 = [ Tl+1 , Tl+2 , . . . ]T .

(30)

where 1 , 2 , . . . , are the rows of i . 4. Obtain the system matrices by solving the following linear least-squares problem: ⎞ ⎛ Xˆ i      ˆ w A N B ⎜ Xi+1 ⎟ = . (31) ⎝ Ui  Xˆ i ⎠ + C 0 0 Yi|i v Ui Known parameterizations of the matrices A, B, N, and C can be easily handled by solving the above problem as a constrained linear least-squares problem. The estimates of A, B and N obtained from subspace identification are further refined through optimization using least-squares. Using kMC simulations, we initially generate an ensemble of trajectories E(x, ¯ M, [0, T ], x(0), u) from a variety of random initial conditions, x(0) ∈ Rk×M , and input actuation, u ∈ Rm×M . Subsequently, A, B, and N are estimated through

1252

A. Varshney, A. Armaou / Chemical Engineering Science 63 (2008) 1246 – 1260

optimization such that the error between model prediction and kMC is minimized. Mathematically, the above problem can be formulated as a nonlinear program (NLP), whose solution can be obtained using standard NLP algorithms such as successive quadratic programming (SQP). The NLP is formulated as min E(x, ¯ M, [0, T ], x(0), u) − E(x, M, [0, T ], x(0), u)

A,B,N

s.t. x i (t + 1) = Ax i (t) + N ui (t) ⊗ x i (t) + Bui (t), x i (0) = x i (0),

i = 1, . . . , M,

(32)

where E(x, M, [0, T ], x(0), u) denotes the ensemble of trajectories obtained from the approximate model. During ensemble construction the manipulated variable should span the entire range space in order to correctly sample all expected process behavior. Also, care must be taken so that the rate of change of manipulated inputs is constrained so that they do not excite the fast dynamics of the system. 4. Receding horizon controller synthesis

(33)

where x ∈ Rn , u ∈ Rm , and the dimensions n and m are different from the variables n and m defined previously. Let (t, x, u) where u = (u(0), u(1), . . .) be the solutions of Eq. (33) with initial condition x such that (0, x, u) = x. Receding horizon control is an optimization-based approach to feedback stabilization. The interested reader may refer to Henson (1998), and Mayne et al. (2000) for excellent reviews on receding horizon control. The open-loop optimal control problem, solved repeatedly at each instant the measurements become available, is formulated as min uM

M−1 

l((i, x, uM ), u(i)) + g((M, x, uM ))

i=0

s.t. x(t + 1) = Ax(t) + N u(t) ⊗ x(t) + Bu(t), u(·) ∈ U,

(34)

where M denotes the length of the prediction horizon or the output horizon, u = (u(0), u(1), . . .), u(M − 1), g: Rn → R  0 is the terminal cost and l: Rn × Rm → R  0 is the stage cost. U denotes the set of admissible input values which is assumed to be compact. In the most simple form, the set U may represent the upper and lower bounds on the manipulated input, such that, U := {u|umin u umax }.

Remark 2. For many systems, exact solution of the above nonlinear optimization problem may not be practical in real time. In such cases, the stability results obtained for the exact and global solution of the optimization problem can be extended to control laws obtained from only the feasible solution of the optimization problem (solutions that not necessarily minimize the objective function). Such controllers, generally referred as suboptimal model predictive controllers, are discussed in Scokaert et al. (1999). 4.1. Robust stability

In this section we employ the approximate low order statespace model identified in the previous section to design a receding horizon controller to regulate the surface roughness of the thin-film at a desired level during the deposition process. We assume that all the states are available for measurement. We recall that the identified dynamic system (assuming full state estimation) is represented as x(t + 1) = Ax(t) + N u(t) ⊗ x(t) + Bu(t),

Similar constraints may exist on the state, however, such constraints are not considered in the current formulation. Solution of the above optimization problem yields an optimal input sequence uM at each sampling instance. Only the first input vector in the sequence is actually implemented, such that M (x) = uM (0) provides the implicit control law. Subsequently, the prediction horizon is moved forward by one time-step, and the above problem is resolved using new process measurements.

Stability analysis for receding horizon control systems has received considerable attention in the past (Mayne et al., 2000). However, there are instances of receding horizon formulations based on nominal systems that have absolutely no robustness to disturbances or uncertainties (Grimm et al., 2005). The concept of robust stability, that is, maintaining the stability properties of the closed-loop system in the presence of disturbances has also been analyzed based on the nominal system (De Nicolao et al., 1996; Teel, 2004). An alternative approach to achieve robustness is to employ min–max feedback formulation, which considers all possible realizations of the disturbance in the open-loop optimal control problem (Scokaert and Mayne, 1998). For the following discussion, let S (x) denotes the set of solutions of the system of equation (33) starting at x with max{d}. The following definition of robust practical asymptotic stability is adopted from Teel (2004): Definition 4.1. Let S be compact and let R contain a neighborhood of S. For the system of equation (33) the set S is said to be semiglobally practically (in the parameter M) robustly asymptotically stable on int(R) if there exists  ∈ KL, for each compact set C and > 0 there exist M ∗ ∈ Z  0 , and for each M M ∗ there exist  > 0 (also depending on C and ) such that, for each M M ∗ , x ∈ C, and  ∈ S (x), we have | (t, x)|S (|x|S , t) + for all t ∈ Z  0 . Definition 4.2. Consider the system of equation (33). A function : Rn → R  0 is said to be detectable from a function l: Rn ×Rm → R  0 with respect to (¯W , W , W ) if W ∈ K∞ ,

W ∈ K∞ , ¯ W ∈ G and there exists a continuous function W : Rn → R  0 such that for all x ∈ Rn and u ∈ U, W (x) ¯ W ((x)), W (x(t + 1)) − W (x(t)) − W ((x)) + W (l(x, u)).

(35)

A. Varshney, A. Armaou / Chemical Engineering Science 63 (2008) 1246 – 1260

The following theorem states the sufficient conditions for the system of equation (33) to be practically robustly asymptotically stable:

1253

Adsorption

Gas-phase atom

Theorem 4.1 (Teel, 2004). Consider the system of Eq. (33) and the optimization problem of equation (34). Assume that the following hold:

Surface species Immobile species

1. The functions g and l are continuous. 2. For each integer M and x ∈ Rn , there exists a feasible ∗ (x) such that J (x, u∗ (x))=V (x). Moreover, for each uM M M M compact set Z ∈ Rn and integer M there exists  > 0 such ∗ (x) . that x ∈ Z implies uM 3. The set {x: (x) = 0} is compact, where : Rn → R  0 is a measure of the state. 4. For the system of equation (33), the function  is detectable from l with respect to some (¯W ∈ G, W ∈ K∞ , W ∈ K∞ ). 5. There exists ¯ ∈ K∞ such that Vi (x) ¯ ((x)) for all x and all i ∈ {1, . . . , ∞}.

Surface diffusion

Then for the receding horizon closed-loop the set {x: (x) = 0} is practically robustly asymptotically stable on Rn . Remark 3. A sufficient formula for M ∗ as given in Tuna et al. (2006) is M ∗ = 1 + L ln( (L − 1)),

(36)

where L is an upper bound for a sequence {Li } such that 1 Li L and Vi (x)Li (x), ∀i ∈ {1, 2, ...}. 0 is such that g(x(t +1))+l(x, u) (1+ )g(x). Eq. (36) is based on the assumptions that L and exist, (x)l(x, u), ∀u ∈ U, and there exist a control sequence u such that VM (x) = JM (x, u) for each horizon M. 5. Application to one-dimensional deposition process 5.1. Process description We consider an epitaxial thin-film growth process where adsorption of molecules from the gas-phase and diffusion of molecules on the film surface are the dominant phenomena responsible for the evolution of the microstructure of the thinfilm. The evolution of microstructure is modeled using the kMC algorithm on a one-dimensional regular lattice. During adsorption, it is assumed that the incoming molecule randomly attaches to any of the lattice sites with equal probability. Surface diffusion is modeled as hopping of surface molecules to adjacent lattice sites. We consider interaction among first-neighbors and the activation barrier for hoping of a molecule to be dependent upon the number of molecules adjacent to it. Desorption is neglected. The schematic of the deposition process is given in Fig. 1. The probability of the lattice being in a specific configuration  is given by the following master equation (Fichthorn and Weinberg, 1991; VanKampen, 1992): jP (, t)  W ( , )P ( , t) − W (,  )P (, t), = jt  

Fig. 1. A schematic of the deposition process showing adsorption and surface diffusion.

where P (, t) is the probability that the system is in state  at time t and W (,  ) is the probability per unit time of transition from configuration  to  . kMC provides a numerical solution of the above master equation through Monte Carlo sampling. The solution of the master equation is achieved computationally by executing an event chosen randomly among various possible events (adsorption and surface diffusion in the current case) based on the instantaneous event probabilities. The probability of an event being executed is proportional to its instantaneous rate. The rates of adsorption and diffusion are calculated as    kB T E + nE ra = j · NT , rd (n) = NT fn exp − , h kB T n (37) where ra and rd (n) are rates of adsorption and surface diffusion, respectively, j is the adsorption flux, fn is the fraction of sites having n nearest neighbors, and NT is the total number of lattice sites. h denotes the Planck’s constant, E is the energy barrier for surface diffusion, E is the interaction energy between two neighboring adsorbed molecules and n ∈ {0, 1, 2} denotes the number of nearest neighbors to a molecule. In this work the values of E and E were assumed to be 2.0 and 0.5 eV, respectively. Upon execution of an event, time is incremented by computing the interevent time as t =

ln  , ra + n rd (n)

(38)

where  is a random variable uniformly distributed between (0, 1).

A. Varshney, A. Armaou / Chemical Engineering Science 63 (2008) 1246 – 1260

14

10

4 2 0.1 0 -0.1 0

φ (1)

2.5

2

3

4

5

3

4

5

Original surface Reconstructed surface

1.5 0.5 0.1 0 -0.1

(39)

0

i

1

2 time, (s)

Fig. 2. Comparison of temporal evolution of (0) and (1) for random variation in adsorption flux and substrate temperature for reconstruction based on ACF and HDCF.

1300

Input

from random surfaces with identical ACF and HDCF under a range of adsorption fluxes and substrate temperatures. Surface roughness, which is a coarse variable, is widely used as a measure characterizing the microstructure during thinfilm deposition (Raimondeau and Vlachos, 2000; Lou and Christofides, 2003; Varshney and Armaou, 2005). An original surface was obtained from kMC simulations and random surfaces, which are statistically identical to the original surface with respect to ACF and HDCF, were constructed following a stochastic reconstruction algorithm presented in Varshney and Armaou (2006a). In Fig. 2, we present the comparison between the evolution of surface roughness for a number of random original and reconstructed surfaces, subject to random adsorption flux in the range 0–10 monolayers/s and random surface temperature in the range 900–1300 K. The agreement between the two indicate that dominant microstructural traits, relevant to the dynamics of coarse variable of interest (i.e., R), are adequately captured by the two correlation functions. During system identification, a data set comprising of 2048 snapshots obtained from kMC simulations on a lattice of size 100 × 100 was initially generated. During snapshot generation the inputs were randomly excited according to Gaussian random sequence with mean and variance equal to the mid-point and the half-width of their interval of variation, respectively. Fig. 3 shows the profile of substrate temperature used for snapshot generation. Snapshots were obtained as the profiles of lattice heights from which the corresponding correlation functions were subsequently extracted. Application of POD on the ensemble resulted in three dominant POD modes which accounted for more than 99.99% of the energy. The subspace identification

1

3.5

Error

For the above deposition process, it was found that the ACF and HDCF were sufficient to capture the dominant traits of the microstructure. Moreover, it was found that, (r) and (r) at r = 0, 1, 2 were sufficient for accurate reconstruction. Moreover, Eq. (5) with n = 4 satisfactorily captured the first-order statistics. Hence, a finite dimensional vector z ≡ [(0), (1), (2), (1), (2), (3), 3 , 4 ]T can be employed to approximately represent the coarse microscopic state. To elaborate on the above, we investigate the dynamic evolution of surface roughness, defined as (Raimondeau and Vlachos, 2000) 1  t t (|Xi+1 −Xit |+|Xi−1 − Xit |) = (1, t) 2NT

8 6

5.2. Identification of the approximate state-space model

R(t)=

Original surface Reconstructed surface

12 ρ (0)

Evolution of thin-film microstructure is accurately predicted through kMC simulations which are directly based on microscopic events provided these events and the corresponding rate parameters are accurately known. However, since kMC probabilistically samples the underlying master equation and does not furnish closed-form dynamic models, it is difficult to design feedback controllers based on them. Due to this limitation, we employ the methodology of the previous section to identify an approximate state-space model, which is subsequently employed to design a receding horizon controller.

Error

1254

1100 time Fig. 3. Profile of temperature used for data generation during system identification.

algorithm was employed, followed by solution of a NLP given by Eq. (32) to obtain the approximate bilinear state-space model as described in Section 3.4. The NLP was solved using SQP to obtain the matrices A, B and N which minimized the prediction error between the approximate model and the kMC data. A

A. Varshney, A. Armaou / Chemical Engineering Science 63 (2008) 1246 – 1260

Surface roughness

Surface roughness

0.75

Approximate Model kMC

0.9

1255

0.6

0.3

0

0.5

Open loop Closed loop

Set Point

0.25

0 0

25

50

0

time (s)

first-order McLaurin series expansion with p = 1 of Eq. (19) followed by a truncation of the infinite series of Eq. (22) with q = 1 provided a sufficiently accurate empirical model. In order to analyze the predictive capability of the identified model, in Fig. 4, we compare the evolution of surface roughness obtained from kMC simulations with the reduced model’s prediction. It is observed that the identified model approximates the true system evolution quite well and the normalized error between the kMC simulations and the model prediction is 4.2%. Note that these kMC simulations are independent of the simulations used for the generation of the model. 5.3. Closed-loop simulation results with receding horizon control

Fig. 5. Open-loop and closed-loop surface roughness profiles during deposition for growth rate j = 1 monolayers/s.

1240

W=1 W=0.25

1175

1100 0

The receding horizon controller described in the previous section was initially implemented to a thin-film growth process with growth rate W =1 monolayer/s. The control variable, i.e., the substrate temperature, was manipulated between 1100 and 1300 K. State measurements were assumed to be available every 1 s. In Eq. (34), the stage cost and terminal costs were chosen to be l(x, u) = x T Qx + uT Ru and g(x) = x T Sx, respectively, where Q > 0, S > 0, and R 0 are diagonal matrices. With the assumption that A is stable, taking (x) = (Q)x T x, where (Q) denotes the minimum eigenvalue of Q and setting ¯ {Li } = L with L = (Z)/ (Q) where Z is such that AT ZA − Z  −Q, S−Z 0, it can be shown that with appropriate choice of , all the assumptions for Eq. (36) are satisfied. Using Eq. (36) resulted in M ∗ = 1 for the current case. Motivated by the above result, the control horizon was set to 5 s. With this choice of the control horizon M, g(x), l(x, u), and U = [1100, 1300], we note that assumptions 1 and 2 of Theorem 1 are satisfied. With (x) = |x|2 , assumption 3 is also satisfied. Moreover, assumptions 4 can be satisfied by taking W (x) ≡ 0, ¯ W ≡ 0,

50

time (s)

Surface temperature (K)

Fig. 4. Expectation of the temporal evolution of surface roughness obtained from kMC simulations and low-order model for random variation in adsorption flux and substrate temperature.

25

25

50

time (s)

Fig. 6. Temporal profiles of the manipulated variable, the surface temperature, during deposition for growth rates j = 1 and j = 0.25 monolayers/s.

W = W = Id, where Id is the identity function defined as Id(x) = x. Hence, using Theorem 1, the closed-loop system is robustly asymptotically practically stable. The objective of the controller was to stabilize the surface roughness of the deposited film to a set-point which was chosen as Rs = 0.3. Fig. 5 presents the open-loop and closed-loop surface roughness profiles and Fig. 6 presents the manipulated input (surface temperature) profile for the closed-loop simulation. It can be seen that the controller successfully regulates the surface roughness of the film to the desired set-point while maintaining robustness with respect to noise inherent to kMC simulations. In contrast, the surface roughness is approximately 100% higher than the set-point for the open-loop operation.

1256

A. Varshney, A. Armaou / Chemical Engineering Science 63 (2008) 1246 – 1260

0.4

0.4

Surface Roughness

Open loop Closed loop

Surface roughness

Set Point

0.2

0.2

0 0 0 0

25

50

time (s) Fig. 7. Open-loop and closed-loop surface roughness profiles during deposition for growth rate j = 0.25 monolayers/s.

Fig. 9. Closed-loop surface roughness profile during deposition for step-disturbance in growth rate as shown in Fig. 8.

reject the disturbance and regulate the surface roughness close to the desired set-point.

1.2

6. Application to anisotropic two-dimensional epitaxial growth Substrate Temperature Growth Rate

0.6

1170

Growth Rate

Substrate Temperature

1240

50

25 time (s)

0

1100 0

25 time (s)

Fig. 8. Input disturbance in growth rate and closed-loop surface temperature profiles.

Subsequently, we analyzed the closed-loop performance by changing the growth rate to 0.25 monolayers/s. The set-point for the surface roughness was kept constant at Rs = 0.3, and, after analysis, the control horizon was again set to 5 s. The corresponding closed-loop surface temperature profile is shown in Fig. 6 and closed-loop and open-loop surface roughness profiles are shown in Fig. 7. As with the previous case, the closed-loop surface roughness was successfully regulated at the set-point. As expected, the corresponding closed-loop values of surface temperature are lower than the previous case owing to a lower value of adsorption flux. Finally, the performance capabilities of the controller were illustrated by introducing a pulse disturbance as shown in Fig. 8. The corresponding closed-loop control and surface roughness trajectories are shown in Figs. 8 and 9. Similarly to the previous cases, the controller was able to successfully

6.1. Process description To further illustrate the applicability of the proposed identification and control methodology, we consider the anisotropic growth of a simple cubic lattice using solid-on-solid model. Such models have been widely utilized to describe the growth of a variety of thin-films including Si(0 0 1) (Clarke et al., 1990), gallium arsenide (Ishii and Kawamura, 1999), indium arsenide, etc. The schematic of the deposition process is shown in Fig. 10. The microscopic processes included in the growth model are random deposition of atoms from the gas-phase and migration of the surface atoms. As before, desorption is neglected, consistent with the fact that surface residence time of molecules in metallic films are comparatively large. Adsorption occurs at each lattice site with equal probability. On the surface, the atoms are allowed to migrate with an Arrhenius-type site-dependent migration rate, rd =

kB T −E/kB T e , h

(40)

where the activation barrier E for an atom depends on the local configuration. In order to introduce anisotropy, we consider the activation barrier to of the following form: E = E 0 + n  E + n ⊥ E⊥ ,

(41)

where E0 is the surface term and E and E⊥ are the nearestneighbor bonds parallel and perpendicular, respectively, to the enhancement direction and n and n⊥ denote the number of such bonds, respectively. For instance, in case of Si(0 0 1) growth, E denotes the barrier in the direction parallel to dimer formation and in case of GaAs(0 0 1) growth (Shitara

A. Varshney, A. Armaou / Chemical Engineering Science 63 (2008) 1246 – 1260

Adsorption

1257

0.7

ρ (0)

Original Surface Reconstructed Surface 0.35

Migration

0 0

10

20

30

time Fig. 10. The schematic of the deposition process showing adsorption and anisotropic surface diffusion.

et al., 1992),  and ⊥ may denote the directions [1¯ 1 0] and [1 1 0], respectively. The parameters E0 , E , E⊥ were chosen to be those reported in Heyn (2001) with E0 = 1.3 eV and E + E⊥ = 1.2 eV, E /E⊥ = 5. 6.2. Identification of the approximate state-space model Since, the surface diffusion is greater in the direction perpendicular to the enhancement direction, the resulting thinfilm will be anisotropic. Hence, separate sets of correlation functions, each for directions parallel and perpendicular to the enhancement direction, are required to macroscopically characterize the deposition surface. To account for the anisotropy of the surface, we define the ACF in the direction parallel and perpendicular to the enhancement direction as follows:  N −r T  1  (r) = (Xi,j − X )(Xi+r,j − X ) , NT (NT − r) i=1

⊥ (r) =

1 NT (NT − r)

N −r T  j =1

j



(Xi,j − X )(Xi,j +r − X ) . i

(42) The definition of ACF and HDCF were modified to account for two-dimensional lattice. The corresponding  and ⊥ are defined analogously. In Fig. 11, we compare the evolution of R from a number of original surfaces and surfaces reconstructed based on  , ⊥ ,  and ⊥ . The agreement between the two is attributed to the fact that the anisotropy of the deposition surface was accurately captured by defining separate correlation functions, parallel and perpendicular to the enhancement direction.

Fig. 11. Comparison of temporal evolution of (0) for random variation in adsorption flux and substrate temperature for reconstruction based on  , ⊥ ,  , and ⊥ .

During system identification, a data set comprising of 2048 snapshots obtained by exciting the inputs using a Gaussian random sequence was initially generated. Application of POD on the ensemble resulted in three dominant POD modes which accounted for more than 99.99% of the energy. The subspace identification algorithm was employed, followed by solution of a NLP given by Eq. (32) to obtain the approximate bilinear state-space model of Section 3.4. The NLP was solved using SQP to obtain the matrices A, B and N which minimized the prediction error between the approximate model and the kMC data. Finally, a first-order McLaurin series expansion with p=1 in Eq. (19) followed by truncation of infinite series of Eq. (22) at q = 1 provided sufficiently a accurate empirical model. 6.3. Closed-loop simulations with receding horizon control The receding horizon controller, described in Section 4, was implemented to regulate the average surface roughness of the film at a desired set-point, which was chosen to be 0.2 monolayers. The growth rate was assumed to be W =1 monolayers/s. The control variable, i.e., substrate temperature, was manipulated between 550 and 750 K. As with the previous case, state measurements were assumed to be available every 1 s and the stage cost and terminal costs were chosen to be l(x, u) = x T Qx + uT Ru and g(x) = x T Sx, respectively, where Q > 0, S > 0, and R 0. Application of Eq. (36) resulted in M ∗ = 1 for the current case also. Motivated by the above result, the control horizon was set equal to 5 s. Similar to the previous case, the current closed-loop system can be shown to be robustly asymptotically practically stable. Fig. 12 presents the open-loop and closed-loop surface roughness profiles and Fig. 13 presents the manipulated input (surface temperature) profile for the closed-loop simulation. It can be seen that the controller was

1258

A. Varshney, A. Armaou / Chemical Engineering Science 63 (2008) 1246 – 1260 1.8

1.8

Open loop Closed loop

1.2

Surface Roughness

Surface Roughness

Open loop Closed loop

0.6

Set Point

1.2

0.6

Set Point

0

0 0

25

50

0

time (s)

25

50

time (s)

Fig. 12. Open-loop and closed-loop surface roughness profiles for growth rate of j = 2 monolayers/s.

Fig. 14. Closed-loop surface roughness profile for step-disturbance in growth rate shown in Fig. 15.

760 2.5

720

660

2 Growth rate

Surface Temperature

Surface Temperature

750

650 1.5

600

550 0

0

25

50

25 time (s)

1.0 50

time (s) Fig. 13. Temporal profiles of the manipulated variable, the surface temperature, during deposition for growth rate of j = 2 monolayers/s.

successful in regulating the surface roughness of the film to the desired set-point while the surface roughness was approximately 100% higher than the set-point for the open-loop operation. The performance of the controller was also analyzed by introducing a pulse disturbance as shown in Fig. 15. The corresponding closed-loop control and surface roughness trajectories are shown in Figs. 14 and 15. It can be seen from the figure that the controller successfully rejected the disturbance and regulated the surface roughness close to the desired set-point.

Fig. 15. Input disturbance in growth rate and closed-loop surface temperature profiles.

7. Conclusions The present work outlined a methodology for feedback control of microstructure during thin-film growth. The issue of non-availability of closed form dynamic models for the evolution of microstructure was addressed by deriving a low-order state-space model for a finite set of observables, identified from spatial correlation functions, which approximates the underlying kinetic Monte Carlo model. Initially, the data obtained from spatial correlation functions was projected onto a (reduced) subspace using an optimal basis computed through proper

A. Varshney, A. Armaou / Chemical Engineering Science 63 (2008) 1246 – 1260

orthogonal decomposition. Subsequently, the unknown dynamic model was transformed into an equivalent bilinear form using Carleman linearization. Finally, the bilinear statespace model was identified using a subspace based procedure followed by nonlinear least-squares minimization. The statespace model was subsequently employed to design a receding horizon controller that regulates the surface roughness of the thin-film at a specified set-point during the growth process by manipulating the substrate temperature. The proposed methodology was applied to two distinct deposition process. Closed-loop simulations at different growth rates and in the presence of a disturbances were performed to demonstrate the effectiveness of the controller. Acknowledgements Financial support from National Science Foundation, CAREER Award #CBET 06-44519, American Chemical Society, Research Initiation Award #PRF 44300-G9 and Pennsylvania State University, Dean’s fund is gratefully acknowledged. References Akiyama, Y., Imaishi, N., Shin, Y.S., Jung, S.C., 2002. Macro- and microscale simulation of growth rate and composition in MOCVD of Yttriastabilized zirconia. Journal of Crystal Growth 241, 352–362. Armaou, A., 2005. Output feedback control of dissipative distributed processes via microscopic simulations. Computers & Chemical Engineering 29, 771–782. Armaou, A., Christofides, P.D., 1999. Plasma-enhanced chemical vapor deposition: modeling and control. Chemical Engineering Science 54, 3305–3314. Armaou, A., Christofides, P.D., 2002. Dynamic optimization of dissipative PDE systems using nonlinear order reduction. Chemical Engineering Science 57, 5083–5114. Armaou, A., Baker, J., Christofides, P.D., 2001. Feedback control of plasma etching reactors for improved etching uniformity. Chemical Engineering Science 56, 257–265. Aroian, L.A., 1937. The type B Gram–Charlier series. The Annals of Mathematical Statistics 8 (4), 183–192. Badgwell, T.A., Breedijk, T., Bushman, S.G., Butler, S.W., Chatterjee, S., Edgar, T.F., Toprac, A.J., Trachtenberg, I., 1995. Modeling and control of microelectronics materials processing. Computers & Chemical Engineering 19, 1–41. Banks, H.T., Beeler, S.C., Kepler, G.M., Tran, H.T., 2002. Reduced order modeling and control of thin film growth in HPCVD reactor. SIAM Journal of Applied Mathematics 62, 1251–1280. Boas, R.P., 1949. Representations of probability distributions by Charlier series. The Annals of Mathematical Statistics 20 (3), 376–392. Christofides, P.D., 2001. Nonlinear and Robust Control of PDE Systems: Methods and Applications to Transport-Reaction Processes. Birkhäuser, Boston. Christofides, P.D., Daoutidis, P., 1997. Finite-dimensional control of parabolic PDE systems using approximate inertial manifolds. Journal of Mathematical Analysis & Applications 216, 398–420. Clarke, S., Wilby, M.R., Vvedensky, D.D., Kuwamura, T., Miki, K., Tokumoto, H., 1990. Step stability, domain coverage, and nonequilibrium kinetics in Si(0 0 1) molecular-beam epitaxy. Physical Review B 41, 10198–10201. De Nicolao, G., Magni, L., Scattolini, R., 1996. On the robustness of recedinghorizon control with terminal constraints. IEEE Transactions on Automatic Control 31, 451–453. Edgar, T.F., Butler, S.W., Campbell, W.J., Pfeiffer, C., Bode, C., Hwang, S.B., Balakrishnan, K.S., Hahn, J., 2000. Automatic control in microelectronics

1259

manufacturing practices, challenges, and possibilities. Automatica 36, 1567–1603. Emami-Naeini, A., Ebert, J.L., de Roover, D., Kosut, R.L., Dettori, M., Porter, L.L., Ghosal, S., 2003. Modeling and control of distributed thermal systems. IEEE Transactions on Control Systems Technology 11, 668–683. Favoreel, W., 1999. Subspace methods for identification and control of linear and bilinear systems. Ph.D. Thesis, Katholieke Universiteit Leuven, Belgium. Favoreel, W., Moor, B.D., Overschee, P.V., 1999. Subspace identification of bilinear systems subject to white inputs. IEEE Transactions on Automatic Control 44, 1157–1165. Favoreel, W., Moor, B.D., Overschee, P.V., 2000. Subspace state space system identification for industrial processes. Journal of Process Control 10, 149–155. Fenner, D.B., 2004. Fractal topography of surfaces exposed to gas-cluster ion beams and modeling simulations. Journal of Applied Physics 95, 5408–5418. Fichthorn, K.A., Weinberg, W.H., 1991. Theoretical foundations of dynamical Monte Carlo simulations. Journal of Chemical Physics 95, 1090–1096. Gallivan, M.A., Murray, R.M., 2004. Reduction and identification methods for Markovian control systems, with application to thin film deposition. International Journal of Robust & Nonlinear Control 14, 113–132. Gneiting, T., Schlather, M., 2004. Stochastic models that separate fractal dimension from the Hurst effect. SIAM Review 46, 269–282. Grimm, G., Messina, M.J., Tuna, S.E., Teel, A.R., 2005. Examples when nonlinear model predictive control is nonrobust. Automatica 40, 1729–1738. Henson, M.A., 1998. Nonlinear model predictive control: current status and future directions. Computers & Chemical Engineering 23, 187–202. Heyn, C., 2001. Anisotropic island growth during submonolayer epitaxy. Physical Review B 63, 033403. Hildebrandt, E.H., 1931. Systems of polynomials connected with Charlier expansions and the Pearson differential and difference equations. The Annals of Mathematical Statistics 2, 379–439. Holmes, P., Lumley, J.L., Berkooz, G., 1996. Turbulence, Coherent Structures, Dynamical Systems and Symmetry. Cambridge University Press, New York. Ishii, A., Kawamura, T., 1999. Monte carlo simulation of homoepitaxial growth on two-component compound semiconductor surfaces. Surface Science 436, 38–50. Kowalski, K., Steeb, W.H., 1991. Nonlinear Dynamical Systems and Carleman Linearization. World Scientific Publishing Co. Pte. Ltd., Singapore. Ljung, L., 1987. System Identification: Theory for the User. Prentice-Hall, Inc., Englewood Cliffs, NJ. Lou, Y., Christofides, P.D., 2003. Estimation and control of surface roughness in thin film growth using kinetic Monte-Carlo methods. Chemical Engineering Science 58, 3115–3129. Lou, Y., Christofides, P.D., 2005. Feedback control of surface roughness using stochastic PDEs. A.I.Ch.E. Journal 51, 345–352. Lou, Y., Christofides, P.D., 2006. Nonlinear feedback control of surface roughness using a stochastic PDE: design and application to a sputtering process. Industrial & Engineering Chemistry Research 45, 7177–7189. Mayne, D.Q., Rawlings, J.B., Rao, C.V., Scokaert, P.O.M., 2000. Constrained model predictive control: stability and optimality. Automatica 36, 789–814. Morari, M., Lee, J.H., 1999. Model predictive control: past, present and future. Computers and Chemical Engineering 23, 667–682. Ni, D., Christofides, P.D., 2005. Multivariable predictive control of thin film deposition using a stochastic PDE model. Industrial & Engineering Chemistry Research 44, 2416–2427. Park, H.M., Yoon, T.Y., Kim, O.Y., 1999. Optimal control of rapid thermal processing systems by empirical reduction of modes. Industrial & Engineering Chemistry Research 38, 3964–3975. Raimondeau, S., Vlachos, D.G., 2000. Low-dimensional approximations of multiscale epitaxial growth models for microstructure control of materials. Journal of Computational Physics 160, 564–576. Rathinam, M., Petzold, L.R., 2003. A new look at proper orthogonal decomposition. SIAM Journal on Numerical Analysis 41, 1893–1925.

1260

A. Varshney, A. Armaou / Chemical Engineering Science 63 (2008) 1246 – 1260

Rozman, M.G., Utz, M., 2002. Uniqueness of reconstruction of multiphase morphologies from two-point correlation functions. Physical Review Letters 89, 135501. Scokaert, P.O.M., Mayne, D.Q., 1998. Min–max feedback model predictive control for constrained linear systems. IEEE Transactions on Automatic Control 43, 1136–1142. Scokaert, P.O.M., Mayne, D.Q., Rawlings, J.B., 1999. Suboptimal model predictive control (feasibility implies stability). IEEE Transactions on Automatic Control 44, 648–654. Shitara, T., Vvedensky, D.D., Wilby, M.R., Zhang, J., Neave, J.H., Joyce, B.A., 1992. Misorientation dependence of epitaxial growth on vicinal GaAs(0 0 1). Physical Review B 46, 6825–6833. Shvartsman, S.Y., Kevrekidis, I.G., 1998. Nonlinear model reduction for control of distributed parameter systems: a computer assisted study. A.I.Ch.E. Journal 44, 1579–1595. Sobczyk, K., Kirkner, D.J., 2001. Stochastic Modeling of Microstructures. Birkhäuser, Basel. Svoronos, S., Papageorgiou, D., Tsiligiannis, C., 1994. Discretization of nonlinear control systems via the Carleman linearization. Chemical Engineering Science 49 (19), 3263–3267. Teel, A.R., 2004. Discrete time receding horizon optimal control: is the stability robust? Lecture Notes in Control and Information Science, 301, 3–27. Theodoropoulou, A., Adomaitis, R.A., Zafiriou, E., 1998. Model reduction for optimization of rapid thermal chemical vapor deposition systems. IEEE Transactions on Semiconductor Manufacturing 11, 85–98.

Torquato, S., 2001. Random Heterogeneous Materials: Microstructure and Macroscopic Properties. Springer, New York. Torquato, S., 2002. Statistical description of microstructures. Annual Review of Materials Science 32, 77–111. Tuna, S.E., Messina, M.J., Teel, A.R., 2006. Shorter horizons for model predictive control. In: Proceedings of the 2006 American Control Conference, Minneapolis, Minnesota, USA, pp. 863–868. VanKampen, N.G., 1992. Stochastic Processes in Physics and Chemistry. Elsevier Science Publishers B.V., Amsterdam. Varshney, A., Armaou, A., 2005. Multiscale optimization of thin-film growth using hybrid PDE/kMC process systems. Chemical Engineering Science 60, 6780–6794. Varshney, A., Armaou, A., 2006a. Identification of macroscopic variables for low-order modeling of thin-film growth. Industrial & Engineering Chemistry Research 45, 8290–8298. Varshney, A., Armaou, A., 2006b. Optimal operation of GaN thin film epitaxy employing control vector parameterization. A.I.Ch.E. Journal 52, 1378–1391. Verdult, V., Verhaegen, M., 2001. Identification of multivariable bilinear state space systems based on subspace techniques and separable least squares optimization. International Journal of Control 74, 1824–1836. Voigtlander, B., 2001. Fundamental processes in Si/Si and Ge/Si studied by scanning tunneling microscopy during growth. Surface Science 43, 127. Zapien, J.A., Messier, R., Collin, R.W., 2001. Ultraviolet-extended real-time spectroscopic ellipsometry for characterization of phase evolution in BN thin films. Applied Physics Letters 78, 1982.