CoCalc Public FilesSPC test.ipynbOpen in with one click!
Author: phonchi chung
Views : 8
Description: Jupyter notebook SPC test.ipynb
In [ ]:
import numpy as np import matplotlib.pyplot as plt %matplotlib inline

Engineering Process Control and Statistical Process Control

  1. Do not know how to estimate AkA_k and TgtyTgty => Ak = ARRK[599] or linear regressioon on past data, TgtY = PreY - TgtPostY
  2. Should I use AkA_k inside ARRkARR_k???? => Ak0 = const
  3. Assume ARRv=0.9*ARRk + variation from AVM, RI = 0.9 from AVM
  4. Use const Tgt in uk instead of tgtk+1 in paper? => Still use tgtk+1 now
  5. Should I add tou to process output? => Now add into processoutput and vm

ARRv is AVM's output, ARRk is atual measurement value, and Ak is used internal by R2R

The basic Runto Run Controller

β0\beta_0 and β1\beta_1are assumed to be constant over time??? They are unknown and are to be estimated from available data.

The W2W control scheme:

  • The yky_k is either form AVM or Metrology tool
  • AkA_k (or A) is typically chosen to be least square estimates of β1β_1 based on historical data.The control valuable is set to nullify the deviation from target.

CMP Example

  1. yky_k is the actual removal amount measured from the metrology tool and PostYkPostY_k is the actual post CMP thickness of run k. The specification of PostYkPostY_k is 2800±150 Angstrom (Å) with 2800 being the target value denoted by TgtPostYTgtPostY Deifned PreYkPreY_k be the one before CMP process, ARRkARR_k the polish rate, μk\mu_k the polish time that we can control!

The material removal model for CMP can be divided into two parts, mechanical model and chemical model. The chemical action of slurry is responsible for continuosly softening the silicon oxide. The fresh silicon oxide or metal surface is then rapidly removed by mechanical part.

  1. The AkA_k is the nominal removal rate, which is empirically simulated by a polynomial curve fitting of parts usage count between PMs (denoted by PU varying from 1 to 600)

The estimation of A(beta_1) and beta_0 is based on linear regression

  1. Process gain (Ak or ARRk):
  1. Simulation parameters:

PU,PU2,PU3PU, PU^2, PU^3 -> empirical parameters

simulated parameters:

In [ ]:
from statsmodels.tsa.arima_model import ARIMA
In [ ]:
Ak0 = 300 # Not sure, does not know PU
In [ ]:
# Simulated environment Virtual = np.random.normal(1, np.sqrt(0.01), 600) Error = np.random.normal(0, np.sqrt(300), 600) PM1 = np.random.normal(0, np.sqrt(100), 600) PM2 = np.random.normal(0, np.sqrt(6), 600) Stress1 = np.random.normal(1000, np.sqrt(2000),600) Stress2 = np.random.normal(0, np.sqrt(20), 600) Rotspd1 = np.random.normal(100, np.sqrt(25), 600) Rotspd2 = np.random.normal(0, np.sqrt(1.2), 600) Sfuspd1 = np.random.normal(100, np.sqrt(25), 600) Sfuspd2 = np.random.normal(0, np.sqrt(1.2), 600) PreY = np.random.normal(3800, np.sqrt(2500),600) ListY = [2800]*600 TgtPostY = np.array(ListY) ### Should I use Ak inside here? def ARR (stress1, stress2, rotspd1, rotspd2, sfuspd1, sfuspd2, pm1, pm2): return Ak0*((stress1+stress2)/1000)*((rotspd1+rotspd2)/100)*((sfuspd1+sfuspd2)/100)+pm1+pm2#+error ARRk = ARR(Stress1, Stress2, Rotspd1, Rotspd2, Sfuspd1, Sfuspd2, PM1, PM2) # beta1 actual process gain #ARRv = ARRk*Virtual ## Hypothetic value of VM #Ak = ARRk*0.95 ## Hypothetic value of VM Ak = np.repeat(ARRk[599],600) TgtY = np.zeros(len(PreY)) TgtY = PreY - TgtPostY
In [ ]:
model = ARIMA(Error, order=(0, 1, 1)) results_AR = model.fit(disp=-1) results_AR.fittedvalues

Round 1

First two rounds, we use actual metrodlogy value

