以光刻机工件台为典型代表的超精密运动系统对高速、高加减速工况下的运动精度有极高要求，前馈控制技术是实现超精密运动系统严苛的轨迹跟踪精度目标的重要手段^[1-3]。

前馈控制的基本思路是拟合被控对象的逆模型，其中比较经典的方法是基于模型(model-based)的前馈^[4-6]。然而，光刻机工件台一般都具有结构复杂、模型不确定性大、柔性振动模态凸显等特征，难以获得精确的模型结构及参数，不可避免地会产生建模误差；基于模型的前馈对系统非最小相位零点的近似处理也会产生求逆误差；上述建模、求逆误差的存在导致基于模型的前馈难以完全拟合被控对象模型的逆，最终难以达到光刻机工件台纳米级运动精度的要求^[7]。

相比之下，数据驱动(data-driven)的方法避免了建模、求逆的过程，将系统视作黑箱，通过直接采集系统的输入、输出等控制过程数据，将控制器的参数整定问题转化为关于某一目标函数的优化问题^{[3, 8]}。根据前馈形式的不同，数据驱动前馈可进一步分为基于信号(signal-based)的前馈和基于控制器结构(structure-based)的前馈。

迭代学习控制(iterative learning control, ILC)是一类典型的基于信号的前馈^[9-10]。对于重复性轨迹跟踪问题，ILC通过利用上次实验的误差信号生成下次实验的前馈补偿量，理论上可以完全消除重复性误差，是轨迹跟踪精度最高的控制算法之一。然而ILC严重依赖初始条件，当输入的参考轨迹发生变化时，控制效果可能会严重变差，此时需要重新学习补偿信号^[11]。

基于控制器结构的前馈具有更强的轨迹泛化性能，基本思路与基于模型的前馈一致，但不再利用模型相关信息来逼近逆模型，而是通过优化一组与模型无关的基函数的系数来逼近逆模型^[12-13]。最常用的是多项式类型前馈，如四阶前馈^[14]，采用微分算子作为基函数，将参考轨迹的速度(velocity)、加速度(acceleration)、加速度一阶导数(jerk)、加速度二阶导数(snap)信号的加权组合作为前馈补偿信号，由于具有物理直观性，该方法被广泛应用在光刻机工件台^[15]、打印机^[14]等精密运动系统的控制中。多项式类型前馈的线性参数化结构保证了前馈控制器一定是稳定的，并且控制器参数整定为凸优化问题也易于获取最优参数，但该线性参数化结构也导致只能拟合含有分母的被控对象，难以补偿系统的高频段柔性谐振模态。针对上述问题，文[16-17]提出了有理分式类型前馈。相比于多项式的有限冲击响应结构，有理分式的无限冲击响应结构使有理分式类型前馈具备更强的逆模型拟合能力，但同时也导致前馈控制器参数整定为非凸问题，一般的参数优化算法难以保证获取最优值^[12]，并且在迭代过程中还需要考虑前馈控制器分母的稳定性问题。

综上所述，无论是基于模型的前馈，还是基于控制器结构的数据驱动前馈，拟合逆模型的精度受限于前馈选择的形式，最终会不可避免地产生拟合残差；在关注的频段范围内，该拟合残差会严重影响系统最终的运动精度。因此需要在上述前馈基础上添加前馈补偿项，进一步提高轨迹跟踪精度^[18-19]。由于ILC训练得到的前馈控制信号中包含丰富的模型信息，有学者提出了通过拟合ILC信号来整定补偿前馈控制器参数的方法^{[18, 20-21]}。

为减少逆模型的拟合残差对光刻机工件台控制性能的影响，本文在现有四阶前馈的基础上，提出一种新的通过数据驱动方法迭代整定参数的前馈补偿结构。通过优化广义残余重复误差，将有理分式结构补偿前馈参数整定过程中的非凸优化问题转化为凸优化问题，进而基于Gauss-Newton法给出了补偿前馈全局最优参数更新率，同时给出了参数整定过程中梯度、Hessian矩阵及有关变量的无偏估计方法。

1 光刻机工件台及其控制方法 1.1 光刻机工件台

在集成电路的制造过程中，光刻机工件台用于承载晶圆，完成测量、曝光等工艺流程。清华大学机械工程系IC装备实验室开发的一款光刻机工件台如图 1所示，采用粗精叠层结构：粗动台由直线电机驱动，该直线电机的磁钢固定在大理石台上，动子线圈固定在粗动台上，通过直线导轨使粗动台在Y向进行大行程微米级运动；微动台叠加在粗动台上，由4个水平音圈电机和4个垂向音圈电机产生六自由度运动，通过固定在机架上的亚纳米分辨率激光干涉仪反馈实时位移，实现小行程纳米级运动。光刻机高产率(≥295硅片/h)、高套刻精度(≤2.5 nm)的整机指标需要工件台在高速、高加速的运动条件下实现纳米级的运动精度，为控制带来极大的挑战。

微动台运动系统通常可由多质量块模型建模，其传递函数一般可表示为^[18-19]

$ P(s)=\frac{1}{M s^{2}}+\frac{1}{M} \sum\limits_{i=1}^{I} \frac{\gamma_{i}}{s^{2}+2 \eta_{i} \omega_{i} s+\omega_{i}^{2}} . $

