A Multivariate EWMA Controller for Linear Dynamic Processes

Most research of run-to-run process control has been based on single-input and single-output processes with static input–output relationships. In practice, many complicated semiconductor manufacturing processes have multiple-input and multiple-output (MIMO) variables. In addition, the effects of previous process input recipes and output responses on the current outputs might be carried over for several process periods. Under these circumstances, using conventional controllers usually results in unsatisfactory performance. To overcome this, a complicated process could be viewed as dynamic MIMO systems with added general process disturbance and this article proposes a dynamic-process multivariate exponentially weighted moving average (MEWMA) controller to adjust those processes. The long-term stability conditions of the proposed controller are derived analytically. Furthermore, by minimizing the total mean square error (TMSE) of the process outputs, the optimal discount matrix of the proposed controller under vector IMA(1, 1) disturbance is derived. Finally, to highlight the contribution of the proposed controller, we also conduct a comprehensive simulation study to compare the control performance of the proposed controller with that of the single MEWMA and self-tuning controllers. On average, the results demonstrate that the proposed controller outperforms the other two controllers with a TMSE reduction about 32% and 43%, respectively.


INTRODUCTION
Run-to-run (RTR) process control techniques have been widely used in semiconductor manufacturing operations, and their common goal is to bring the process output as close to a desired target as possible (Chen and Guo 2001;Hamby, Kabamba, and Khargonekar 1998;Moyne et al. 2002;Sachs, Hu, and Ingolfsson 1995). In particular, the model-based feedback control scheme has received great attention. This technique consists of a predicted model and the feedback control scheme. First, the predicted model is constructed using off-line stage data. Then, at the end of each production run, to adjust the process output response to a desired target, the feedback control scheme determines the new input recipe for the next run by comparing the output with the target.
Traditionally, statistical process control (SPC, including exponentially weighted moving average, EWMA, and cumulative sum, CUSUM, charts) and engineering process control (EPC, including feedback control schemes) are two well-known techniques that have been commonly used for the purpose of quality control Qiu 2014;Qiu and Li 2011;Qiu and Xiang 2014;Zou and Tsung 2011). Recently, the EWMA statistic has been more proactive to be used in semiconductor manufacturing for adjusting the output of a dynamic RTR process to a desired target. For a single-input and single-output (SISO) system, Ingolfsson and Sachs (1993) and Butler and Stefani (1994) proposed single EWMA (sEWMA) and double EWMA (dEWMA) controllers for adjusting nondrifted and drifted processes, respectively. Del Castillo (1999) and Tseng, Chou, and Lee (2002a) further investigated the long-term sta-bility conditions and transient performance and determined the optimal discount factor of the dEWMA controller.
All of these RTR studies assumed that the dynamic behavior of the output is only due to disturbance dynamics, not to process dynamics. However, in practical processes, the previous input recipes and output responses might have carry-over effects on the current process outputs. Hence, process dynamics and disturbance dynamics might occur simultaneously, which has been demonstrated in relevant literature (Capilla et al. 1999;Jen and Jiang 2008). Pan and Del Castillo (2001) investigated the identification and fine tuning of closedloop dynamic models under discrete sEWMA and proportionalintegral (PI) adjustments. Fan et al. (2002) proposed a triple EWMA controller to compensate the positive autocorrelation that is present in the observed response. Using the sEWMA controller to adjust dynamic processes, Tseng and Lin (2009) investigated the long-term stability conditions of the controller. Furthermore, Tseng and Mi (2014) proposed a quasi minimum mean square error (qMMSE) controller and comprehensively investigated its stability conditions and control performance. The results showed that the qMMSE controller has the ability to adjust the process efficiently when process dynamics exist.
Although these studies reveal the importance of process dynamics, they are restricted to SISO systems. However, many manufacturing processes are multiple-input and multiple-output (MIMO) systems (Del Castillo and Rajagopal 2002;Lee, Chou, and Tseng 2008;Tseng, Chou, and Lee 2002b). Jen, Jiang, and Fan (2004) proposed an MIMO RTR control framework using a self-tuning (ST) controller, which is capable of compensating for chemical mechanical polishing (CMP) processes that exhibit process dynamics and noise disturbance simultaneously. Generally speaking, the long-term stability conditions and shortterm control performance of the ST controller with a recursive least squares (RLS) algorithm are difficult to be analytically addressed. In this study, focusing on a parsimonious transfer function with added general vector autoregressive integrated moving average (ARIMA) disturbance, a dynamic-process multivariate EWMA (dpMEWMA) controller with minimal adjustment is proposed. Next, the theoretical results of the long-term stability conditions are derived and the optimal discount matrix of the proposed controller is investigated. Finally, to assess the contribution of the proposed controller, the control performance of the proposed controller is compared with those of existing controllers, such as the single MEWMA (sMEWMA) and ST controllers, under various combinations of parameter settings. Compared to the sMEWMA controllers, on average, the proposed controller has the ability to reduce the total mean square error (TMSE) by as much as 32%. This demonstrates that the sMEWMA controller, which completely ignores dynamic effects, induces poor control performance. Furthermore, the proposed controller also outperforms the ST controller with a TMSE reduction as high as 43%. This may be because the parameter estimates in ST controller is not able to converge quickly to their true values within a short transient period.
The remainder of this article is organized as follows: Section 2 presents the problem formulation. Section 3 addresses the analytical expression of the process output and the stability conditions of the proposed controller. Section 4 discusses the derivation of the optimal discount matrix. Section 5 illustrates the performance study of the proposed method. Section 6 uses a real example to illustrate the proposed procedure and also provides the results of comparing the control performances of the proposed, sMEWMA, and ST controllers. Finally, Section 7 offers concluding remarks and suggests future research topics. Note that all the appendices are given in the online supplementary material.