In [ ]:
#TgtY = 1000 # Not sure assume etch 1000A, since the mean of PreY is 3800A and target PostY is 2800A stated in the paper alpha1 = 0.35 Toub = np.zeros(len(PreY)) mu = np.zeros(len(PreY)) PostY = np.zeros(len(PreY)) yz = np.zeros(len(PreY)) #Ak = ARRk[599]
In [ ]:
Ak = np.array([300]*600) #>300 diverge? mu[0] = (TgtY[0] - Toub[0])/Ak[0] # k=1 yz[0] = ARRk[0]*mu[0]+200 + results_AR.fittedvalues[0] # Actual value PostY[0] = PreY[0] - (ARRk[0]*mu[0]+200) Toub[1] = alpha1*(yz[0]-Ak[0]*mu[0])+(1-alpha1)*Toub[0] mu[1] = (TgtY[1] - Toub[1])/Ak[1]

Round2

In [ ]:
# k=2[0] # Toub[1] =0 yz[1] = ARRk[1]*mu[1]+200 + results_AR.fittedvalues[1] # Actual value PostY[1] = PreY[1] - (ARRk[1]*mu[1]+200) Toub[2] = alpha1*(yz[1]-Ak[1]*mu[1])+(1-alpha1)*Toub[1] mu[2] = (TgtY[2] - Toub[2])/Ak[2]

CASE1 R2R with in-situ metrology

set to 0.35, and all actual metrology data are available

In [ ]:
for k in range(2,599): yz[k] = ARRk[k]*mu[k]+200 + results_AR.fittedvalues[k] ## ARRk[k], Toub[k] approximate beta_1, beta_0(not tou!!) #PostY[k] = PreY[k] - yz[k] PostY[k] = PreY[k] - (ARRk[k]*mu[k]+200) Toub[k+1] = alpha1*(yz[k]-Ak[k]*mu[k])+(1-alpha1)*Toub[k] #mu[k+1] = (1000 - Toub[k+1])/Ak[k+1] mu[k+1] = (TgtY[k+1] - Toub[k+1])/Ak[k+1] x = range(10, 590) U = np.empty(580) L = np.empty(580) UCL = 2950 LCL = 2650 U.fill(2950) L.fill(2650) # plot(x, PostY[3:600], type="o") plt.plot(x,U) plt.plot(x,L) plt.plot(x,PostY[10:590]) np.mean(PostY[10:590]) #plt.ylim([2600, 3000])
In [ ]:
x = range(10, 590) plt.plot(x,Toub[10:590])
In [ ]:
x = range(10, 590) plt.plot(x,yz[10:590])

CpK(Process Capability)

In [ ]:
Cpk=np.zeros(len(PreY)) MAPEp =np.zeros(len(PreY))
In [ ]:
for k in range(3,599): Cpk[k] = min((UCL-np.mean(PostY[1:k]))/(3*np.std(PostY[1:k], dtype=np.float64)), (np.mean(PostY[1:k])-LCL)/(3*np.std(PostY[1:k], dtype=np.float64))) plt.plot(x,Cpk[10:590])

Mean-absolute-percentage error (MAPEp)

In [ ]:
for k in range(3,599): MAPEp[k] = sum(np.absolute((PostY[1:k]-2800)/2800))/k*100 plt.plot(x,MAPEp[10:590])

CASE2 R2R+VM without RI

, apply to following wafers