(1)

其中: $P(s)$表示$P$是关于$s$的函数, 为被控对象, $s$为Laplace变换中的复变数, 为书写简单, 后文会省略传递函数中的$s ; M$为系统总质量; $\eta_{i} 、\omega_{i} 、\gamma_{i}$分别为第$i$个谐振的阻尼比、谐振频率和常数系数, $i=1, 2, \cdots, I$。在低频段精密运动系统主要呈现双积分特性, 而在高频段柔性谐振模态开始凸显。

1.2 工件台控制方法

工件台运动系统通常采用“前馈+反馈”的二自由度控制结构^{[2, 8]}，如图 2所示。

图 2中：$r$为参考轨迹; $y$为系统的实际输出轨迹; $e=r-y$为轨迹跟踪误差(简称“误差”); v为测量噪声; $u$为$P(s)$的控制力信号; $u_{\mathrm{ff}}$为前馈(feedforward) 控制力, 对应前馈控制器$F(s)$的输出; $u_{\mathrm{fb}}$为反馈(feedback) 控制力, 对应反馈控制器$C(s)$的输出; $u_{\mathrm{c}}$为额外补偿(compensation) 力信号, 易知有$u=u_{\mathrm{fb}}+u_{\mathrm{ff}}+u_{\mathrm{c}}$。

当不考虑$u_{\mathrm{c}}$时, $e$与$r 、v$之间的传递函数关系可表示为

$ e=\underbrace{S(s) r-S_{\mathrm{p}}(s) F r}_{e_{F}}-\underbrace{S(s) v}_{e_{v}} . $

(2)

其中: $S(s)$为灵敏度函数, $S(s)=1 /(1+P(s) C(s))$; $S_{\mathrm{p}}(s)$为过程(process) 灵敏度函数, $S_{\mathrm{p}}(s)=S(s) P(s)$; $e_{F}$为使用前馈后的轨迹导致的重复性误差分量; $e_{v}$为由噪声导致的非重复性误差分量, $e_{v}=S(s) v$。$C(s)$在$S(s)$的分母上, 主要用来保证闭环系统稳定, 同时减少噪声干扰; $F(s)$利用轨迹信息用于消除$e_{F}$。显然, 当$F(s)$等于$P(s)$的逆, 即

$ F(s)=P^{-1}(s) \text {. } $

(3)

此时, $e_{F}=0$。

由于光刻机工件台微动台的无阻尼特性，通常$F(s)$采用$r$的加速度、加速度一阶导数、加速度二阶导数信号加权组合的形式，也即四阶前馈^[14]，此时$F(s)$可表示为

$ F(s)=\varepsilon_{\text {acc }} s^{2}+\varepsilon_{\text {jerk }} s^{3}+\varepsilon_{\text {snap }} s^{4} . $

(4)

其中: $\varepsilon_{\text {acc }}$为加速度权重系数, $\varepsilon_{\text {jerk }}$为加速度一阶导数权重系数, $\varepsilon_{\text {snap }}$为加速度二阶导数权重系数。文[18]指出, 式(1) 所示$P(s)$的逆通过Taylor展开可表示为

$ P^{-1}(s)=\mu_{1} s^{2}+\mu_{2} s^{3}+\mu_{3} s^{4}+\Delta G(s) . $

(5)

其中: $\mu_{1} 、\mu_{2} 、\mu_{3}$为Taylor展开的前3项系数, $\Delta G(s)$为5次及以上高阶展开项。如果$r$的能量主要集中在低频段, 那么此时可以忽略$\Delta G(s)$; 当$\varepsilon_{\text {acc }}=\mu_{1} 、\varepsilon_{\text {jerk }}=\mu_{2} 、\varepsilon_{\text {snap }}=\mu_{3}$时, 式(4) 中$F(s)$可以很好地拟合$P(s)$的逆, 这也是四阶前馈被广泛使用的原因。此时前馈控制器$F=P^{-1}-\Delta G, e_{F}$可表示为

$ e_{F}=S_{\mathrm{p}} \cdot \Delta G \cdot r . $

(6)

可以看出, $r$的高阶导数$\Delta G \cdot r$会通过$S_{\mathrm{p}}(s)$最终作用在$e_{F}$上。对于图 1所示的光刻机工件台, 由$C(s)$将$S_{\mathrm{p}}(s)$低频段的幅值特性调整为$-140 \mathrm{~dB}$左右, 在准确拟合四阶前馈参数的情况下，如果$\Delta G \cdot r$的存在导致控制力出现$0.1 \mathrm{~N}$的偏差, 则会产生约$10 \mathrm{~nm}$的重复性误差, 导致最终难以实现纳米级运动精度。因此需要在基本四阶前馈的基础上进一步设计补偿力信号以提高最终的运动精度。

1.3 额外补偿前馈

额外补偿力信号可通过以下2种方式获取。

1) 基于信号的补偿控制。

通过ILC获取补偿力信号, 即$u_{\mathrm{c}}=u_{\mathrm{ILC}}$, 而迭代学习信号$u_{\mathrm{ILC}}$可通过以下迭代进行表示^[10]：

$ u_{\mathrm{ILC}}^{k+1}=L_{1}\left(u_{\mathrm{ILC}}^{k}+L_{2} e^{k}\right) \text {. } $

(7)

其中：$k$表示迭代次数; $L_{1}$为低通滤波器, $L_{2}$为迭代学习控制器

2) 基于控制器的补偿控制。

通过额外补偿前馈控制器$\Delta F$生成补偿力信号，

$ u_{\mathrm{c}}=\Delta F \cdot r^{p} . $

(8)

其中: $r^{p}$表示$r$的$p$阶导数。$\Delta F$选择离散域有理分式结构, 可表示为

$ \begin{gathered} \Delta F=z^{d} \cdot \frac{B\left(z^{-1}, \boldsymbol{\theta}_{B}\right)}{A\left(z^{-1}, \boldsymbol{\theta}_{A}\right)}= \\ z^{d} \cdot \frac{b_{0}+b_{1} z^{-1}+\cdots+b_{m} z^{-m}}{1+a_{1} z^{-1}+\cdots+a_{n} z^{-n}} . \end{gathered} $

(9)

其中: $z$表示时移算子, $z^{d}$表示将$r^{p}$向前时移$d$个采样周期, $z^{-1}$表示向后时移1个采样周期; $A\left(z^{-1}, \boldsymbol{\theta}_{A}\right)$表示$\Delta F$的分母多项式, 其变量为$z^{-1}$, 多项式系数为向量$\boldsymbol{\theta}_{A}$, 阶次为$n$; 有$A\left(z^{-1}\right.$, $\left.\boldsymbol{\theta}_{A}\right)=1+a_{1} z^{-1}+\cdots+a_{n} z^{-n}, \boldsymbol{\theta}_{A}=\left[a_{1}, a_{2}, \cdots\right.$, $\left.a_{n}\right]^{\mathrm{T}} \in \bf{R}^{n}$。同理, $B\left(z^{-1}, \boldsymbol{\theta}_{B}\right)$表示$\Delta F$的分子多项式, 其变量为$z^{-1}$, 多项式系数为向量$\boldsymbol{\theta}_{B}$, 阶次为$m$; 有$B\left(z^{-1}, \boldsymbol{\theta}_{B}\right)=b_{0}+b_{1} z^{-1}+\cdots+b_{m} z^{-m}$, $\boldsymbol{\theta}_{B}=\left[b_{0}, b_{1}, \cdots, b_{m}\right]^{\mathrm{T}} \in \bf{R}^{m+1}$。可知$\Delta F$的待整定参数列向量为$\boldsymbol{\theta}=\left[\boldsymbol{\theta}_{A}^{\mathrm{T}}, \boldsymbol{\theta}_{B}^{\mathrm{T}}\right]^{\mathrm{T}}$, 此时控制结构框图如图 3所示。

由于基于信号的补偿前馈的变轨迹适应能力较差，故而通常采用基于控制器结构的补偿前馈。此时，问题的核心变为如何整定$\Delta F$的参数$\boldsymbol{\theta}$, 使其能够补偿式(6)中由于基本四阶前馈逆模型拟合能力有限导致的重复性误差分量。

2 数据驱动额外补偿前馈参数整定方法

假设$F$已经调节完毕, 由图 3可知, 此时$e$关于向量$\boldsymbol{\theta}$的关系式可表示为

$ e(\boldsymbol{\theta})=\underbrace{S r-S_{\mathrm{p}} F r}_{e_{F}}-\frac{B\left(z^{-1}, \boldsymbol{\theta}_{B}\right)}{A\left(z^{-1}, \boldsymbol{\theta}_{A}\right)} S_{\mathrm{p}} \cdot z^{d} \cdot r^{p}-S v $

(10)

由式(10)知，由于$\boldsymbol{\theta}_{A}$在分母部分, $e$关于$\boldsymbol{\theta}$是非线性的。当采用$e$的加权2范数的平方作为目标函数进行参数整定时，该问题为非凸优化问题，难以实现全局优化。

2.1 全局最优参数整定方法

借鉴文[12]和系统辨识的相关方法，可将上述非凸优化问题转化为凸优化问题。设计ΔF的目的是补偿基本四阶前馈没有消除的残余重复性误差，使下式成立：

$ e_{F}-\frac{B\left(z^{-1}, \boldsymbol{\theta}_{B}\right)}{A\left(z^{-1}, \boldsymbol{\theta}_{A}\right)} S_{\mathrm{p}} z^{d} r^{p}=0 . $

(11)

定义广义残余重复误差$\eta_{r}(\boldsymbol{\theta})$，具体可表示为

$ \eta_{r}(\boldsymbol{\theta})=A\left(z^{-1}, \boldsymbol{\theta}_{A}\right) e_{F}-B\left(z^{-1}, \boldsymbol{\theta}_{B}\right) S_{\mathrm{p}} z^{d} r^{p} . $

(12)

其中$\eta_{r}(\boldsymbol{\theta})$表示关于$\boldsymbol{\theta}$的函数。定义基函数行向量$\boldsymbol{\phi}_{n-1}(z)=\left[1, z^{-1}, \cdots, z^{-(n-1)}\right]$和$\boldsymbol{\phi}_{m}(z)=[1$, $\left.z^{-1}, \cdots, z^{-m}\right]$。式(12)可进一步表示为