PROBLEM FORMULATION
In the following section, matrices, vectors, and scalars are denoted by bold capital letters (such as A), bold lower case letters (such as a), and italic lower case letters (such as a), respectively.
Assuming that y t = (y 1,t , . . . , y n,t ) and x t = (x 1,t , . . . , x m,t ) denote the vectors of n output responses and m input recipes, respectively, a dynamic m × n (m ≥ n) MIMO process is postulated in terms of the transfer function model: where α = (α 1 , . . . , α n ) and B = (b ij ) n×m are the vector and matrix of intercept and slope parameters of the process, respectively; = (ψ ij ) n×n and C = (c ij ) n×m are the matrices that represent the carry-over effects for the output responses and input recipes, respectively; for t = 1, 2, . . . , η t = (η 1,t , . . . , η n,t ) is the n-dimensional disturbance vector, which is assumed to be a general vector ARIMA series. In other words, where p B p and (B) = I n×n − 1 B − · · · − q B q denote the matrix polynomials with autoregressive and moving average parameters, 1 , . . . , p and 1 , . . . , q , respectively; the roots of det( (B)) = 0 are outside the unit circle; ε t is independent of each time point with mean vector 0 n×1 and covariance matrix = (σ ij ) n×n . In addition, all the eigenvalues of are less than 1 to ensure the stationarity of the underlying model.
At the offline stage, let a 0 ,B = (b ij ) n×m ,Ĉ = (ĉ ij ) n×m , and = (ψ ij ) n×n are the estimates of α, B, C, and , respectively. These estimates can be obtained by a time-series transfer function model and hereafter, we treat these estimates as fixed constants under RTR control setting. Similar to the ordinary EWMA controller of Ingolfsson and Sachs (1993), the following controller is proposed to recursively update the estimates of the process intercepts: where = diag (ω 1 , . . . , ω n ) denotes a diagonal discount matrix. Furthermore, to bring the process outputs to the desired targets with the smallest adjustment between consecutive input recipes, the new input recipes are updated as follows: where τ = (τ 1 , . . . τ n ) is the vector of the desired targets. Using the Lagrange multiplier method, a recursive updating formula for the input recipes is obtained: A controller satisfying Equations (3) and (4) is called a dp-MEWMA controller. From Equation (4), the following result is obtained: The proofs of Equations (4) and (5) are given in Appendix 1 (see online supplementary materials). Recursively iterating x t − x t−1 in Equation (5) shows that the proposed controller uses not only the current data, y t , but also the previous data, {y j } t−1 j =1 , to update the input recipes. The proposed controller reduces to the conventional sMEWMA controller ifˆ = 0 n×n andĈ = 0 n×m . Compared to the sMEWMA controller, Equations (3) and (4) show that the proposed controller can efficiently eliminate dynamic effects if the true I-O model is dynamic. Moreover, as n = m = 1, Equations (1) and (4) can be expressed as and where B is a backshift operator. Equation (7) is essentially equivalent to the qMMSE controller (Tseng and Mi 2014), which is constructed using the MMSE approach. For MIMO processes, a controller is asymptotically stable if the expectations of the process outputs converge to the desired targets and the corresponding asymptotic covariance matrix remains bounded: and Equations (8) and (9) are known as stability conditions. Issues regarding the proposed dpMEWMA controller include: (a) What are the conditions that lead the process outputs controlled by the dpMEWMA controller to be convergent to the desired targets? (b) How to determine optimal discount matrix of the dp-MEWMA controller so that its TMSE can be minimized?
Both issues will be discussed in details in Sections 3 and 4, respectively, and also be illustrated through two examples in Sections 5 and 6.

STABILITY CONDITIONS OF THE PROPOSED CONTROLLER
Before studying the stability conditions of the dpMEWMA controller, the following notations must be considered: Let where Furthermore, for all t ≥ 1, let where E 0 = I n×n , F 0 = 0 n×n , and G 0 = 0 n×m .
The following lemma shows that the deviances between the process outputs and the targets, y t − τ , can be partitioned into two parts: a nonrandom part, γ t−1 , and a random part, ξ t .
Lemma 1. The process outputs controlled by the dpMEWMA controller can be expressed as follows: where and γ 0 = α + Bx 0 − τ is the vector of the initial biases. Appendix 2 (see online supplementary materials) provides the proof of Lemma 1. Lemma 1 provides a very useful result for studying the stability conditions of the dpMEWMA controller. It reveals that the stability conditions can be attained when lim t→∞ γ t = 0 n×1 and the limit of the covariance matrix of ξ t is bounded.
Define the maximum absolute values of the eigenvalues of matrix M as follows: and then the stability conditions are derived and stated in the following theorem: Theorem 1. Assuming that η t satisfies Equation (2), the dp-MEWMA controller is asymptotically stable if max 1≤i≤n d i ≤ 1 and ρ (M) < 1.
Appendix 3 (see online supplementary materials) provides the proof of Theorem 1.
Note that settingˆ = 0 n×n andĈ = 0 n×m in M to obtain M s leads Then, substituting M s for M in Equation (16) can directly obtain the stability conditions of the sMEWMA controller. In Section 5, we will use an example to demonstrate how to apply Theorem 1 for finding the feasible regions of dpMEWMA and sMEWMA controllers, respectively. Also, the discussion of required sample size of historical data would also be investigated using the same example.
Except for deriving the stability of dpMEWMA controller, we may also interest in deriving the exact distribution of ξ t . Generally speaking, it is not easy to derive an analytical expression of ξ t in Equation (12). However, if we imposeĈ = kB, with |k| < 1, it allows us to have a closed form on ξ t . Note that this condition is often reasonable in practice. To ensure a convergent pattern in Equation (16), we usually expect the dynamic carry-over effect in input variables will decay gradually. In that case, we can assume that the element inĈ is roughly proportional to the element inB. Now, underĈ = kB, by using matrix operation techniques, we have the following result: Appendix 4 (see online supplementary materials) provides the proof of Equation (17). This result will be very useful in deriving the optimum discount matrix for the dpMEWMA controller which will state in the next section.

OPTIMAL DISCOUNT MATRIX OF THE PROPOSED CONTROLLER
In relevant literature, the TMSE of the process outputs is usually used to evaluate the control performance of controllers, and is defined as follows: where N is assumed to be moderately large. Based on Lemma 1, Equation (18) can be expressed as follows: The first term of Equation (19) comes from the nonrandom part, which is irrelevant to the process disturbances. As the stability conditions are attained, then For the random part, from Equation (17), ξ t follows a stationary vector ARMA (3 + p, 2 + q − d) process if d = max{d 1 , . . . , d n } ≤ 1. Therefore, Equation (19) can be further expressed as follows: where υ = trace cov ξ t .

PERFORMANCE STUDY
In the following, one simulated dataset is used to illustrate how to derive the stability conditions and the optimal discount matrix of the proposed controller, respectively.

Illustration for the Stability of the dpMEWMA Controller
We slightly modify example 9.6 of Del Castillo (2002, p. 317) with the settings for considering only a single dynamic term on both the input and output variables. That is, the process I-O model comes from Equation (1) For example, if (ω 1 , ω 2 ) = (0.4, 0.6), then we have ρ (M) = 0.62 < 1. From Theorem 1, the process outputs adjusted by dpMEWMA controller with discount matrix diag(0.4, 0.6)will satisfy the stability conditions. Furthermore, let F = {(ω 1 , ω 2 )|ρ (M) < 1} denote the feasible region for the discount matrix of the dpMEWMA controller and it can be shown that F = (0, 1) × (0, 1) . Thus, it achieves global stability. By using the same argument, witĥ . This result demonstrates that the dpMEWMA controller is globally stable, but the feasible region of the sMEWMA controller will be relatively narrow to control the dynamic process and it thus takes much more runs to adjust the output to the desired target.

How
Large the Sample Size Is Needed to Achieve Global Stability?
In previous Section 5.1, the predicted model is developed based on the offline data. Therefore, the sample size and the strength of I-O relationship play major roles in determining the stability of the process. For the static SISO and MIMO models, Tseng and Hsu (2005) and Tseng, Tang, and Lin (2007) derived simple formulas for the adequate sample size to achieve stability with a guaranteed probability, respectively. Generally speaking, it is not easy to obtain an analytical result for the case of dynamic MIMO models. Instead, we use a simulation study to address this decision problem. In the following study, we assume that all the elements in matrices , B, and (stated in Section 5.1) are proportional to its original scale with a uniform random variable. That is, where f ∼ unif(0.9, 1.1). Under various settings of the sample size at the offline stage N 0 = 30, 40, . . . , 100,we implement 1000 simulation trials from the process with this new parameter settings. Figure 1 shows the relationship of the proportions of achieving global stability (among 500 simulation trials) with respect to sample size.
For example, when N 0 = 50, we can see 430 trials achieved global stability. Figure 1 also demonstrates the monotonic pattern that the higher the proportion of achieving global stability, the larger sample size is required. Roughly speaking, we will have a 97.5% confidence to achieve the global stability when N 0 ≥ 70.

Illustration for the Constrained Optimal Discount
Matrix of the dpMEWMA Controller Under the same parameter settings mentioned in Section 5.1 and further imposing the condition ofĉ ij = k ×b ij , for all 1 ≤ i ≤ n, and 1 ≤ j ≤ m, the corresponding estimators of these parameters are: Assume that the desired targets for output responses are τ 1 = 0 andτ 2 = 0. By using the formulas stated in Appendices 5 and 6 (see online supplementary materials) and for the total run size N = 50, Figure 2 shows the contour plot of TMSE under various combinations of (ω 1 , ω 2 ) , where 0 < ω 1 < 1 and 0 < ω 2 < 1.
The cross-point of the two solid lines in the figure represents the optimal TMSE * (= 889.48, in which the nonrandom part equals 763.14 and the random part equals126.34) and the corresponding constrained optimal discount factors (ω * 1 = 0.40, ω * 2 = 0.92) can be traced directly along the lines. That is, * = 0.40 0 0 0.92 .
Furthermore, a scale discount matrix, 0 = 0.20 × I 2×2 , is also shown in Figure 2 as the cross-point of two dotted lines. The corresponding nonrandom and the random parts of the TMSE are 2026.59 and 122.14, respectively, resulting in the corresponding TMSE = 2148.74. Comparing these two TMSE values, the proposed controller with the constrained optimal discount factors has the ability of reducing TMSE roughly up to 59%. Figure 3 plots the sequences of process outputs controlled by 0 = 0.2 × I 2×2 and * . The result demonstrates that the dpMEWMA controller with the constrained optimal discount matrix regulates the process outputs more efficiently than that of using nonoptimal discount matrix.

REAL EXAMPLE
A real application of CMP process appearing in RTR literature is used to illustrate the control performance of the proposed dp-MEWMA controller. CMP process is a planarization technique that removes the roll from the surface of a wafer, which is an important processing step in semiconductor manufacturing processes. Del Castillo and Yeh (1998) proposed a 4 × 2 MIMO process aimed to maintain process outputs on targets. In this dataset, four input recipes: plate speed, back pressure, polishing down-force, and profile of conditioning system at time t are denoted by x 1,t , x 2,t , x 3,t , and x 4,t , respectively; while two output responses: the within wafer nonuniformity and removal rate at time t are denoted by y 1,t and y 2,t , respectively. Assume that the desired targets for output responses are τ 1 = 100 and τ 2 = 2000.
For illustrating our proposed procedure, we adopt the dynamic process (see Appendix 8 in the online supplementary Figure 3. The sequences of process outputs controlled by the dpMEWMA controller with two discount matrices 0 and * .  (22) are actually unknown and we need to estimate these 30 unknown parameters at the offline stage. Obviously, with too small observations at the offline stage, it is impossible to get a reasonable estimation of the prediction model. Therefore, we first need to address the effect of the offline observations on the controller performance of sMEWMA and dpMEWMA controllers.

The Effect of the Prediction Model on the Control Performance
In the following study, three different settings, say T = 60, 80, 100 observations were simulated from the uncontrolled process in Equation (22) and the prediction models for the static and dynamic systems arê and y t = ⎧ ⎨ ⎩ˆ 1 y t−1 + a 01 +B 1 x t−1 +Ĉ 1 x t−2 , T = 60 2 y t−1 + a 02 +B 2 x t−1 +Ĉ 2 x t−2 , T = 80 3 y t−1 + a 03 +B 3 x t−1 +Ĉ 3 x t−2 , T = 100, whereˆ i , a 0i ,B i ,Ĉ i , i = 1, 2, 3are shown in Table 1. The estimates of unknown parameters in Table 1 were essentially computed by the subroutine "pem" of the system identifi- cation toolbox of Matlab (Del Castillo 2002, chap. 9.3). However, due to the typical "nonidentifiable" problem of MIMO transfer function model, it may lead the estimates ofˆ and C to be confounded and unstable. Instead of incorporating the complex identifiability constraints on parameters, a two-stage approach is proposed to solve this problem. First, we use regression technique to obtain a reasonable estimation forˆ and then the other parameters can be estimated by the Matlab subroutine "pem" directly. Furthermore, the norm in Table 1 (Horn and Johnson 1990) is defined as It can be used to measure the estimated accuracy ofĤ with respect toH. From Table 1, it demonstrates that the estimated accuracies under T = 80 and T = 100 are significantly better than that of T = 60.
With the above prediction models, we can use Equations (3) and (4) to recursively adjust the process outputs for the sMEWMA and dpMEWMA controllers. The ST controller is also commonly used for controlling dynamic processes in literature (Del Castillo and Hurwitz 1997). Therefore, it is sensible to compare the performance of those controllers with that of the dpMEWMA controller. To ensure a fair comparison of control performances, all three controllers adopt their corresponding optimal discount matrices (by using the complete enumeration search over all of their feasible regions). For example, when T = 60 and the total production run size N = 50, we adopt the best discount matrices, * dp = [ 0.45 0 0 0.15 ] and * s = [ 0.21 0 0 0.01 ], for the dpMEWMA and sMEWMA controllers, and (λ * , ρ * ) = (0.51, 0.94) for the ST controller (derived in Appendix 9; see online supplementary materials). We define the relative efficiencies (REs) to measure the control performances of the proposed controller with respect to sMEWMA and ST controllers as follows: and  Under T = 60, 80, 100, the values of R 1 are 0.76, 0.77, and 0.72; and the values of R 2 are 0.24, 0.67, and 0.61, respectively. It means that the dpMEWMA controller has the ability to reduce the process TSSE at least 23%, compared with the sMEWMA and the ST controllers. In the meantime, it also demonstrates that the offline sample size does not have a significant effect on the control performance in this study. Figure 4 also plots the corresponding sequences of process outputs of all three controllers under T = 60, 80, 100. It demonstrates that the process outputs controlled by the sMEWMA and ST controllers have larger variations than that of the dpMEWMA controller.

Comprehensive Comparisons With Two Competing Controllers
In the above CMP example, we compare the control performances of the proposed controller with respect to those competing controllers that are only based on a single simulation trial. Note that the control performance strongly depends on the estimated accuracy of the process I-O prediction model, which depends not only on the sample size at the offline stage, but also on the covariance structure ( ) of ε t and the dynamic matrix ( ). Therefore, we will further discuss these issues via a com-prehensive simulation study. In the following, besides for T = 60, 80, 100, two different levels for and are considered as follows: The main differences in these settings are that the off-diagonal elements of 2 have larger effects than that of 1 , while the elements of 2 have smaller covariance than that of 1 . For each combination, we conduct 500 simulation trials to address the REs of the proposed dpMEWMA controllers with respect to sMEWMA and ST controllers. In each simulation trial, we follow the procedure stated in Section 6.1 to compute the REs of these three controllers. For all 1 ≤ k ≤ 500, let R 1 (k) and R 2 (k) (Equations (25) and (26)) denote the REs of the proposed controller with respect to sMEWMA and ST controllers under the kth simulation trial. Their overall means and standard deviations (SD) can be computed by 500 k=1 R i (k)/500 and 500 , and the offline observations with {T = 60, T = 80, T = 100}. Compared with the sMEWMA controller, it demonstrates that the overall RE performances of the dpMEWMA with respect to sMEWMA controllers (i.e., under = 1 , and = 1 ) are very similar to the results of a single simulation trial in Section 6.1. Moreover, it also demonstrates that the REs can be slightly improved when 2 is adopted. Furthermore, compared with ST controller, the means of REs are ranging from 0.53 to 0.59 for all combinations of , and the offline observations with {T = 60, T = 80, T = 100}.
Figures 5(a) and 5(b) also show the boxplots of these REs performances in Table 2. The plots demonstrate that the RE performances are relatively stable for all combinations in this study. In addition, all potential outliers (with only one exception) in these plots occurred at very smaller values of REs. From Equations (25) and (26), it means that the process outputs controlled by sMEWMA and ST controllers may cause relatively large variations, compared with dpMEWMA controller. The reasons may be mainly due to the facts that the dynamic effects of process I-O variables are not well-treated by sMEWMA controller; while the parameter estimates in ST controller may not be able to converge quickly to their true values within a short transient period.

CONCLUSION
This study proposed a dpMEWMA controller for adjusting linear dynamic processes with added general process disturbance. Under the assumed-true process, the stability conditions of the proposed controller are derived, and the control performance and the optimal discount matrix are analytically investigated. Furthermore, based on the simulation study, the control performance of the proposed controller is compared with that of the sMEWMA and ST controllers. On average, the results demonstrate that the proposed controller outperforms the other two controllers with a TMSE reduction about 32% and 43%, respectively.
Further remarks and opportunities for future research are listed as follows: 1. What's the performance if the true process does not follow the assumed model?
In practice, the true process I-O model may be more general than Equation (1). Adding second-order dynamic term and linear drift to Equation (1), we assume the true process I-O model as follows: Obviously, how to address the effects of model misspecification on control performance (of the proposed controller) turns out to be an important research topic. Tseng and Mi (2014) addressed this issue under the case of the SISO model. They compared the performance of the proposed controller with the competing controllers and the results demonstrate that the proposed controller outperforms the competing controllers even when the process I-O model is misspecified. For the MIMO case, it shall be an interesting, yet challenging, research topic for further research. 2. For an MIMO system, the proposed controller is constructed using the EWMA statistics with minimal recipe adjustment. When the model parameters are perfectly estimated and = I n×n − , Equation (4) reduces to Note that this controller has a different expression with the MMSE controller, which is stated in Reinsel (2003), that is, The difference is induced by using the constraint of minimal adjustment. However, from the derived TMSE formula in Equation (20), it demonstrates that the dp-MEWMA controller is essentially an MMSE controller.
3. The MEWMA controller (stated in Equation (3)) can be easily reduced to the following equation: a t = a t−1 + (y t − τ ) .
It has a very good interpretation that the amount of adjustment in the intercept vector is proportional to the difference from the output vector to a desired target vector. This may be the reason that the MEWMA controller is a popular model-based controller in RTR control. Note that CUSUM and EWMA are two popular SPC charts. Qiu (2014) made a comprehensive comparison between CUSUM charts and EWMA charts in his new book. However, there is no suitable controller which is constructed by CUSUM statistic currently. Therefore, how to construct a suitable CUSUM-based controller shall be a challenging issue for future research. 4. Similar to Tseng and Jou (2009), the proposed controller can be easily constructed by considering the estimate of the covariance matrix of input variables (S x ) when searching for the smallest adjustment between consecutive input recipes, that is, subject toˆ y t + a t +Bx t +Ĉx t−1 = τ .
Furthermore, given a specific value ofρ 0 , we can also consider the following formulation: subject to ˆ y t + a t +Bx t +Ĉx t−1 − τ < ρ 0 .
These two approaches would be benefit for saving the offline-stage sample size, which are worthy of future investigation. 5. The results of this study are based on the assumption that the process response variables can be measured immediately by metrology systems. However, many critical processes might encounter the situations of metrology delay (Lin, Tseng, and Wang 2013;Wu et al. 2008). The effects of metrology delay on the stability conditions and transient behavior of the proposed controller are of great practical importance and worthy of further research.

SUPPLEMENTARY MATERIALS
1. All the proofs and the process I-O Model (used in Section 6) are given in Appendices. (Appendices.pdf) 2. Matlab source code for the outputs of Sections 5.3 and 6.1. (matlab codes.pdf)