import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
The basic Runto Run Controller
The W2W control scheme:
1) $y_k$ is the actual removal amount measured from the metrology tool and $PostY_k $ is the actual post CMP thickness of run k. The specification of $PostY_k$ is 2800±150 Angstrom (Å) with 2800 being the target value denoted by $TgtPostY$ Deifned $PreY_k$ be the one before CMP process, $ARR_k$ is the polish rate, $\mu_k$ is 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.
2) The $A_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)
3) Process gain:
4) Simulation parameters:
5) Assumption:
yz is atual measurement value, and A is used internal by R2R
from statsmodels.tsa.arima_model import ARIMA
import sklearn.linear_model as skl_lm
# Simulated environment
Error = np.random.normal(0, 25, 600)
PM1 = np.random.normal(0, 50, 600)
PM2 = np.random.normal(0, 3, 600)
Stress1 = np.random.normal(1000, 55,600)
Stress2 = np.random.normal(0, 4, 600)
Rotspd1 = np.random.normal(100, 5.5, 600)
Rotspd2 = np.random.normal(0, 0.4, 600)
Sfuspd1 = np.random.normal(100, 5.5, 600)
Sfuspd2 = np.random.normal(0, 0.4, 600)
PreY = np.random.normal(3800, 70, 600)
TgtPostY = np.array([2800]*600)
# Assumption
Virtual = np.random.normal(1, np.sqrt(0.01), 600)
Ak0 = 1200 # Not sure, does not know PU
Ak = np.array([1200]*600)
RI = 0.9
etak0 = 0
etak = 500
TgtY = PreY - TgtPostY
def ARR (stress1, stress2, rotspd1, rotspd2, sfuspd1, sfuspd2, pm1, pm2, error):
return Ak0*((stress1+stress2)/1000)*((rotspd1+rotspd2)/100)*((sfuspd1+sfuspd2)/100)+pm1+pm2+error
ARRk = ARR(Stress1, Stress2, Rotspd1, Rotspd2, Sfuspd1, Sfuspd2, PM1, PM2, Error) # beta1 actual process gain
RI=RI*np.sqrt(np.sqrt(Virtual))
# ARIMA (0 1 1) assume error
model = ARIMA(Error, order=(0, 1, 1))
results_AR = model.fit(disp=-1)
alpha1 = 0.35
etak = np.zeros(len(PreY))
mu = np.zeros(len(PreY))
PostY = np.zeros(len(PreY))
yz = np.zeros(len(PreY))
mu[0] = (TgtY[0] - etak0)/Ak[0]
yz[0] = ARRk[0]*mu[0]+ results_AR.fittedvalues[0] # Actual value
PostY[0] = PreY[0] - yz[0]
etak[1] = alpha1*(yz[0]-Ak[0]*mu[0])+(1-alpha1)*etak0
mu[1] = (TgtY[1] - etak[1])/Ak[1]
yz[1] = ARRk[1]*mu[1]+ results_AR.fittedvalues[1] # Actual value
PostY[1] = PreY[1] - yz[1]
etak[2] = alpha1*(yz[1]-Ak[1]*mu[1])+(1-alpha1)*etak[1]
mu[2] = (TgtY[2] - etak[2])/Ak[2]
for k in range(2,599):
yz[k] = ARRk[k]*mu[k]+ 700 + results_AR.fittedvalues[k] ## ARRk[k], Toub[k] approximate beta_1, beta_0(not tou!!)
PostY[k] = PreY[k] - yz[k]
etak[k+1] = alpha1*(yz[k]-Ak[k]*mu[k])+(1-alpha1)*etak[k]
mu[k+1] = (TgtY[k+1] - etak[k+1])/Ak[k+1]
# Regression coefficients (Ordinary Least Squares)
mu = mu.reshape((600,1))
yz = yz.reshape((600,1))
regr = skl_lm.LinearRegression()
regr.fit(mu,yz)
print(regr.intercept_)
print(regr.coef_)
print(regr.score(mu,yz))
# Simulated environment
Error = np.random.normal(0, 25, 600)
PM1 = np.random.normal(0, 50, 600)
PM2 = np.random.normal(0, 3, 600)
Stress1 = np.random.normal(1000, 55,600)
Stress2 = np.random.normal(0, 4, 600)
Rotspd1 = np.random.normal(100, 5.5, 600)
Rotspd2 = np.random.normal(0, 0.4, 600)
Sfuspd1 = np.random.normal(100, 5.5, 600)
Sfuspd2 = np.random.normal(0, 0.4, 600)
PreY = np.random.normal(3800, 70, 600)
TgtPostY = np.array([2800]*600)
# Assumption
Virtual = np.random.normal(1, np.sqrt(0.01), 600)
Ak0 = 1200 # Not sure, does not know PU
Ak = np.array([1200]*600)
RI = 0.9
etak0 = 0
etak = 500
TgtY = PreY - TgtPostY
def ARR (stress1, stress2, rotspd1, rotspd2, sfuspd1, sfuspd2, pm1, pm2, error):
return Ak0*((stress1+stress2)/1000)*((rotspd1+rotspd2)/100)*((sfuspd1+sfuspd2)/100)+pm1+pm2+error
ARRk = ARR(Stress1, Stress2, Rotspd1, Rotspd2, Sfuspd1, Sfuspd2, PM1, PM2, Error) # beta1 actual process gain
# ARIMA (0 1 1) assume error
model = ARIMA(Error, order=(0, 1, 1))
results_AR = model.fit(disp=-1)
RI=RI*np.sqrt(np.sqrt(Virtual))
First two rounds, we use actual metrodlogy value
alpha1 = 0.35
etak = np.zeros(len(PreY))
mu = np.zeros(len(PreY))
PostY = np.zeros(len(PreY))
yz = np.zeros(len(PreY))
mu[0] = (TgtY[0] - etak0)/Ak[0]
yz[0] = ARRk[0]*mu[0]+ results_AR.fittedvalues[0] # Actual value
PostY[0] = PreY[0] - yz[0]
etak[1] = alpha1*(yz[0]-Ak[0]*mu[0])+(1-alpha1)*etak0
mu[1] = (TgtY[1] - etak[1])/Ak[1]
print etak[1]
print yz[0]
print mu[0]
print mu[1]
print PostY[0]
yz[1] = ARRk[1]*mu[1]+ results_AR.fittedvalues[1] # Actual value
PostY[1] = PreY[1] - yz[1]
etak[2] = alpha1*(yz[1]-Ak[1]*mu[1])+(1-alpha1)*etak[1]
mu[2] = (TgtY[2] - etak[2])/Ak[2]
print etak[2]
print yz[1]
print mu[2]
print PostY[1]
$𝛼_1$ set to 0.35, and all actual metrology data are available
for k in range(2,599):
yz[k] = ARRk[k]*mu[k]+ results_AR.fittedvalues[k] ## ARRk[k], Toub[k] approximate beta_1, beta_0(not tou!!)
PostY[k] = PreY[k] - yz[k]
etak[k+1] = alpha1*(yz[k]-Ak[k]*mu[k])+(1-alpha1)*etak[k]
mu[k+1] = (TgtY[k+1] - etak[k+1])/Ak[k+1]
for k in range(0, 7):
print etak[k+1]
for k in range(0, 7):
print mu[k+1]
for k in range(0, 7):
print yz[k]
for k in range(0, 7):
print PostY[k]
x = range(2, 590)
plt.plot(x,etak[2:590])
# Plotting
x = range(10, 590)
U = np.empty(580)
L = np.empty(580)
UCL = 2950
LCL = 2650
U.fill(2950)
L.fill(2650)
plt.plot(x,U)
plt.plot(x,L)
plt.plot(x,PostY[10:590])
plt.ylim([2600, 3000])
np.mean(PostY[10:590])
CpK(Process Capability)
Cpk=np.zeros(len(PreY))
MAPEp =np.zeros(len(PreY))
for k in range(2,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])
plt.ylim([0.3, 2.0])
Mean-absolute-percentage error (MAPEp)
for k in range(2,599):
MAPEp[k] = sum(np.absolute((PostY[1:k]-2800)/2800))/k*100
plt.plot(x,MAPEp[10:590])
plt.ylim([0.7, 2.5])
$𝛼_2=𝛼_1=0.35$, apply to following wafers
# Reinitialized
alpha1 = 0.35
betak = 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))
# Round1
mu[0] = (TgtY[0] - betak0)/Ak[0]
yz[0] = ARRk[0]*mu[0]+beta0+results_AR.fittedvalues[0] # Actual value
PostY[0] = PreY[0] - yz[0]
betak[1] = alpha1*(yz[0]-Ak[0]*mu[0])+(1-alpha1)*betak0
mu[1] = (TgtY[1] - betak[1])/Ak[1]
# Round2
yz[1] = ARRk[1]*mu[1]+beta0+results_AR.fittedvalues[1] # Actual value
PostY[1] = PreY[1] - yz[1]
betak[2] = alpha1*(yz[1]-Ak[1]*mu[1])+(1-alpha1)*betak[1]
mu[2] = (TgtY[2] - betak[2])/Ak[2]
for k in range(2,599):
yz[k] = ARRk[k]*mu[k]+results_AR.fittedvalues[k]+beta0 # Metrology data
if(k%25)==1:
y[k] = yz[k]
else:
y[k] = yz[k]*Virtual[k] # Assume VM value ok ARR=ARR_bar
PostY[k] = PreY[k] - yz[k]
betak[k+1] = alpha1*(y[k]-Ak[k]*mu[k])+(1-alpha1)*betak[k]
mu[k+1] = (TgtY[k+1] - betak[k+1])/Ak[k+1]
x = range(2, 590)
plt.plot(x,betak[2:590])
x = range(10, 590)
U = np.empty(580)
L = np.empty(580)
U.fill(2950)
L.fill(2650)
plt.plot(x,U)
plt.plot(x,L)
plt.plot(x,PostY[10:590])
np.mean(PostY[10:590])
plt.ylim([2600, 3000])
Cpk=np.zeros(len(PreY))
MAPEp =np.zeros(len(PreY))
for k in range(2,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])
plt.ylim([0.3, 2.0])
for k in range(2,599):
MAPEp[k] = sum(np.absolute((PostY[1:k]-2800)/2800))/k*100
plt.plot(x,MAPEp[10:590])
plt.ylim([0.7, 2.5])
# Reinitialized
alpha1 = 0.35
betak = 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))
# Round1
mu[0] = (TgtY[0] - betak0)/Ak[0]
yz[0] = ARRk[0]*mu[0]+beta0+results_AR.fittedvalues[0] # Actual value
PostY[0] = PreY[0] - yz[0]
betak[1] = alpha1*(yz[0]-Ak[0]*mu[0])+(1-alpha1)*betak0
mu[1] = (TgtY[1] - betak[1])/Ak[1]
# Round2
yz[1] = ARRk[1]*mu[1]+beta0+results_AR.fittedvalues[1] # Actual value
PostY[1] = PreY[1] - yz[1]
betak[2] = alpha1*(yz[1]-Ak[1]*mu[1])+(1-alpha1)*betak[1]
mu[2] = (TgtY[2] - betak[2])/Ak[2]
for k in range(2,599):
yz[k] = ARRk[k]*mu[k]+beta0+results_AR.fittedvalues[k] # process output how to filter out
if(k%25)==1:
y[k] = yz[k]
alpha = alpha1
else:
y[k] = yz[k]*Virtual[k]
alpha = RI[k]*alpha1
PostYv[k] = PreY[k] - yz[k]*Virtual[k]
if PostYv[k]>2950 or PostYv[k]<2650:
alpha = 0
betak[k+1] = alpha*(y[k]-Ak[k]*mu[k])+(1-alpha)*betak[k] #Toub canceled out
mu[k+1] = (TgtY[k+1] - betak[k+1])/Ak[k+1]
PostY[k] = PreY[k] - yz[k] # Process output
x = range(2, 590)
plt.plot(x,betak[2:590])
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])
plt.ylim([2600, 3000])
Cpk=np.zeros(len(PreY))
MAPEp =np.zeros(len(PreY))
for k in range(2,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])
plt.ylim([0.3, 2.0])
for k in range(2,599):
MAPEp[k] = sum(np.absolute((PostY[1:k]-2800)/2800))/k*100
plt.plot(x,MAPEp[10:590])
plt.ylim([0.7, 2.5])
# Reinitialized
alpha1 = 0.35
betak = 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))
# Round1
mu[0] = (TgtY[0] - betak0)/Ak[0]
yz[0] = ARRk[0]*mu[0]+beta0+results_AR.fittedvalues[0] # Actual value
PostY[0] = PreY[0] - yz[0]
betak[1] = alpha1*(yz[0]-Ak[0]*mu[0])+(1-alpha1)*betak0
mu[1] = (TgtY[1] - betak[1])/Ak[1]
# Round2
yz[1] = ARRk[1]*mu[1]+beta0+results_AR.fittedvalues[1] # Actual value
PostY[1] = PreY[1] - yz[1]
betak[2] = alpha1*(yz[1]-Ak[1]*mu[1])+(1-alpha1)*betak[1]
mu[2] = (TgtY[2] - betak[2])/Ak[2]
for k in range(2,599):
yz[k] = ARRk[k]*mu[k]+beta0+results_AR.fittedvalues[k] # process output how to filter out
if(k%25)==1:
y[k] = yz[k]
alpha = alpha1
else:
y[k] = yz[k]*Virtual[k]
alpha = (1-RI[k])*alpha1
PostYv[k] = PreY[k] - yz[k]*Virtual[k]
if PostYv[k]>2950 or PostYv[k]<2650:
alpha = 0
betak[k+1] = alpha*(y[k]-Ak[k]*mu[k])+(1-alpha)*betak[k] #Toub canceled out
mu[k+1] = (TgtY[k+1] - betak[k+1])/Ak[k+1]
PostY[k] = PreY[k] - yz[k] # Process output
x = range(2, 590)
plt.plot(x,betak[2:590])
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])
plt.ylim([2600, 3000])
Cpk=np.zeros(len(PreY))
MAPEp =np.zeros(len(PreY))
for k in range(2,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])
plt.ylim([0.3, 2.0])
for k in range(2,599):
MAPEp[k] = sum(np.absolute((PostY[1:k]-2800)/2800))/k*100
plt.plot(x,MAPEp[10:590])
plt.ylim([0.7, 2.5])
# Reinitialized
alpha1 = 0.35
betak = 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))
# Round1
mu[0] = (TgtY[0] - betak0)/Ak[0]
yz[0] = ARRk[0]*mu[0]+beta0+results_AR.fittedvalues[0] # Actual value
PostY[0] = PreY[0] - yz[0]
betak[1] = alpha1*(yz[0]-Ak[0]*mu[0])+(1-alpha1)*betak0
mu[1] = (TgtY[1] - betak[1])/Ak[1]
# Round2
yz[1] = ARRk[1]*mu[1]+beta0+results_AR.fittedvalues[1] # Actual value
PostY[1] = PreY[1] - yz[1]
betak[2] = alpha1*(yz[1]-Ak[1]*mu[1])+(1-alpha1)*betak[1]
mu[2] = (TgtY[2] - betak[2])/Ak[2]
for k in range(2,599):
yz[k] = ARRk[k]*mu[k]+beta0+results_AR.fittedvalues[k]
PostYv[k] = PreY[k]- yz[k]*Virtual[k]
if k>= 25:
if PostYv[k]>2950 or PostYv[k]<2650: # should be replace with RI, GSI
alpha = 0
else:
alpha = (1-RI[k])*alpha1
else:
if PostYv[k]>2950 or PostYv[k]<2650:
alpha = 0
else:
alpha = RI[k]*alpha1
if(k%25)==1:
y[k] = yz[k]
alpha = alpha1
else:
y[k] = yz[k]*Virtual[k]
betak[k+1] = alpha*(y[k]-Ak[k]*mu[k])+(1-alpha)*betak[k] #Toub canceled out
mu[k+1] = (TgtY[k+1] - betak[k+1])/Ak[k+1]
PostY[k] = PreY[k] - yz[k] # Process output
x = range(2, 590)
plt.plot(x,betak[2:590])
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])
plt.ylim([2600, 3000])
Cpk=np.zeros(len(PreY))
MAPEp =np.zeros(len(PreY))
for k in range(2,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])
plt.ylim([0.3, 2.0])
for k in range(2,599):
MAPEp[k] = sum(np.absolute((PostY[1:k]-2800)/2800))/k*100
plt.plot(x,MAPEp[10:590])
plt.ylim([0.7, 2.5])