$ \begin{gathered} \eta_{r}(\boldsymbol{\theta})=e_{F}+z^{-1} \boldsymbol{\phi}_{n-1}(z) e_{F} \boldsymbol{\theta}_{A}- \\ \boldsymbol{\phi}_{m}(z) S_{\mathrm{p}} z^{d} r^{p} \boldsymbol{\theta}_{B}=e_{F}-\mathit{\boldsymbol{ \boldsymbol{\varPhi} }} \boldsymbol{\theta} . \end{gathered} $

(13)

其中$ \mathit{\boldsymbol{ \boldsymbol{\varPhi} }} $为与$\Delta F$有关的向量, 其构成可表示为

$ \left\{\begin{array}{l} \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}=\left[\begin{array}{ll} -\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{A} & \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{B} \end{array}\right], \\ \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{A}=z^{-1} \boldsymbol{\phi}_{n-1}(z) e_{F}, \\ \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{B}=\boldsymbol{\phi}_{m}(z) S_{\mathrm{p}} z^{d} r^{p} . \end{array}\right. $

(14)

其中: $\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{A} 、\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{B}$分别为与$\Delta F$分母、分子有关的向量。显然, 式(11) 成立等价于令式(13) 等于零, 且式(13) 中$\eta_{r}(\boldsymbol{\theta})$关于$\boldsymbol{\theta}$是线性的。选择$\eta_{r}(\boldsymbol{\theta})$的加权2范数的平方作为参数优化的目标函数$J(\boldsymbol{\theta})$, 具体可表示为

$ J(\boldsymbol{\theta})=\frac{1}{2}\left\|\eta_{r}(\boldsymbol{\theta})\right\|_{Q}^{2}=\frac{1}{2}\left(\eta_{r}(\boldsymbol{\theta})\right)^{\mathrm{T}} \boldsymbol{Q}_{\eta_{r}}(\boldsymbol{\theta}) . $

(15)

其中Q为正定的加权矩阵。此时参数整定问题变为凸优化问题，采用经典的Gauss-Newton优化方法即可获得较优的参数整定效果，参数迭代率可表示为

$ \boldsymbol{\theta}^{k+1}=\boldsymbol{\theta}^{k}-\alpha^{k}\left(\nabla^{2} J\left(\boldsymbol{\theta}^{k}\right)\right)^{-1} \nabla J\left(\boldsymbol{\theta}^{k}\right) . $

(16)

其中: $\alpha^{k}$为第$k$次迭代的步长; $\nabla J\left(\boldsymbol{\theta}^{k}\right) 、\nabla^{2} J\left(\boldsymbol{\theta}^{k}\right)$分别为第$k$次迭代时的梯度和Hessian矩阵, 具体表示如下:

$ \left\{\begin{array}{l} \nabla J\left(\boldsymbol{\theta}^{k}\right)=(-\mathit{\boldsymbol{ \boldsymbol{\varPhi} }})^{\mathrm{T}} \boldsymbol{Q}_{\eta_{r}}\left(\boldsymbol{\theta}^{k}\right), \\ \nabla^{2} J\left(\boldsymbol{\theta}^{k}\right)=\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}^{\mathrm{T}} \boldsymbol{Q} \mathit{\boldsymbol{ \boldsymbol{\varPhi} }} . \end{array}\right. $

(17)

2.2 梯度及Hessian矩阵的无偏估计

采用式(16) 所示算法整定参数时, 除了预先选定$\alpha^{k} 、Q$外, 还需要获取$\eta_{r}\left(\boldsymbol{\theta}^{k}\right) 、\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}$这2个中间变量的值, 其中$\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}$的获取如式(14) 所示, 需要知道$e_{F} 、S_{\mathrm{p}}$的信息。有关变量的估计方法如下:

1) $e_{F}$可直接通过$F(s)$整定后的$e$来估计, 由式(2) 可得其估计值$\hat{e}_{F}$为

$ \hat{e}_{F}=S r-S_{\mathrm{p}} F r-S v . $

(18)

其中^表示对应变量的估计值。可知当v为白噪声时，$\hat{e}_{F}$为$e_{F}$的无偏估计。

2) 由于$A\left(z^{-1}, \boldsymbol{\theta}_{A}^{k}\right)$已知, $\eta_{r}\left(\boldsymbol{\theta}^{k}\right)$可直接通过式(10) 测量第$k$次迭代的误差$e\left(\boldsymbol{\theta}^{k}\right)$来估计, 其估计值为

$ \begin{gathered} \hat{\eta}_{r}\left(\boldsymbol{\theta}^{k}\right)=A\left(z^{-1}, \boldsymbol{\theta}_{A}^{k}\right) e\left(\boldsymbol{\theta}^{k}\right)= \\ \eta_{r}\left(\boldsymbol{\theta}^{k}\right)-A\left(z^{-1}, \boldsymbol{\theta}_{A}^{k}\right) S v . \end{gathered} $

(19)

显然, $\hat{\eta}_{r}\left(\boldsymbol{\theta}^{k}\right)$也为$\eta_{r}\left(\boldsymbol{\theta}^{k}\right)$的无偏估计。

3) $S_{p}$的估计。

实际计算中需要将离散域形式转化为代数形式，此时离散传递函数变为其单位脉冲响应信号对应的Toeplitz矩阵^[1]。后文以相同符号的粗体表示传递函数对应的Toeplitz矩阵。

采用图 4所示脉冲响应(impulse response)实验即可无偏估计$S_{\mathrm{p}}$对应的Toeplitz矩阵$\boldsymbol{S}_{\mathrm{p}}$ ^[1]。

给定参考轨迹$r=0$; 单位脉冲输人为$u_{\mathrm{imp}}$, 其对应的离散采样向量为$\boldsymbol{u}_{\mathrm{imp}}=[1, 0, \cdots, 0]^{\mathrm{T}}$; 采集控制系统的输出为$y$, 得到其离散采样向量$\boldsymbol{y}=[y$ $(0), y(1), \cdots, y(N-1)]^{\mathrm{T}}$, 其中$N$为离散采样点数。因此, 有

$ \underbrace{\left[\begin{array}{c} y(0) \\ y(1) \\ \vdots \\ y(N-1) \end{array}\right]}_{\boldsymbol{y}}=\boldsymbol{S}_{\mathrm{p}} \underbrace{\left[\begin{array}{c} 1 \\ 0 \\ \vdots \\ 0 \end{array}\right]}_{\boldsymbol{u}_{\text {imp }}}+\boldsymbol{S}\left[\begin{array}{c} v(0) \\ v(1) \\ \vdots \\ v(N-1) \end{array}\right] . $

(20)

其中$\boldsymbol{S}$表示灵敏度函数$S$对应的Toeplitz矩阵。利用$\boldsymbol{y}$即可得到$\boldsymbol{S}_{\mathrm{p}}$的无偏估计, $\hat{\boldsymbol{S}}_{\mathrm{p}}$可表示为

$ \begin{aligned} & \underbrace{\left[\begin{array}{cccc} y(0) & 0 & \cdots & 0 \\ y(1) & y(0) & \cdots & 0 \\ \vdots & \vdots & & \vdots \\ y(N-1) & y(N-2) & \cdots & y(0) \end{array}\right]}_{\hat{s}_{\mathrm{p}}}=\boldsymbol{S}_{\mathrm{p}}+ \\ & \boldsymbol{S}\underbrace{\left[\begin{array}{cccc} v(0) & 0 & \cdots & 0 \\ \mid v(1) & v(0) & \cdots & 0 \\ \vdots & \vdots & & \vdots \\ v(N-1) & v(N-2) & \cdots & v(0) \end{array}\right]}_{\boldsymbol{v}} . \end{aligned} $

(21)

其中$\boldsymbol{V}$为$v$对应的Toeplitz矩阵。显然, 在白噪声假设下, $E\left\{\hat{\boldsymbol{S}}_{\mathrm{p}}\right\}=E\left\{\boldsymbol{S}_{\mathrm{p}}+\boldsymbol{S} \boldsymbol{V}\right\}=\boldsymbol{S}_{\mathrm{p}}$, 表明$\hat{\boldsymbol{S}}_{\mathrm{p}}$是无偏的。

获取上述3个中间变量的无偏估计之后，得出如下结论：

定理   通过2次基本四阶前馈轨迹跟踪实验和2次脉冲响应实验, 即可获得$\nabla J\left(\boldsymbol{\theta}^{k}\right) 、\nabla^{2} J\left(\boldsymbol{\theta}^{k}\right)$的无偏估计。

证明   2次基本四阶前馈轨迹跟踪实验$\mathrm{H} 1$、$\mathrm{H} 2$的有关变量使用下标$\mathrm{H} 1 、\mathrm{H} 2$进行标注, 2次脉冲响应实验H3、H4的有关变量使用下标H3、$\mathrm{H} 4$进行标注。假设每次实验的噪声是独立同分布的白噪声，由式(18)和(21) 可以得到

$ \left\{\begin{array}{l} \hat{e}_{F, \mathrm{H} 1}=e_{\mathrm{H} 1}=e_{F}-S v_{\mathrm{H} 1}, \\ \hat{e}_{F, \mathrm{H} 2}=e_{\mathrm{H} 2}=e_{F}-S v_{\mathrm{H} 2}, \\ \hat{\boldsymbol{S}}_{\mathrm{p}, \mathrm{H} 3}=\boldsymbol{S}_{\mathrm{p}}+\boldsymbol{S} \boldsymbol{V}_{\mathrm{H} 3}, \\ \hat{\boldsymbol{S}}_{\mathrm{p}, \mathrm{H} 4}=\boldsymbol{S}_{\mathrm{p}}+\boldsymbol{S} \boldsymbol{V}_{\mathrm{H} 4} . \end{array}\right. $

(22)

将式(22)代人式(14) 中, 分别使用实验H1和H3与实验$\mathrm{H} 2$和$\mathrm{H} 4$计算$\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}$的2个估计$ \mathit{\boldsymbol{ \boldsymbol {\hat\varPhi} }}_{1}$和$ \mathit{\boldsymbol{ \boldsymbol {\hat\varPhi} }}_{2}$, 具体可表示为

$ \left\{\begin{array}{l} \mathit{\boldsymbol{ \boldsymbol {\hat\varPhi} }}_{1}=\left[-\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{A, \mathrm{H} 1} \quad \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{B, \mathrm{H} 3}\right]=\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}+ \\ \underbrace{\left[-z^{-1} \boldsymbol{\phi}_{n-1}(z) S v_{\mathrm{H} 1} \boldsymbol{\phi}_{m}(z) \boldsymbol{S} \boldsymbol{V}_{\mathrm{H} 3} z^{d} r^{p}\right]}_{\Delta \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{1}}, \\ \mathit{\boldsymbol{ \boldsymbol {\hat\varPhi} }}_{2}=\left[\begin{array}{ll} -\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{A}, \mathrm{H} 2 & \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{B}, \mathrm{H} 4 \end{array}\right]=\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}+ \\ \underbrace{\left[-z^{-1} \boldsymbol{\phi}_{n-1}(z) S v_{\mathrm{H} 2} \quad \boldsymbol{\phi}_{m}(z) \boldsymbol{S} \boldsymbol{V}_{\mathrm{H} 4} z^{d} r^{p}\right]}_{\Delta \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{2}} . \end{array}\right. $

(23)

易知式(23) 中的2个估计均是关于$\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}$的无偏估计。令Hessian矩阵的估计为$\hat{\nabla}^{2} J\left(\boldsymbol{\theta}^{k}\right)=\hat{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}}_{1}^{\mathrm{T}} \boldsymbol{Q} \hat{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}}_{2}$, 则有

$ \begin{gathered} E\left\{\hat{\nabla}^{2} J\left(\boldsymbol{\theta}^{k}\right)\right\}=E\left\{\hat{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}}_{1}^{\mathrm{T}} \boldsymbol{Q} \hat{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}}_{2}\right\}= \\ E\left\{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}^{\mathrm{T}} \boldsymbol{Q} \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}+\left(\Delta \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{1}\right)^{\mathrm{T}} \boldsymbol{Q} \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}+\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}^{\mathrm{T}} \boldsymbol{Q} \Delta \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{2}+\right. \\ \left.\Delta \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{1}^{\mathrm{T}} \boldsymbol{Q} \Delta \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{2}\right\}=\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}^{\mathrm{T}} \boldsymbol{Q} \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}=\nabla^{2} J\left(\boldsymbol{\theta}^{k}\right) . \end{gathered} $

(24)

式(24) 表明所提方法得到的Hessian矩阵是无偏的。同样, 令梯度的估计为$\hat{\nabla} J\left(\boldsymbol{\theta}^{k}\right)=-\hat{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}}_{1}^{\mathrm{T}} \boldsymbol{Q} \hat\eta_{r}\left(\boldsymbol{\theta}^{k}\right)$, 则有

$ \begin{gathered} E\left\{\hat{\nabla} J\left(\boldsymbol{\theta}^k\right)\right\}=E\left\{-\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_1^{\mathrm{T}} \boldsymbol{Q} \hat{\eta}_r\left(\boldsymbol{\theta}^k\right)\right\}=-\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}^{\mathrm{T}} \boldsymbol{Q} \eta_r\left(\boldsymbol{\theta}^k\right)- \\ E\left\{\Delta \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_1^{\mathrm{T}} \boldsymbol{Q} \eta_r\left(\boldsymbol{\theta}^k\right)-\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}^{\mathrm{T}} \boldsymbol{Q} A\left(z^{-1}, \boldsymbol{\theta}_A^k\right) S v-\right.\\ \left.\Delta \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_1^{\mathrm{T}} \boldsymbol{Q} A\left(z^{-1}, \boldsymbol{\theta}_A^k\right) S v\right\} . \end{gathered} $

(25)

由于$\hat{\eta}_{r}\left(\boldsymbol{\theta}^{k}\right)$中的噪声$v$是整定额外补偿前馈控制器$\Delta F$时的实验噪声, 与$\mathrm{H} 1 、\mathrm{H} 2 、\mathrm{H} 3 、\mathrm{H} 4$实验中的噪声不相关, 因此$E\left\{\hat{\nabla} J\left(\boldsymbol{\theta}^k\right)\right\}=\nabla J\left(\boldsymbol{\theta}^k\right)$, 表明梯度的估计也是无偏的。

2.3 参数收敛性能及额外补偿前馈控制器的稳定性

理想的最优(optimal) 补偿前馈参数$\boldsymbol{\theta}_{\mathrm{opt}}$能够完全消除$e_{F}$, 即使式(13) 等于零。$\eta_{r}(\boldsymbol{\theta})$可表示为

$ \eta_{r}(\boldsymbol{\theta})=\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}\left(\boldsymbol{\theta}_{\mathrm{opt}}-\boldsymbol{\theta}\right) . $

(26)

代入式(16)和(17)中，得到的参数迭代率可表示为

$ \boldsymbol{\theta}^{k+1}-\boldsymbol{\theta}_{\mathrm{opt}}=\left(1-\alpha^{k}\right)\left(\boldsymbol{\theta}^{k}-\boldsymbol{\theta}_{\mathrm{opt}}\right) \text {. } $

(27)

显然，当步长满足$0<\alpha^{k}<2$时，上述迭代过程是收敛的。当考虑2.2节中的无偏估计方法引入噪声时，可通过类似分析证明$\boldsymbol{\theta}$的期望是收敛的^[22]。

此外, 在参数整定过程中, 有理分式结构中分母$A\left(z^{-1}, \boldsymbol{\theta}_{A}\right)$的参数选取不当可能导致额外补偿前馈控制器$\Delta F$不稳定。对于离散域形式的$\Delta F$, 可通过分析极点与单位圆的位置关系判断$\Delta F$的稳定性; 如果不稳定, 采用稳定反转(stableinversion）的技术手段获取有界的额外补偿力信号, 具体实施可参考文[23]。

3 实验结果与分析 3.1 实验设置

本文以图 1中光刻机工件台的微动台为实验对象。由于预先解耦降低了各自由度间的影响，其Y向可视为单自由度问题，因此本文只关注Y向的轨迹跟踪误差。

反馈控制器由经典的比例-积分-微分(proportional-integral-derivative, PID)控制器、陷波滤波器构成^[8]。通过调节反馈控制器，使工件台定位误差为$\pm 7 \mathrm{~nm}$。在此基础上, 进行轨迹跟踪实验, 基本前馈$F$选用四阶前馈, 参数可通过文[18]中所提方法进行整定。在上述前馈基础上，对比以下2种补偿前馈$\Delta F$的$\boldsymbol{\theta}$整定方法。

1) M1:   基于数据驱动的$\Delta F$的$\boldsymbol{\theta}$整定方法。

2) M2:  文[18]方法，首先进行迭代学习补偿实验，通过式(7)获得$u_{\mathrm{ILC}}$; 然后将$\Delta F$的$\boldsymbol{\theta}$整定问题视为一个输人为$r^{p}$、输出为$u_{\mathrm{ILC}}$的系统辨识问题, 通过最小二乘法求解$\Delta F$的$\boldsymbol{\theta}$。

关于何时停止迭代, 可通过定义收敛准则进行判定, 可参考文[19]; 但是对于ILC或凸优化问题, 收玫速度很快, 在本文实验中, $F 、\Delta F$和ILC都选择5次迭代。

3.2 轨迹跟踪实验结果

光刻机工件台的参考轨迹常用四阶轨迹^{[1, 12]}，本文实验中所使用的四阶轨迹$R_{1}$和$R_{2}$的参数信息如表 1所示。其中, 轨迹$R_{1}$用来训练前馈控制器参数, 轨迹$R_{2}$用来检验该前馈控制器的变轨迹泛化能力。

本文实验中, 选择$\Delta F$的分子阶次等于分母阶次, 即$n=m$, 超前时移量$d=2$。针对M1、M2这2种方法, 分别整定阶次$n=5 、8 、10$的$\Delta F$。

$\boldsymbol{\theta}$整定过程中误差2范数的平方随迭代次数的变化曲线如图 5所示。初始无前馈状态经过5次迭代, 基本四阶前馈已经收玫, 此时误差2范数的平方在$2 \times 10^{5} \mathrm{~nm}^{2}$量级，加速段(加速度大于0)、减速段(加速度小于0)误差幅值为±35 nm，刚进入匀速段(加速度等于0)的误差幅值高达26 nm，导致匀速段的误差振荡较大，振动最大幅值为18 nm，如图 6所示。在基本四阶前馈的框架下，误差难以进一步减小，不能满足高端光刻机纳米级运动精度需求，因此需要通过添加$\Delta F$以进一步提高运动精度。

在基本四阶前馈的基础上，分别添加的ILC补偿、M1方法补偿，结果如图 5和6所示，可以看出：

1) 图 5中，M1方法3种阶次的$\Delta F$的误差2范数的平方几乎收敛到相同值，证明所提参数整定算法具有收敛性。

2) 图 5中，使用M1方法进行补偿时，误差2范数的平方由2×10⁵下降至2×10⁴，接近ILC补偿的水平；图 6中，添加$\Delta F(n=10)$后，加减速段的误差大幅降低，匀速段的误差减小至±10 nm，轨迹跟踪误差接近添加ILC补偿后的误差。上述结果说明添加补偿前馈$\Delta F$有效。

ILC得到的补偿信号通过M2方法同样可以得到3种阶次的$\Delta F$。M1、M2具体轨迹跟踪控制效果对比如图 7所示，其中每个控制器重复进行3次实验。可以看出：

1) 无论是对比误差2范数的平方, 还是对比误差的幅值, M1整定得到的$\Delta F$的轨迹跟踪控制效果(蓝线)优于M2 (黑线)。

2) 对M1而言，图 7c和7d中轨迹跟踪误差的幅值与定位误差基本相当，可认为其主要呈现噪声特性，残余重复性误差分量基本被消除。

3) 对M2而言，当整定阶次较低时，通过M2整定得到的$\Delta F$的控制效果较差，如图 7b和7c所示，误差中仍然存在明显的低频分量，这是由于选择的前馈控制器阶次较低，拟合$u_{\text {ILC }}$不准确导致的；随着阶次升高，其轨迹跟踪精度也越来越高，逐渐接近M1的控制效果。

因此，为了实现良好的补偿效果，M2通常会给出一个阶次较高的前馈控制器。但是对于工件台运动系统而言，其伺服周期较短，留给计算控制补偿量的时间有限，而阶次较高的前馈控制器会增加计算时间。

此外，匀速段是工件台用于曝光的阶段，也是其性能被严格要求的阶段。工件台具体的运动性能指标为移动平均值(moving average, MA)和移动标准差(moving standard deviation, MSD)^[8]，具体可表示为

$ \left\{\begin{array}{l} M_{\mathrm{a}}(k)=\frac{1}{l} \sum\limits_{i=k-l / 2}^{k+l / 2-1} e(i), \\ M_{\mathrm{sd}}(k)=\sqrt{\frac{1}{l} \sum\limits_{i=k-l / 2}^{k+l / 2-1}\left(e(i)-M_{\mathrm{a}}(k)\right)^{2}} . \end{array}\right. $

(28)

其中: $M_{\mathrm{a}}(k) 、M_{\mathrm{sd}}(k)$分别表示第$k$个采样时刻的$\mathrm{MA}$和MSD值; $l$为曝光狭缝宽度; $e(i)$为曝光狭缝内的第$i$个采样误差, $i=k-l / 2, \cdots, k+l / 2-$ 1。由式(28)可知, $M_{\mathrm{a}}(k)$和$M_{\mathrm{sd}}(k)$等于$l$内误差的均值和标准差, 因此$M_{\mathrm{a}}(k) 、M_{\mathrm{sd}}(k)$可分别衡量误差的低频和高频特性。表 2给出了M1、M2这2种方法匀速段误差的$M_{\mathrm{a}} 、M_{\mathrm{sd}}$最大值。

由表 2可知，M1的控制效果明显优于M2。此外，随着$\Delta F$阶次升高，M1匀速段误差的$\left\|M_{\mathrm{a}}\right\|_{\infty}$的平均值不但没有降低，反而略有升高，但是都在3.8~4.0 nm，这可能是因为接近定位精度上限使误差难以进一步降低。

将$R_{1}$下得到的$u_{\mathrm{ILC}} 、\Delta F$应用在$R_{2}$的控制上，每种补偿方法同样重复进行3次实验，结果如图 8所示。由图 8a可知，相比于四阶前馈，添加$R_{1}$轨迹得到的ILC补偿信号$u_{\mathrm{ILC}}$之后，加速段误差没有改善，反而在进入匀速段之后的误差明显恶化，可见迭代学习补偿不具备变轨迹适应能力。由图 8b—8d可知，无论是M2得到的$\Delta F$，还是本文所提M1得到的ΔF，在变换轨迹后依旧能保持较好的控制精度。此外，在同阶次时与M2相比，M1整定得到的控制器的轨迹跟踪精度更好。

3.3 结果分析

前述实验验证了本文所提算法有效，本节就上述结果展开进一步的分析。

1) 补偿前馈$\Delta F$中超前时移量d的选择原则。

理想的补偿前馈可表示为

$ \Delta F=P^{-1}-F \text {. } $

(29)

当选用四阶前馈时, $F$的超前阶次为4, 故需要获取$P^{-1}$的超前阶次。可通过零相位误差跟踪控制(zero phase error tracking control, ZPETC) 等模型求逆方法^[6]估计其稳定逆模型的超前阶次$d_{\text {ZPETC }}$, 所以补偿前馈中的超前时移量$d=$ $\max \left(d_{\text {ZPETC }}, 4\right)-p$, 其中$p$为参考轨迹的导数阶次。但是对于实际被控对象, 考虑到存在建模误差以及实际系统的延迟, 选择略高的阶次$d \geqslant$ $\max \left(d_{\mathrm{ZPETC}}, 4\right)-p$, 控制效果可能会更好。此外亦可通过Monte Carlo实验的方式确定$d$，详见文[18]。

2) 基于ILC信号拟合的$\Delta F$参数整定方法在低阶次拟合时控制效果较差的原因。

从信号拟合的角度看, 式(7) 中理想的迭代学习信号$u_{\mathrm{LLC}}^{*}$可表示为

$ u_{\mathrm{ILC}}^{*}=\left(S_{\mathrm{p}}\right)^{-1}\left(S r-S_{\mathrm{p}} F r\right) . $

(30)

M2相当于一个输人为$r^{p}$、输出为$u_{\mathrm{ILC}}^{*}$的最小二乘拟合, 即最优化的目标是信号拟合残差。但是该拟合残差不能直接反映最终的轨迹跟踪误差, 从图 3可知, 该拟合残差最终通过$S_{\mathrm{p}}$作用在轨迹跟踪误差上。本文所提M1方法, 如式(11) 所示, 相当于输人为$S_{\mathrm{p}} r^{p}$、输出为$S r-S_{\mathrm{p}} F r$的拟合，拟合残差直接等于轨迹跟踪误差(严格地讲，变为凸优化问题后是广义残余重复误差)，因此本文所提M1要优于M2。

4 结论

针对光刻机工件台前馈控制，本文提出了一种基本前馈加额外有理分式补偿前馈的控制结构，减少了前馈控制中经典的四阶前馈拟合模型不准确对控制性能造成的影响。在此控制结构基础上，本文提出了一种基于数据驱动的额外补偿前馈的参数整定方法，通过优化广义残余重复误差将原非凸优化问题转化为凸优化问题，并给出了迭代过程中梯度、Hessian矩阵的无偏估计。通过光刻机工件台实验，证明这种前馈控制结构有效，并验证该参数整定方法具有可行性。相比于已有的基于ILC信号拟合的补偿前馈，本文研究方法在选用阶次较低补偿控制器时依然能够获得较好的轨迹跟踪控制效果。