import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

System Prototype¶

The basic Runto Run Controller

The W2W control scheme:

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

CMP Example¶

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:

$A_k$ is set to 300, and A is set to mean of $ARR_k$ (also 300)
$\beta_0$ is set to 700 and $b_0$ is 500
RI is simulated with 0.9 plus some variation, and assume controller can detect any value that is beyond UCL or LCL
yk from AVM is virtual times yz(actual metrodlogy data), where virtual is a normal distribution variable
$\eta_k$ is simulated as ARIMA(0 0 1)
Control $\beta$ not $\eta$

　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))

[ 755.61616784]
[[ 962.33970185]]
0.531339007585

# 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))

Round 1¶

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]

9.4087832029
1064.34565718
0.864552849547
0.872436677123
2773.11776228

Round2¶

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]

45.0780102323
1158.24487298
0.790866989964
2698.08792277

CASE1 R2R with in-situ metrology¶

$𝛼_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]

9.4087832029
45.0780102323
24.864613292
5.37253882486
8.19373790403
38.2825699152
65.735506947

for k in range(0, 7):
    print mu[k+1]

0.872436677123
0.790866989964
0.881409276899
0.822494223586
0.760932003591
0.769812147481
0.858734964209

for k in range(0, 7):
    print yz[k]

1064.34565718
1158.24487298
936.365835502
1026.86410424
1000.42617593
1007.28023367
1040.49410984

for k in range(0, 7):
    print PostY[k]

2773.11776228
2698.08792277
2857.75256269
2855.69164133
2791.9394312
2714.03190854
2721.56303705

x = range(2, 590)
plt.plot(x,etak[2:590])

[<matplotlib.lines.Line2D at 0x7f2d3dbdc690>]

# 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])

2800.5189283424916

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])

(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])

(0.7, 2.5)

CASE2 R2R+VM without RI¶

$𝛼_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])

[<matplotlib.lines.Line2D at 0x7f68a7814c50>]

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])

(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])

(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])

(0.7, 2.5)

CASE3: R2R+VM with RI¶

# 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])

[<matplotlib.lines.Line2D at 0x7f68a74fb850>]

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])

(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])

(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])

(0.7, 2.5)

CASE4: R2R+VM with (1-RI)¶

# 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])

[<matplotlib.lines.Line2D at 0x7f68a77adbd0>]

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])

(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])

(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])

(0.7, 2.5)

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

# 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])

[<matplotlib.lines.Line2D at 0x7f68a76ad550>]

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])

(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])

(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])

(0.7, 2.5)