In [ ]:
#TgtY = 1000 # Not sure assume etch 1000A, since the mean of PreY is 3800A and target PostY is 2800A stated in the paper alpha1 = 0.35 Toub = np.zeros(len(PreY)) mu = np.zeros(len(PreY)) PostY = np.zeros(len(PreY)) PostYv = np.zeros(len(PreY)) y = np.zeros(len(PreY)) yz = np.zeros(len(PreY)) #Ak = ARRk[599] Ak = np.array([280]*600) #>300 diverge? mu[0] = (TgtY[0] - Toub[0])/Ak[0] # k=1 yz[0] = ARRk[0]*mu[0]+results_AR.fittedvalues[0]+200 # Actual value PostY[0] = PreY[0] - (ARRk[0]*mu[0]+200) Toub[1] = alpha1*(yz[0]-Ak[0]*mu[0])+(1-alpha1)*Toub[0] mu[1] = (TgtY[1] - Toub[1])/Ak[1] # k=2[0] # Toub[1] =0 yz[1] = ARRk[1]*mu[1]+results_AR.fittedvalues[1]+200 # Actual value PostY[1] = PreY[1] - (ARRk[1]*mu[1]+200) Toub[2] = alpha1*(yz[1]-Ak[1]*mu[1])+(1-alpha1)*Toub[1] mu[2] = (TgtY[2] - Toub[2])/Ak[2]
In [ ]:
for k in range(2,599): yz[k] = ARRk[k]*mu[k]+results_AR.fittedvalues[k]+200 # Metrology data if(k%25)==1: y[k] = yz[k] else: y[k] = yz[k]*Virtual[k] # Assume VM value ok ARR=ARR_bar # Actual PostY[k] = PreY[k] - (ARRk[k]*mu[k]+200) PostYv[k] = PreY[k] - (yz[k]*Virtual[k]) Toub[k+1] = alpha1*(y[k]-Ak[k]*mu[k])+(1-alpha1)*Toub[k] mu[k+1] = (TgtY[k+1] - Toub[k+1])/Ak[k+1]
In [ ]:
x = range(10, 590) U = np.empty(580) L = np.empty(580) U.fill(2950) L.fill(2650) # plot(x, PostY[3:600], type="o") plt.plot(x,U) plt.plot(x,L) plt.plot(x,PostY[10:590]) np.mean(PostY[10:590])
In [ ]:
Cpk=np.zeros(len(PreY)) MAPEp =np.zeros(len(PreY))
In [ ]:
for k in range(3,599): Cpk[k] = min((UCL-np.mean(PostY[1:k]))/(3*np.std(PostY[1:k], dtype=np.float64)), (np.mean(PostY[1:k])-LCL)/(3*np.std(PostY[1:k], dtype=np.float64))) plt.plot(x,Cpk[10:590])
In [ ]:
for k in range(3,599): MAPEp[k] = sum(np.absolute((PostY[1:k]-2800)/2800))/k*100 plt.plot(x,MAPEp[10:590])

CASE3: R2R+VM with RI

In [ ]:
#TgtY = 1000 # Not sure assume etch 1000A, since the mean of PreY is 3800A and target PostY is 2800A stated in the paper alpha1 = 0.35 Toub = np.zeros(len(PreY)) mu = np.zeros(len(PreY)) PostY = np.zeros(len(PreY)) PostYv = np.zeros(len(PreY)) y = np.zeros(len(PreY)) yz = np.zeros(len(PreY)) #Ak = ARRk[599] Ak = np.array([280]*600) #>300 diverge? mu[0] = (TgtY[0] - Toub[0])/Ak[0] # k=1 yz[0] = ARRk[0]*mu[0]+results_AR.fittedvalues[0]+200 # Actual value PostY[0] = PreY[0] - (ARRk[0]*mu[0]+200) Toub[1] = alpha1*(yz[0]-Ak[0]*mu[0])+(1-alpha1)*Toub[0] mu[1] = (TgtY[1] - Toub[1])/Ak[1] # k=2[0] # Toub[1] =0 yz[1] = ARRk[1]*mu[1]+results_AR.fittedvalues[1]+200 # Actual value PostY[1] = PreY[1] - (ARRk[1]*mu[1]+200) Toub[2] = alpha1*(yz[1]-Ak[1]*mu[1])+(1-alpha1)*Toub[1] mu[2] = (TgtY[2] - Toub[2])/Ak[2]
In [ ]:
RI = 0.9 # Assume value of VM alpha2 = RI*alpha1 for k in range(2,599): yz[k] = ARRk[k]*mu[k]+results_AR.fittedvalues[k]+200 # process output how to filter out???? if(k%25)==1: y[k] = yz[k] alpha = alpha1 else: y[k] = (ARRk[k]*mu[k]+200)*Virtual[k] alpha = alpha2 PostYv[k] = PreY[k] - (ARRk[k]*mu[k]+200)*Virtual[k] if PostYv[k]>2950 or PostYv[k]<2650: alpha = 0 Toub[k+1] = alpha*(y[k]-Ak[k]*mu[k])+(1-alpha)*Toub[k] #Toub canceled out mu[k+1] = (TgtY[k+1] - Toub[k+1])/Ak[k+1] PostY[k] = PreY[k] - (ARRk[k]*mu[k]+200) # Process output else: Toub[k+1] = alpha*(y[k]-Ak[k]*mu[k])+(1-alpha)*Toub[k] #Toub canceled out mu[k+1] = (TgtY[k+1] - Toub[k+1])/Ak[k+1] PostY[k] = PreY[k] - (ARRk[k]*mu[k]+200) # Process output
In [ ]:
x = range(10, 590) U = np.empty(580) L = np.empty(580) U.fill(2950) L.fill(2650) # plot(x, PostY[3:600], type="o") plt.plot(x,U) plt.plot(x,L) plt.plot(x,PostY[10:590]) np.mean(PostY[10:590])
In [ ]:
Cpk=np.zeros(len(PreY)) MAPEp =np.zeros(len(PreY))
In [ ]:
for k in range(3,599): Cpk[k] = min((UCL-np.mean(PostY[1:k]))/(3*np.std(PostY[1:k], dtype=np.float64)), (np.mean(PostY[1:k])-LCL)/(3*np.std(PostY[1:k], dtype=np.float64))) plt.plot(x,Cpk[10:590])
In [ ]:
for k in range(3,599): MAPEp[k] = sum(np.absolute((PostY[1:k]-2800)/2800))/k*100 plt.plot(x,MAPEp[10:590])

CASE4: R2R+VM with (1-RI)

In [ ]:
#TgtY = 1000 # Not sure assume etch 1000A, since the mean of PreY is 3800A and target PostY is 2800A stated in the paper alpha1 = 0.35 Toub = np.zeros(len(PreY)) mu = np.zeros(len(PreY)) PostY = np.zeros(len(PreY)) PostYv = np.zeros(len(PreY)) y = np.zeros(len(PreY)) yz = np.zeros(len(PreY)) #Ak = ARRk[599] Ak = np.array([280]*600) #>300 diverge? mu[0] = (TgtY[0] - Toub[0])/Ak[0] # k=1 yz[0] = ARRk[0]*mu[0]+results_AR.fittedvalues[0]+200 # Actual value PostY[0] = PreY[0] - (ARRk[0]*mu[0]+200) Toub[1] = alpha1*(yz[0]-Ak[0]*mu[0])+(1-alpha1)*Toub[0] mu[1] = (TgtY[1] - Toub[1])/Ak[1] # k=2[0] # Toub[1] =0 yz[1] = ARRk[1]*mu[1]+results_AR.fittedvalues[1]+200 # Actual value PostY[1] = PreY[1] - (ARRk[1]*mu[1]+200) Toub[2] = alpha1*(yz[1]-Ak[1]*mu[1])+(1-alpha1)*Toub[1] mu[2] = (TgtY[2] - Toub[2])/Ak[2]
In [ ]:
RI = 0.9 # Assume value of VM alpha2 = (1-RI)*alpha1 for k in range(2,599): yz[k] = ARRk[k]*mu[k]+results_AR.fittedvalues[k]+200 # process output how to filter out???? if(k%25)==1: y[k] = yz[k] alpha = alpha1 else: y[k] = (ARRk[k]*mu[k]+200)*Virtual[k] alpha = alpha2 PostYv[k] = PreY[k] - (ARRk[k]*mu[k]+200)*Virtual[k] if PostYv[k]>2950 or PostYv[k]<2650: alpha = 0 Toub[k+1] = alpha*(y[k]-Ak[k]*mu[k])+(1-alpha)*Toub[k] #Toub canceled out mu[k+1] = (TgtY[k+1] - Toub[k+1])/Ak[k+1] PostY[k] = PreY[k] - (ARRk[k]*mu[k]+200) # Process output else: Toub[k+1] = alpha*(y[k]-Ak[k]*mu[k])+(1-alpha)*Toub[k] #Toub canceled out mu[k+1] = (TgtY[k+1] - Toub[k+1])/Ak[k+1] PostY[k] = PreY[k] - (ARRk[k]*mu[k]+200) # Process output
In [ ]:
x = range(10, 590) U = np.empty(580) L = np.empty(580) U.fill(2950) L.fill(2650) # plot(x, PostY[3:600], type="o") plt.plot(x,U) plt.plot(x,L) plt.plot(x,PostY[10:590]) np.mean(PostY[10:590])
In [ ]:
Cpk=np.zeros(len(PreY)) MAPEp =np.zeros(len(PreY))
In [ ]:
for k in range(3,599): Cpk[k] = min((UCL-np.mean(PostY[1:k]))/(3*np.std(PostY[1:k], dtype=np.float64)), (np.mean(PostY[1:k])-LCL)/(3*np.std(PostY[1:k], dtype=np.float64))) plt.plot(x,Cpk[10:590])
In [ ]:
for k in range(3,599): MAPEp[k] = sum(np.absolute((PostY[1:k]-2800)/2800))/k*100 plt.plot(x,MAPEp[10:590])

CASE5: R2R+VM with RI.S.(1-RI)

