CoCalc Public FilesSPC_4.htmlOpen in with one click!
Author: phonchi chung
Views : 11
Description: Jupyter html version of SPC_4.ipynb
SPC_4
In [8]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

System Prototype

The basic Runto Run Controller

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 is the polish rate, μk\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 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)

3) Process gain:

4) Simulation parameters:

5) Assumption:

  1. AkA_k is set to 300, and A is set to mean of ARRkARR_k (also 300)
  2. β0\beta_0 is set to 700 and b0b_0 is 500
  3. RI is simulated with 0.9 plus some variation, and assume controller can detect any value that is beyond UCL or LCL
  4. yk from AVM is virtual times yz(actual metrodlogy data), where virtual is a normal distribution variable
  5. ηk\eta_k is simulated as ARIMA(0 1 1)
  6. Control β\beta not η\eta

 yz is atual measurement value, and A is used internal by R2R

In [9]:
from statsmodels.tsa.arima_model import ARIMA
import sklearn.linear_model as skl_lm
In [10]:
# Simulated environment
Error   = np.random.normal(0, 25,    600)
PM1     = np.random.normal(0, 50,    600)
PM2     = np.random.normal(0, 3.5,      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 = 1100 # Not sure, does not know PU
Ak = np.array([1100]*600)
RI = 0.9
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) 
beta0 = 0
RI=RI*np.sqrt(np.sqrt(Virtual))

alpha1 = 0.35
betak = np.zeros(len(PreY))
mu = np.zeros(len(PreY))
PostY = np.zeros(len(PreY))
yz = np.zeros(len(PreY))

mu[0] = (TgtY[0] - beta0)/Ak[0]
yz[0] = ARRk[0]*mu[0]#+ results_AR.fittedvalues[0]  # Actual value
PostY[0] = PreY[0] - yz[0]
betak[1] = alpha1*(yz[0]-Ak[0]*mu[0])+(1-alpha1)*beta0 
mu[1] = (TgtY[1] - betak[1])/Ak[1]

yz[1] = ARRk[1]*mu[1]#+ 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]  ## ARRk[k], Toub[k] approximate beta_1, beta_0(not tou!!)
   PostY[k] = PreY[k] - yz[k]
   betak[k+1] = alpha1*(yz[k]-Ak[k]*mu[k])+(1-alpha1)*betak[k] 
   mu[k+1] = (TgtY[k+1] - betak[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)

# Plotting
x = range(1, 20)
U = np.empty(19)
L = np.empty(19)
UCL = 2950
LCL = 2650
#U.fill(2950)
#L.fill(2650)
plt.plot(x,U)
plt.plot(x,L)
plt.plot(x,ARRk[1:20])
#plt.ylim([2600, 3000])
np.mean(PostY[10:590])

print(regr.intercept_)
print(regr.coef_)
[-17.34643741]
[[ 1115.39468816]]
In [12]:
x = range(1, 20)
U = np.empty(19)
L = np.empty(19)
UCL = 2950
LCL = 2650
U.fill(2950)
L.fill(2650)
plt.plot(x,U)
plt.plot(x,L)
plt.plot(x,PostY[1:20])
Out[12]:
[<matplotlib.lines.Line2D at 0x7fbbc0808910>]
In [4]:
x = range(2, 590)
plt.plot(x,betak[2:590])
Out[4]:
[<matplotlib.lines.Line2D at 0x7fbbc0c4bc50>]
In [ ]:
x = range(2, 590)
plt.plot(x,mu[2:590])
In [ ]:
RI=RI*np.sqrt(np.sqrt(Virtual))
RI

Round 1

First two rounds, we use actual metrodlogy value

In [ ]:
alpha1 = 0.35
betak = np.zeros(len(PreY))
mu = np.zeros(len(PreY))
PostY = np.zeros(len(PreY))
yz = np.zeros(len(PreY))
In [ ]:
mu[0] = (TgtY[0] - beta0)/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)*beta0 
mu[1] = (TgtY[1] - betak[1])/Ak[1]

Round2

In [ ]:
yz[1] = ARRk[1]*mu[1]+ betak[1] + 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]

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]+ betak[k] + results_AR.fittedvalues[k] ## ARRk[k], Toub[k] approximate beta_1, beta_0(not tou!!)
   PostY[k] = PreY[k] - yz[k]
   betak[k+1] = alpha1*(yz[k]-Ak[k]*mu[k])+(1-alpha1)*betak[k] 
   mu[k+1] = (TgtY[k+1] - betak[k+1])/Ak[k+1]
In [ ]:
x = range(2, 590)
plt.plot(x,betak[2:590])
In [ ]:
# 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)

In [ ]:
Cpk=np.zeros(len(PreY))
MAPEp =np.zeros(len(PreY))
In [ ]:
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)

In [ ]:
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])

CASE2 R2R+VM without RI

, apply to following wafers

In [ ]:
# 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]
In [ ]:
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]                 
In [ ]:
x = range(2, 590)
plt.plot(x,betak[2:590])
In [ ]:
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])
In [ ]:
Cpk=np.zeros(len(PreY))
MAPEp =np.zeros(len(PreY))
In [ ]:
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])
In [ ]:
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])

CASE3: R2R+VM with RI

In [ ]:
# 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]
In [ ]:
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 
In [ ]:
x = range(2, 590)
plt.plot(x,betak[2:590])
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])
plt.ylim([2600, 3000])
In [ ]:
Cpk=np.zeros(len(PreY))
MAPEp =np.zeros(len(PreY))
In [ ]:
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])
In [ ]:
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])

CASE4: R2R+VM with (1-RI)

In [ ]:
# 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]
In [ ]:
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 
In [ ]:
x = range(2, 590)
plt.plot(x,betak[2:590])
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])
plt.ylim([2600, 3000])
In [ ]:
Cpk=np.zeros(len(PreY))
MAPEp =np.zeros(len(PreY))
In [ ]:
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])
In [ ]:
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])

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

In [ ]:
# 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]
In [ ]:
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
In [ ]:
x = range(2, 590)
plt.plot(x,betak[2:590])
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])
plt.ylim([2600, 3000])
In [ ]:
Cpk=np.zeros(len(PreY))
MAPEp =np.zeros(len(PreY))
In [ ]:
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])
In [ ]:
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])
In [ ]:
 
In [ ]:
 
In [ ]: