Ind. Eng. Chem. Res. 2001, 40, 1661-1672
1661
Control of Wafer Temperature Uniformity in Rapid Thermal Processing Using an Optimal Iterative Learning Control Technique Kwang S. Lee,* Jinho Lee, Insik Chin, and Jinhoon Choi Department of Chemical Engineering, Sogang University, Shinsoo-1, Mapogu, Seoul 121-742, Korea
Jay H. Lee School of Chemical Engineering, Georgia Institute of Technology, 778 Atlantic Drive, Atlanta, Georgia 30332-0100
An iterative learning control (ILC) technique based on a quadratic optimal criterion is proposed for temperature uniformity control of a wafer under rapid thermal processing (RTP). The proposed technique is based on a time-varying linear state space model which approximates a nonlinear system along a reference trajectory. For effective use of feedback measurements, it employs a periodically time-varying Kalman filter which is based on an augmented model capable of describing the transition behaviors from one run to next, as well as those from one sample time to another. By retaining the relevant information from previous runs in the Kalman state, the control technique is able to make improvements in the control performance from one run to the next and eventually to converge to a minimum-achievable tracking error despite model error. A simple but effective method for identifying a time-varying linear state space model for the RTP system is also proposed. The performance of the proposed technique is evaluated through a simulation study involving the RTP model processing 8-in. silicon wafers. 1. Introduction After a decade of intensive research, rapid thermal processing (RTP) has emerged as an indispensable semiconductor manufacturing processes. The main advantage of RTP over the conventional furnace process is that various wafer manufacturing operations such as annealing, oxidation, nitridation, chemical vapor deposition, and cleaning can be conducted in a single system with reduced thermal budget, enhanced wafer granularity, and cluster compatibility.1 Despite the great promise, the widespread application of RTP has been slowed by the difficulty of attaining a required level of uniformity in the temperature distribution across the wafer surface. To overcome this barrier, research efforts have been made on both the design and control sides.2-4 As the standard size of the silicon wafer has grown from 6 to 8 in. and is now moving to 12 in., the temperature control requirement will continue to pose new challenges. From the viewpoint of control system design, the RTP system has several unique traits. First, it is a highly interacting, nonlinear, multivariable batch system requiring precise temperature control during the heating, processing, and cooling of a wafer, as well as during transitions from one phase to another. Second, whereas its fast dynamics, the sampling interval is limited to less than 0.2 s and can be as small as 0.02 s. Henceforth, any viable algorithm must be computationally efficient enough to allow this level of execution time. Certainly, on-line dynamic optimization based on a high-order fundamental model would not be feasible. Third, the temperature measurements are corrupted significantly by noise. The thermocouples, which are usually em* Author to whom all correspondence should be addressed. Phone: (82-2)705-8477. Fax: (82-2)3272-0319. E-mail: kslee@ ccs.sogang.ac.kr.
ployed for off-line calibration experiments, are susceptible to the electromagnetic noise from the high-power electric circuits as well the thermal noise from the tungsten-halogen lamps. The pyrometric measurement is sensitive to small changes in the surface state of the wafer. This means that filtering should be an integral part of the overall control algorithm. Finally, because a single RTP system is used for different wafer processing steps, e.g., annealing, oxidation, and nitridation, and because the associated reference temperature trajectory varies with the type of processing, it is necessary to use a nonlinear model or multiple linear models. Along with the control algorithm, an accompanying identification method must be developed to allow for easy, efficient model construction. In early studies of RTP system control, multiloop PID control was attempted but found to be unsatisfactory5 because of the strong interactions in the system. After that, the focus shifted to model-based multivariable control techniques, including model predictive control,6,7 linear quadratic Gaussian (LQG) control,8 LQG control with an optimally designed feedforward input signal,9 and internal model control.10,11 The performance of a model-based feedback control technique inevitably depends on the quality of the model. Given the fact that the characteristics of a RTP system change slowly over repeated operations by contamination of the chamber wall, the performance of the conventional model-based approaches might be limited unless the model is updated continuously. Integral control, a conventional method of fighting model errors through feedback, is not as effective here because the control problem is highly transient. Motivated by the above considerations, we propose an iterative learning control (ILC) technique specifically suited for temperature uniformity control in RTP systems. The ILC12-14 is a special control technique for
10.1021/ie0005553 CCC: $20.00 © 2001 American Chemical Society Published on Web 03/13/2001
1662
Ind. Eng. Chem. Res., Vol. 40, No. 7, 2001
Figure 1. Schematic diagram of the RTP system.
batch or repetitive processes that updates the current batch operation by feeding back the previous batch information. Indeed, wafer temperature control based on the exploitation of the repetitive nature of RTP operation has been attempted by other researchers. Chen et al.15 considered a so-called PI learning rule, which does not rely on a process model, and applied it to a single-input single-output (SISO) RTP model. Zafiriou et al.16,17 devised a nonlinear model-based runto-run scheme, which is a direct implementation of a numerical algorithm of optimal control to a real system. They demonstrated the performance with a SISO RTP model. The above approaches are rather preliminary in that they conduct only run-to-run improvement (lacking real-time feedback control) in a deterministic framework. Considering the above general background, the proposed technique is based on a time-varying linear state space model identified along the reference trajectory. A unique aspect of the technique is that it combines runto-run control with real-time feedback control so that the two work in an optimal and complementary manner. The integral action we embed in the run-to-run part of the control algorithm enables the controller to steer the system to a minimum-achievable tracking error despite model errors and repetitive disturbances. The real-time feedback control part fights random disturbances as they occur and also accelerates the convergence to an optimal trajectory. The deterministic formulation of the proposed control was originally put forward by Amann et al.,18 but the result is extended to a stochastic formulation in this paper for the more effective handling of measurement noises and disturbances that vary randomly from run to run. In a sense, the proposed technique can be considered as an extension of the linear model-based version of Zafiriou et al.’s run-torun scheme16,17 by the addition of real-time optimal control in a stochastic framework. In addition to the control algorithm, a simple method for identifying a time-varying linear state space model for a RTP system is also presented. The proposed control technique is evaluated in a numerical simulation study involving an 8-in. wafer RTP system. 2. Description of the RTP System The RTP system considered in this research is of the same type as studied by Norman.19 The system has been scaled up for 8-in. wafer processing. Figure 1 shows a schematic diagram of the system. It has three concentric tungsten-halogen lamp arrays with maximum powers of 2, 15, and 50 kW. The silicon wafer is assumed to be a cylinder with radius R ) 101 mm and thickness Z )
0.68 mm. It is placed on a rotating support to avoid an azimuthal temperature distribution. The wafer temperatures are assumed to be measured at four radial positions (r ) 0, R/3, 2R/3, R) on the wafer’s bottom surface. The reference temperature trajectory was designed to linearly ramp up at a rate of 60 K/s from 673 K, to hold at 1273 K for up to 30 s, and to linearly cool to 873 K at a rate of 26.7 K/s for 15 s. It is assumed that temperature uniformity is important up to 873 K in the cooling phase. The cooling rate was chosen so that all of three lamp arrays are not extinguished during cooling. The corners of the trajectory were smoothed to avoid abrupt changes in input powers. This trajectory was applied equally to the four measured outputs. Numerical RTP System. The heat conduction within the wafer was modeled by the following parabolic partial differential equation:
F
∂cp(T)T 1 ∂ ∂T ∂ ∂T ) k(T)r + k(T) ∂t r ∂r ∂r ∂z ∂z
[
]
[
]
(1)
The temperature dependencies of the heat capacity cp and thermal conductivity k of the silicon wafer were represented as quadratic functions of T. At the boundaries, the heat flux by conduction is balanced by the radiational heat flux plus the convective loss. The resulting boundary conditions are
at r ) 0
∂T )0 ∂r
at r ) R k(T)
∂T )q j rad e + he(Tgas - T) ∂r
at z ) 0 k(T)
∂T ) -q j rad b - hb(Tgas - T) ∂z (2)
∂T )q j rad + ht(Tgas - T) at z ) Z k(T) t ∂r In the above, subscripts e, b, and t denote the edge, bottom, and top of the wafer, respectively. q j rad represents the heat flux by radiation. The radiant heat input is again decomposed into radiant absorption and emissive loss. For numerical simulations, the silicon wafer was divided into 40 × 3 cylindrical annulus cells such that all of the annuli have a same area. The partial derivatives along the spatial coordinates are approximated using the centered differences. Denoting the wafer temperatures at the internal and boundary nodes (arranged in vector form) as Tin and Tbd, resepctively, the resulting discretized version of eq 1 is given by
dTin(t) ) f[Tin(t),Tbd(t)] dt
(3)
In the above, Tbd is obtained from the boundary equations. To discretize eq 2, the view factors between the wafer surface and other radiant surfaces need to be estimated first. For this, the chamber ceiling and floor are divided into 18 and 36 equal-width concentric annuli, and the chamber wall is divided into 16 equalheight cylinders. The lamp arrays are assumed to be located at the 2nd, 6th, and 15th concentric rings of the ceiling. By taking the backward difference on the derivatives in eq 2 and substituting the complex radiant heat transfer formula into the resulting equation,19 eq
Ind. Eng. Chem. Res., Vol. 40, No. 7, 2001 1663
2 can be rearranged to
0 ) Tbd(t) + A1T4bd(t) + A2p(t - d) + A3Twall + A4Tgas + A5Tin(t) (4) where T4bd is a vector of which the entries are the fourth-power terms of the boundary temperatures. p ) [p1 p2 p3]T represents the powers (in kilowatts) from each lamp array zone. We introduced 0.2 s of delay (d ) 0.2) intto each input to take account of the tardy activation of the tungsten-halogen lamps. As shown above, the RTP model is approximated by a set of differential-algebraic equations comprising eqs 3 and 4. The MATLAB source code for the two-dimensional RTP model is available at http:// procd.sogang.ac.kr/. As can be imagined from Figure 1, the RTP system is highly interacting, as the power from a lamp array has simultaneous effects on the wafer temperature at all positions. Also, as the temperature rises, the net heat exchange between the lamps and the wafer is reduced. This effect, together with the temperature dependent thermal conductivity and heat capacity of the silicon wafer, results in severe nonlinearity. From preliminary open-loop simulations, the sampling time for control was chosen as 0.1 s, which translates to a total of 450 samples in each run. 3. Model Identification The proposed control technique requires a timevarying linear state space model identified along a reference trajectory. Techniques for identifying timevarying linear models are not well-established. In this research, we employ a simple technique based on interpolating time-invariant state space models identified at two or more operating points. This technique is considered to be useful for RTP systems in that multiple time-varying models for different reference trajectories can easily be obtained in a constructive way. The following describes the overall procedure. Step 1. The entire range of operating temperatures is decomposed into a finite number of subranges. At a representative temperature of each subrange, an identification experiment is conducted, and a linear state space model is found using a subspace identification method, e.g., N4SID.20 Step 2. The state space models for different temperature ranges are transformed to a common state basis with respect to the input-output coordinates. Step 3. A time-varying linear state space model is constructed by interpolating among the transformed models according to the time versus temperature relationship of the reference trajectory. N4SID employed in step 1 provides the following stochastic state space model in a balanced form
x(t+1) ) Ax(t) + Bu(t) + Kv(t) y(t) ) Cx(t) + v(t) with Rv } E{v(t)v(t)T}
(5)
where u ∈ R3 and y ∈ R4 represent the lamp power inputs p (watts) and the temperature measurements T (K), respectively, and x(t) and v(t) denote the Kalman state and the innovation, respectively. Because of the balanced realization, the state bases for different temperatures are defined differently (with respect to the inputs and outputs) at each temperature. The similarity
transformation in step 2 is needed for this reason. Interpolation in step 3 yields the following time-varying state space model, which is valid along a certain reference trajectory:
x(t+1) ) A(t)x(t) + B(t)u(t) + K(t)v(t) y(t) ) C(t)x(t) + v(t) with Rv(t) } E{v(t)v(t)T}
(6)
The implementation of the proposed technique requires a more detailed explanation, especially about the basis adjustment and model interpolation. Because the state basis is determined from subspace identification, we begin with a brief review on N4SID. Review of N4SID and Basis Adjustment. N4SID is a subspace identification technique that allows us to find a minimal, balanced state space model using input-output data obtained from an open-loop experiment. Only the key steps leading up to the definition of state coordinates will be described here because the main purpose of the discussion is to explain how the state bases must be adjusted. Let t be the present time. With respect to t, we define the past input and output vectors and the future input and output vectors as follows:
up(t) }
[ ] [ ] [ ] [
u(t-n j) y(t-n j) , yp(t) } , l l u(t-1) y(t-1) u(t) y(t) uf(t) } , yf(t) } l l u(t+n j -1) y(t + n j -1)
]
(7)
Here, n j is assumed to be chosen larger than the state dimension n. In an optimal linear predictor, the future output is represented by a linear combination of the present Kalman state and the future inputs. Also, according to the Kalman filter equation, the present Kalman state is represented by a linear combination of the past inputs and the past outputs. This leads to an equation of the form
yf(t) ) L1up(t) + L2yp(t) + L3uf(t) + ef(t) ) Ox(t) + L3uf(t) + ef(t)
(8)
where ef represents the residual term (containing the future innovations) and O is the extended observability matrix defined by
[ ]
C l O} CAnj -1
(9)
In the case when the data are obtained from an openloop experiment, ef is uncorrelated with the other regressors, and hence, consistent estimates of L1, L2, and L3 can be obtained using the least-squares method. For the least-squares estimation, we stack the data vectors in eq 8 at different time points columnwise and write eq 8 in terms of the stacked data matrices. Assuming that we have s data points for the vectors, we can define
Yp } [yp(1) ... yp(s)], Up } [up(1) ... up(s)], X } [x(1) ... x(s)], etc.
(10)
where s . p. Then, eq 8 becomes
Yf ) L1Up + L2Yp + L3Uf + Ef ) OX + L3Uf + Ef (11)
1664
Ind. Eng. Chem. Res., Vol. 40, No. 7, 2001
Consequently
[ ]
U OX ) [L1L2] Y p } LP p
(12)
From the facts that the choice of the state can be arbitrary (within similarity transformations) and the only requirement is that the range spaces of OX and LP be same, O and X can be determined from the singular-value decomposition of LP.
LP ) [U1U2]
[ ][ ]
T S 0 V1 ) U1SV1T ) OX 0 0 V2T
(13)
It is obvious that the minimal system order n is equal to the rank of LP and, hence, can be chosen as the number of diagonal elements in S. Also, the extended observability matrix O and the state X can be determined as
O ) U1S1/2,
X ) S1/2V1T
(14)
In practice, because of the effects of various disturbances and the nonlinearity of the process, the singular values are nonvanishing, and the system order is determined as the number of larger singular values. During this process, insignificant system modes (in terms of the future output prediction) are automatically cut off, and ill-conditioning of the identified model is prevented. In eq 13, LPs obtained with data sets for different nominal temperatures share a common basis because these matrices are defined in terms of the physical inputs and outputs. On the other hand, the states X could be defined in terms of different bases for different nominal temperatures. Let the subscripts i and j denote two temperatures, and let us assume that the state dimensions for the two nominal temperatures are same. For the respective temperatures, eq 13 can be written as
OiXi ) LiPi,
OjXj ) LjPj
(15)
From the fact that both LiPi and LjPj have the same basis, the similarity transform matrix Vij between xi and xj is given as
xi ) Vijxj ) Oi+Ojxj f Vij ) Oi+Oj ) Si-1/2UiTUjSj1/2 (16) Obviously, Vij-1 ) Vji. Let (Ai, Bi, Ci, Ki) be the system matrices identified for the ith nominal temperature. Also, let the subscript r denote some reference temperature. Then the transformations needed to express (Ai, Bi, Ci, Ki) in terms of the basis used for the reference temperature case are given by
h i, C h i, K h i) ) (Vir-1AiVir, Vir-1Bi, CiVir, Vir-1Ki) (A h i, B (17) Model Interpolation. The model matrices in eq 17 are identified for only a limited number of nominal temperatures. To construct a time-varying model from these models, we generate additional models at various intermediate temperatures encountered during the operation. For this, we use linear interpolation. Here, it is logical to perform the interpolation in terms of the
dynamic gains. The interpolation of each system matrix done independently might produce a meaningless model. Let T1 and T2 be two adjacent temperatures, and let h i, C h i, K h i), i ) 1, 2, be the associated basis-adjusted (A h i, B system matrices. Let TR be a temperature between T1 and T2 such that
TR ) (1 - R)T1 + RT2, 0 e R e 1
(18)
Note that the frequency gain from the input to the state h )-1B h , and that of the dual is represented by (ejωI - A h )-1]T. Then, the system matrices system is [C h (ejωI - A at TR that linearly interpolate these two gains can be obtained by solving the following two relationships:
h R)-1[B hR K h R] ) (1 - R)(ejωI - A h 1)-1[B h 1K h 1] (ejωI - A h 2)-1[B h2 K h 2] + R(ejωI - A h R)-1 ) (1 - R)C h 1(ejωI - A h 1)-1 C h R(ejωI - A
(19)
+ RC h 2(ejωI - A h 2)-1 h R, The above equations are overdetermined., i.e., (A h R, B h R) that satisfy the above for all ω ∈[0, π] do not C h R, K exist. Instead, one might require the above to hold at selected finite frequencies and determine the interpolated matrices in the least-squares sense. By pre- and post-multiplying each equation by ejωI A h R, eq 19 can be rearranged to
h R] ) [R1(ω;R) + jI1(ω;R)] + A h R[R2(ω) + [B hR K jI2(ω;R)] (20) C h R ) [R3(ω;R) + jI3(ω;R)] + [R4(ω;R) + jI4(ω;R)]A hR where both the Ri and the Ii denote real matrices. The above leads to
h R] ) R1(ω;R) + A h RR2(ω;R) [B hR K hR C h R ) R3(ω;R) + R4(ω;R)A 0 ) I1(ω;R) + A h RI2(ω;R)
(21)
0 ) I3(ω;R) + I4(ω;R)A hR First, A h R can be determined by the least-squares method using the third and fourth relationships at various h R, K h R) are frequencies. Then, the estimates of (B h R, C directly given by the first and second relationships. In the above, C h R can be obtained by directly interpolating C h as C h defines the state-to-output gain. Nonetheless, the reason to interpolate C h (qI - A h )-1 using the h. dual system concept is for symmetry with (qI - A h )-1B -1 h represents the map of the past Just as (qI - A h) B inputs to the present state, C h (qI - A h )-1 can be interpreted as the map of the present state to the future outputs. 4. Optimal Iterative Learning Control Stochastic System Model for ILC Design. To design a controller that utilizes the information from the previous run, it is necessary to transform the underlying model in eq 6 to contain information about the disturbances in the previous runs that are related to those in the future runs. The first step in this process
Ind. Eng. Chem. Res., Vol. 40, No. 7, 2001 1665
is to introduce appropriate assumptions about the nature of the stochastic disturbances and then to construct a run-to-run transition model accordingly. Throughout this section, the subscript k represents the run index. In eq 6, v(t) is an iid (independent and identically distributed) sequence in time, but vk(t) for different run indices can show significant correlation. In fact, in terms of k, the error term vk(t) can exhibit “persistent” or “drifting” behavior in addition to random fluctuations. Such behavior can be reasonably modeled by the equation
vk(t) ) vj k(t) + vˆ k(t),
vj k - vj k-1(t) ) nk(t) (22)
where vˆ k(t) and nk(t) are independent random sequences with respect to both the k and t indices. Roughly speaking, vj k(t) represents the part of the error that will persist through all future runs, while vˆ k(t) is the part that will disappear. Now, we decompose eq 6 into two parts, one that is driven by uk(t) and vj k(t) and the other driven by vˆ k(t).
xjk(t+1) ) A(t)xjk(t) + B(t)uk(t) + K(t)vj k(t) yjk(t) ) C(t)xjk(t) + vj k(t) xˆ k(t+1) ) A(t)xˆ k(t) + K(t)vˆ k(t)
(23)
(24)
yˆ k(t) ) C(t)xˆ k(t) + vˆ k(t) Of course, yk(t) ) yjk(t) + yˆ k(t). Also, if we denote the reference trajectory for yk(t) by yd(t), the error signal becomes ek(t) ) yk(t) - yd(t). To put eq 23 into the standard form in which the external noise is an independent sequence in terms of k as well as t, we first take the difference between the equations for two consecutive runs to obtain
∆xjk(t+1) ) A(t)∆xjk(t) + B(t)∆uk(t) + K(t)nk(t) (25) yjk(t) ) C(t)∆xjk(t) + yjk-1(t) + nk(t) where ∆xjk } xjk - xjk-1 and ∆uk } uk - uk-1. Note that the output equation can be equivalently written as ejk(t) ) C(t)∆xjk(t) + ejk-1(t) + nk(t), where ejk(t) } yjk(t) - yd(t). Also, from the output eqs 23 and 24, it can be seen that
ek(t) ) C(t)xˆ k(t) + C(t)∆xjk(t) + ejk-1(t) + vˆ k(t) + nk(t) (26) To include ejk-1(t) in the state, we first define e jk } ... ejk(N)T]T. Then, from eq 25, we have
[ ][ ]
C(0)∆xjk(0) nk(0) C(1)∆xjk(1) n (1) j k-1 + + k e jk ) e l l C(N)∆xjk(N) nk(N)
[ejk(0)T
(27)
augmented state space equation:
[ ] [ ][ ] [ ] [ ] [ ] xˆ k(t+1) A(t) 0 0 xˆ k(t) ∆xjk(t+1) ) 0 A(t) 0 ∆xjk(t) + e j k(t+1) j k(t) 0 C(t) I e
K(t)vˆ k(t) 0 K(t)n B(t) ∆uk(t) + k(t) T 0 H (t)nk(t)
xˆ k(t) ˆ (t) C(t) H(t) ] ∆xˆ k(t) + vˆ k(t) + nk(t) ek(t) ) [C e j k(t)
where C(t) } [0 ... CT(t) ... 0]T and H(t) is a matrix that extracts the (t + 1)st element from e j k(t). For transition from one run to next, we can use
[ ][
assuming ∆xjk(i) ) 0 and n(i) ) 0 for i ) t, ..., N (28)
According to the above definition, the (t + 1)st element in e j k(t), which corresponds to the sample time t, is equal j k(t) in the state while to ejk-1(t). Now, by including e combining eqs 24 and 25, we obtain the following
][
][ ] [ ]
xˆ k+1(0) xˆ k(N+1) 0 0 0 xˆ k(N) ∆xjk+1(0) } ∆xjk(N+1) ) 0 0 0 ∆xjk(N) + 0 C(N) I e e j k+1(0) e j k(N+1) j k(N) ˆ k ˆ k (30) HT(N)nk(N)
Notice that ˆ k and jk represent the variations in the initial state from run to run, again differentiated as the randomly fluctuating type and the persistent type, respectively. The above state space equation is denoted as
h (t)zjk(t) + Γ h (t)∆uk(t) + ζhk(t), zjk(t+1) ) Φ h zjk-1(N) + ξhk zk(0) ) Ψ ek(t) ) Σ h (t)zjk(t) + η j k(t) t ) 0, ..., N
(31)
It can be viewed as a periodically time-varying state space system of period N + 1 since we defined zjk-1(N+1) to be zjk(0). BLQG Formulation. The input change ∆uk(t) for the kth run is determined according to the following linear quadratic Gaussian (LQG) objective: N-1
minJk ) E{ekT(N)Mek(N) +
∆uk(‚)
ekT(t)Qek(t) + ∑ t)0 ∆ukT(t)R∆uk(t)} (32)
subject to eq 31. The above is a standard LQG problem, but we named the resulting technique batch LQG (BLQG) control in accordance with the original motivation of the problem. The optimal solution is standard21 and given by a combination of the Kalman filter and the linear quadratic regulator (LQR)
h (t)zjk(t|t) ∆uk(t) ) -L
Now, we define
jk e j k(t) } e
(29)
(33)
where
L h (t) ) [Γ h (t)TP h (t+1)Γ h (t) + R]-1Γ h (t)TP h (t+1)Φ h (t) P h (t) ) Φ h (t)TP h (t+1)[Φ h (t) - Γ h (t)L h (t)] + Σ h (t)TQΣ h (t), h (N) P h (N) ) Σ h (N)TMΣ
(34)
1666
Ind. Eng. Chem. Res., Vol. 40, No. 7, 2001
In the above, zjk(t|t) denotes the state estimate that is obtained using the Kalman filter equation22 applied to eq 31
zjk(t+1|t) ) Φ h (t)zjk(t|t) + Γ h (t)∆uk(t) h k(t)[ek(t) - Σ h (t)zk(t|t-1)], zjk(t|t) ) zjk(t|t-1) + K zjk(0|-1) ) Ψ h zjk-1(N|N)
(35)
yk(t) - yjk-1(t|t) ) Σ(t)zk(t) + ηk(t)
Obviously, eq 39 has a state dimension considerably lower than that of eq 31. The Kalman filter equation is given by
h k(t)Σ h (t)T + Rζhηh ][Σ h (t)S h k(t)Σ h (t)T + Rηh ]-1 K h k(t) ) [Φ(t)S
˜ (t)[yk(t) - yk-1(t|t) zk(t|t) ) zk(t|t-1) + K Σ(t)zk(t|t-1)]
S h k(t+1) ) Φ h (t)S h k(t)Φ h (t)T +
yjk(t|t) ) C(t)∆xjk(t|t) + yjk-1(t|t) T
h k(t)[Σ h (t)S h k(t)Σ h (t) + Rηh ]K h k(t) Rζh - K
T
(36)
S h k(0) ) Ψ hS h k-1(N)Ψ h T + Rξh
5. A Suboptimal BLQG Formulation To enhance the computational efficiency while retaining the basic structure of the underlying model, a suboptimal algorithm is formulated. The key to the suboptimal formulation is to approximate eq 25 as
∆xjk(t+1) ) A(t)∆xjk(t) + B(t)∆uk(t) + K(t)nk(t) (37) yjk(t) ) C(t)∆xjk(t) + yjk-1(t|t) + nk(t) Here, we have replaced yjk-1(t) with yjk-1(t|t), which represents an estimate based on information available at time t of the (k - 1)st run. By combining eqs 24 and 37, we obtain
][ ] [ ]
xˆ k(t+1) xˆ k(t) A(t) 0 0 ) + ∆uk(t) + ∆xjk(t+1) B(t) 0 A(t) ∆xjk(t) K(t)vˆ k(t) K(t)nk(t) yk(t) ) [C(t) C(t) ]
[ ]
with
[
S(t+1) ) Φ(t)S(t)Φ(t)T + ˜ (t)[Σ(t)S(t)Σ(t)T + Rη]K ˜ (t)T Rζ - K
]
(41)
S(0) ) E{[zk(0) - zk(0|-1)][zk(0) - zk(0|-1)]T} where
Rζη )
[
]
K(t) 0 R and 0 K(t) η Rζ )
[
][
K(t) 0 K(t)T 0 Rη 0 K(t) 0 K(t)T
]
The Kalman filter can be initiated at each run with zk(0|-1) ) [00]. The Kalman gains spacings are runinvariant and can be determined off-line based on the covariance matrices of xˆ k(0), ∆xjk(0), nk(t), and vˆ k(t). The information feedback from one run to another still occurs because of the feedforward term yjk-1(t|t). The input signal is calculated by solving the same objective as in 32. By taking yjk(t) - yjk-1(t|t) as the output of the state space model, the output error can be rearranged as ek(t) ) yk(t) - yd(t) ) [yk(t) - yjk-1(t|t)] - [yd(t) - yjk-1(t|t)]. Consequently, the LQG problem becomes a tracking problem for the output to follow yjd(t) - yjk-1(t|t). The solution is standard21 and takes the form of
∆uk(t) ) -Lfb(t)zk(t|t) + Lff(t)bk(t+1)
(42)
where
Lfb(t) ) [Γ(t)TP(t+1)Γ(t) + R]-1Γ(t)TP(t+1)Φ(t) Lff(t) ) [Γ(t)TP(t+1)Γ(t) + R]-1(t)T P(t) ) Φ(t)TP(t+1)[Φ(t) - Γ(t)Lfb(t)] + Σ(t)TQΣ(t), P(N) ) Σ(N)TMΣ(N)
xˆ k(t) + yjk-1(t|t) + vˆ k(t) + nk(t) ∆xjk(t) (38)
which is written in a simplified form as
(40)
K ˜ (t) ) [Φ(t)S(t)Σ(t)T + Rζη][Σ(t)S(t)Σ(t)T + Rη]-1
In the above, Rηj } E{η j k(t)η j k(t)T}, Rζhηj } E{ζhk(t)η j k(t)T}, so forth. Note that, because we made the model periodically time-varying, there is flow of information from one run to another during the transition. Specifically, zjk(N|N) is transformed to zjk(N+1|N), which serves as the initial state for the Kalman filter in the next run [zjk+1(0|-1)]. Because the model is periodically time-varying, the Kalman filter gains eventually become periodic, and hence, the steady-state periodic solution can be used. Because the LQR gains are also periodic, one obtains the same time-varying controller for every run, which can conveniently be computed off-line and implemented before the operation begins. Because the model state contains e j k, which collects the errors during all of the sample points during the kth run, if there are a lot of samples in a single run, the Kalman filter design and the LQR design can involve Riccati equations of very large dimensions. This motivates the development of a suboptimal algorithm, which is presented next.
][
(39)
zk(t+1|t) ) Φ(t)zk(t|t) + Γ(t)∆uk(t)
where
[
zk(t+1) ) Φ(t)zk(t) + Γ(t)∆uk(t) + ζk(t)
T
(43)
T
bk(t) ) [Φ(t) - Γ(t)Lfb(t)] bk(t+1) + Σ(t) Q[yd(t) yjk-1(t|t)], bk(N) ) Σ(N)TM[r(N) - yjk-1(N|N)] In the above, Lfb(t) and Lff(t) are run-invariant and can be estimated off-line; bk(t) is calculated after the (k 1st run. In summary, the only equations that need to be computed in real-time are eqs 40 and 42. The
Ind. Eng. Chem. Res., Vol. 40, No. 7, 2001 1667 Table 1. Parameters of the RTP Model parameter
value/formula
description
F cp (J kg-1 K-1) k (W m-2 K-1) ht, hb (W m-2 K-1) he (W m-2 K-1) Si ch σ (W m-2 K-4) Tgas (°K) Tch (°K) R (mm) Z (mm) r1, r2, r3 (mm)
2330 748.38 + 0.1678T 184.74 - 0.2548T + 9.988 × 10-5T2 7.1 + 4.3(r/R)4 11.4 0.6 0.5 5.68 × 10-8 300 300 100.8 0.675 38.3, 89.4, 204.4
density of Si specific heat of Si, valid over 600-1400 K thermal conductivity of Si, valid over 600-1400 K heat transfer coefficients at z ) 0 and Z heat transfer coefficient at r ) R emissivity of the Si wafer emissivity of the chamber wall Stefan-Boltzmann constant ambient gas temperature chamber wall temperature wafer radius wafer thickness lamp array positions
(Kg/m3)
suboptimal BLQG technique was employed for the numerical study of RTP control. By scrutinizing the input calculations, one can see that uk(t) from suboptimal, as well as optimal, BLQG control is determined as a summation of the output error over the run index. This run-wise integral action enables us to attain the minimum-achievable tracking error as k f ∞ (just as the integral action over time removes steady-state offsets) such that N-1
Jk f min E{eT(N)Me(N) +
eT(t)Qe(t)} ∑ t)0
as k f ∞
despite model uncertainty (within certain limits). 6. Results and Discussion Simulation Conditions. The major emphasis of the simulation study was placed on investigating the robust convergence behavior of the proposed control technique. For this, the (suboptimal) BLQG controller was designed for four different linear models, which will be referred to as model I through model IV. Model I is a timevarying model constructed by interpolating two timeinvariant state space models identified at 773 and 1273 K. Model II is a more refined model that was constructed with seven time-invariant state space models identified at every 100 K from 673 to 1273 K. Interpolation was conducted at every 10 K for both models. Models III and IV are time-invariant state space models identified at 773 and 1273 K, respectively. All of the time-invariant state space models were identified using N4SID. A state dimension of four was found to be appropriate for all temperatures where the identification was conducted. The RTP system was assumed to be operated under two different cases. The first is the nominal case where all of the parameters of the RTP system have nominal values (see Table 1), whereas the second one is a perturbed (or contaminated) case where the view factors in the RTP system are changed randomly by -5 to 0% from their nominal values. To avoid complexity in the presentation, the results for the perturbed case are shown only for model I. In all cases, a zero-mean Gaussian noise of variance 1 K2 was added independently to each of the temperature measurements. For comparison, the optimal BLQG technique was applied only to model I for both cases. The proposed BLQG technique has seven tuning factors: the weighting matrices R, Q, and M in the quadratic objective and the covariance matrices of xjk(0), ∆xjk(0), n(t), and vˆ (t). Of these matrices, R was chosen as diag[1.3, 0.03, 0.012] × 10-3 to reflect the different
scales of the input lamp powers; Q was given as diag[1, 1, 1, 1.2] after some trials to obtain a uniform temperature distribution; and M was set to zero as there is no motivation to emphasize the end-point control. The covariance matrices were chosen as cov{xˆ k(0)} ) cov{∆xjk(0)} ) 5I and cov{n(t)} ) cov{vˆ (t)} ) 0.5I. The larger value for the first two terms reflects the fact that, in a real system, the initial states (e.g., the temperature and the derivatives at the initial time) can deviate significantly from batch to batch. The covariance of vˆ k(t) was chosen to be a rough approximation of the covariance of the measurement noise, and the covariance of nk(t) was introduced to reflect the model error effects but determined without any quantitative ground. To initialize the BLQG controller for the nominal case, multiloop PID control was tried at the first run. The BLQG controller was applied subsequently until the input profiles converge. From a convergent run in the nominal case, the perturbed case simulation was initiated. For performance comparison, the conventional feedback LQG controller with no iterative learning was also designed and applied. The LQG controller was designed on the basis of the following state space model derived from eq 6 to include the integral action:
[
][
][ ] ] [ [ ]
δx(t+1) A(t) 0 δx(t) ) + z(t+1) C(t+1)A(t) 1 z(t) B(t) K(t) δu(t) + δv(t) C(t+1)B(t) C(t+1)K(t)
[
y(t) ) [0 I ]
δx(t) + v(t) z(t)
]
(45)
where δ denotes the time-difference operator such that δx(t) } x(t) - x(t-1). In the above derivation, we introduced the approximation A(t)x(t) - A(t-1)x(t-1) ≈ A(t)δx(t) and similarly for other system matrices. The LQG controller was designed to solve the following objective: N-1
minJ ) E{eT(N)Me(N) + u(‚)
eT(t)Qe(t) + ∑ t)0 δuT(t)Rδu(t)} (46)
The same weighting matrices as in BLQG were applied. The results of LQG control is shown only for model I. Nominal Case. Figure 2 shows the progressive improvement with the run number of the temperature tracking and the associated input powers from suboptimal BLQG for model I. For simplicity, only the temperature at r ) R/3 and the lamp power p3 are
1668
Ind. Eng. Chem. Res., Vol. 40, No. 7, 2001
Figure 2. Performance of suboptimal BLQG designed with model I for the nominal case.
plotted. Optimal BLQG was implemented for comparison, but the results are not shown as a figure because they are very similar to Figure 2. One can see that very tight tracking can be attained up to cooling only in several runs. For comparison, performance of LQG control is shown in Figure 3. As can be observed, the tracking performance of LQG is inferior to BLQGs, although it outperforms the multiloop PID controller. In Table 2, the filtered maximum tracking errors are compared at each measurement point and also at the point where the maximum error occurs over the whole wafer surface. One can obviously see that BLQGs perform far better than LQGs, whereas the optimal BLQG algorithm shows only a marginal improvement over the suboptimal one at the expense of increased computation and memory usage. For BLQGs, no considerable improvement could be made after 10th run. As the theory dictates in eq 44, the resulting tracking offsets are (close to) the minimum-achievable errors under the chosen quadratic objective. If it is desired to reduce the offsets more, an improvement should be sought from the design side. Perturbed Case. In this case, it was assumed that the view factors of the RTP system are reduced randomly by -5 to 0% by contamination or other process change, and the once-perturbed situation is kept for the remaining runs. Figure 4 shows the results of suboptimal BLQG for two selected runs. The controller designed for the nominal case was used without modification in the simulation. As is observed, the tracking
Figure 3. Performance of LQG designed with model I for the nominal case. Table 2. Comparison of the Filtered Maximum Tracking Errors, maxt |yd(t) - yk(t|t)|K, between BLQG and LQG Based on Model I for the Nominal Case output error during
r)0
r ) R/3
r ) 2R/3
r)R
max0er