In [ ]:
#TgtY = 1000 # Not sure assume etch 1000A, since the mean of PreY is 3800A and target PostY is 2800A stated in the paper alpha1 = 0.35 Toub = np.zeros(len(PreY)) mu = np.zeros(len(PreY)) PostY = np.zeros(len(PreY)) PostYv = np.zeros(len(PreY)) y = np.zeros(len(PreY)) yz = np.zeros(len(PreY)) #Ak = ARRk[599] Ak = np.array([280]*600) #>300 diverge? mu[0] = (TgtY[0] - Toub[0])/Ak[0] # k=1 yz[0] = ARRk[0]*mu[0]+results_AR.fittedvalues[0]+200 # Actual value PostY[0] = PreY[0] - (ARRk[0]*mu[0]+200) Toub[1] = alpha1*(yz[0]-Ak[0]*mu[0])+(1-alpha1)*Toub[0] mu[1] = (TgtY[1] - Toub[1])/Ak[1] # k=2[0] # Toub[1] =0 yz[1] = ARRk[1]*mu[1]+results_AR.fittedvalues[1]+200 # Actual value PostY[1] = PreY[1] - (ARRk[1]*mu[1]+200) Toub[2] = alpha1*(yz[1]-Ak[1]*mu[1])+(1-alpha1)*Toub[1] mu[2] = (TgtY[2] - Toub[2])/Ak[2]
In [ ]:
RI = 0.9 alpha2 = RI*alpha1 for k in range(2,599): yz[k] = ARRk[k]*mu[k]+results_AR.fittedvalues[k]+200 if(k%25)==1: y[k] = yz[k] alpha = alpha1 else: y[k] = (ARRk[k]*mu[k]+200)*Virtual[k] alpha = alpha2 PostYv[k] = PreY[k]- (ARRk[k]*mu[k]+200)*Virtual[k] if k>= 25: if PostYv[k]>2950 or PostYv[k]<2650: # should be replace with RI, GSI alpha = 0 Toub[k+1] = alpha*(y[k]-Ak[k]*mu[k])+(1-alpha)*Toub[k] #Toub canceled out mu[k+1] = (TgtY[k+1] - Toub[k+1])/Ak[k+1] PostY[k] = PreY[k] - (ARRk[k]*mu[k]+200) # Process output else: alpha = (1-RI)*alpha2 Toub[k+1] = alpha*(y[k]-Ak[k]*mu[k])+(1-alpha)*Toub[k] #Toub canceled out mu[k+1] = (TgtY[k+1] - Toub[k+1])/Ak[k+1] PostY[k] = PreY[k] - (ARRk[k]*mu[k]+200) # Process output else: if PostY[k]>2950 or PostY[k]<2650: alpha = 0 Toub[k+1] = alpha*(y[k]-Ak[k]*mu[k])+(1-alpha)*Toub[k] #Toub canceled out mu[k+1] = (TgtY[k+1] - Toub[k+1])/Ak[k+1] PostY[k] = PreY[k] - (ARRk[k]*mu[k]+200) # Process output else: alpha = alpha2 Toub[k+1] = alpha*(y[k]-Ak[k]*mu[k])+(1-alpha)*Toub[k] #Toub canceled out mu[k+1] = (TgtY[k+1] - Toub[k+1])/Ak[k+1] PostY[k] = PreY[k] - (ARRk[k]*mu[k]+200) # Process output
In [ ]:
x = range(10, 590) U = np.empty(580) L = np.empty(580) U.fill(2950) L.fill(2650) # plot(x, PostY[3:600], type="o") plt.plot(x,U) plt.plot(x,L) plt.plot(x,PostY[10:590]) np.mean(PostY[10:590])
In [ ]:
x = range(10, 590) plt.plot(x,Toub[10:590])
In [ ]:
Cpk=np.zeros(len(PreY)) MAPEp =np.zeros(len(PreY))
In [ ]:
for k in range(3,599): Cpk[k] = min((UCL-np.mean(PostY[1:k]))/(3*np.std(PostY[1:k], dtype=np.float64)), (np.mean(PostY[1:k])-LCL)/(3*np.std(PostY[1:k], dtype=np.float64))) plt.plot(x,Cpk[10:590])
In [ ]:
for k in range(3,599): MAPEp[k] = sum(np.absolute((PostY[1:k]-2800)/2800))/k*100 plt.plot(x,MAPEp[10:590])
In [ ]:
In [ ]:
In [ ]: