王义冬，李庆春，赵 慧
(1.中国人民解放军 92228部队，北京 100074; 2.海军装备部装备项目管理中心，北京 100074；
3.北京振兴计量测试研究所，北京 100074)
微结构光学阵列(micro array structure)能够实现普通元件难以实现的阵列、微小、集成、成像和波前转换等新功能，因此获得越来越广泛运用[1]。惯性约束聚变(ICF,inertial confinement fusion)和人工复眼(ACE,artificial compound eye)成像系统中均存在微结构光学阵列的应用：ICF是以激光作为驱动源，实现受控热核聚变的一项创新技术[2-6]；
ACE是高精度导弹常用的导引头。在ICF中微结构光学阵列实现对高能激光束的调制控制，在ACE中微结构光学阵列实现校正像差和色差，提高导弹精确制导能力和突防能力[7]。
自由曲面(freeform surface)的应用使得光学成像系统朝着超视距、高分辨、大视场、宽光谱和灵巧化的方向发展[8]。光学成像系统在军事侦查、导弹导航、地质遥感、救灾避险等方面得到极大地运用，成为各国角逐的重要领域。麦克斯韦望远镜(JCMT, James Clerk Maxwell telescope)中的远红外相机[9-10]，其采用9块光学复杂曲面金属铝反射镜， 与上一代SCUBA相机比，该相机视场角提升了12倍，扫描速度提升了上千倍，结构却更加紧凑，充分证明了自由曲面在光学成像系统应用中的巨大优势。
在多周期、大振幅微结构光学阵列的加工中以及自由曲面高频非对称误差的补偿加工中，一般情况下要求刀具的带宽大于100 Hz，行程在1 mm以内，故传统的单点金刚石车削(SPDT， single point diamond turning)难以保证制造精度和效率，快速刀具伺服系统(FTS, fast tool servo)提供了一种高精度、高效率的加工方法。
FTS是安装在超精车机床上，其由超精车机床的X轴、C轴、Z轴和微进给轴组成，FTS的系统性能指标主要取决于微进给轴，该轴输出位移取决于超精车床的主轴位置和X轴的位置。
FTS加工出来的面形是由C轴、X轴和微进给轴联动切削而成，理论上可以加工出任意的微结构光学阵列和自由曲面，这是因为FTS推动刀具运动部分的质量较超精车Z轴小，可以实现高频往复运动，FTS的行程通常小于1 mm，带宽一般在1 kHz以内。超精车机床的X轴、C轴和Z轴联动加工面形就是慢刀伺服系统(STS,slow tool servo)[11]，STS相比较于FTS其运动频响较低，其频响在100 Hz以下，存在动态响应能力不足、加工效率低、加工时间长等缺点。
随着微结构光学阵列以及光学自由曲面的在军事国防、国民经济中的应用，对兼顾大行程、高频响和高跟踪精度FTS的需求较为迫切，而FTS的控制能力是制约FTS精度的关键技术之一。压电陶瓷驱动的FTS，其控制难点主要是由压电陶瓷本身特性造成的，压电陶瓷具有输入电压和输出位移的迟滞非线性、蠕变、多值性和记忆性[12-13]等特性，这给压电陶瓷驱动的FTS伺服控制带来了挑战。
压电陶瓷的迟滞非线性是制约压电陶瓷FTS控制精度的最重要因素[14]。迟滞非线性可以通过电荷电流控制来校正[15]，但该方案不仅可能导致漂移和饱和问题，还可能大幅降低工作范围。另一种解决方案是建立迟滞模型，以便建立前馈补偿算法来校正迟滞非线性[16]。现在建立迟滞模型的方法主要有两种，一种是基于压电陶瓷工作原理的物理建模，另一种是不考虑其工作原理将迟滞曲线看成数学问题的现象建模。现象建模得到了广泛应用，其简单性和准确性使得在工程应用上具有物理建模不具备的优势。常见的数学模型有Preisach模型[17]、广义麦克斯韦滑移滞回模型[18]、Duhem模型[19]、Dahl模型[20]、Chua-Stromsmoe模型[21]、Prandtl-Ishlinski模型[22]、Bouc-Wen模型[23]、Jiles-Atherton模型[24]、椭圆拟合模型[25]。当压电陶瓷运动频率增加，其迟滞非线性表现出和频频率相关的特性，Xiao等人[26]、Deng等人[27]、Tian等人[28]在频率相关迟滞模型方面进行了较多研究。为拟合迟滞模型的曲线，需要估计出迟滞模型中的参数，常用的参数辨识方法有最小二乘法、粒子群算法、神经网络算法、遗传算法、免疫算法和蚁群算法，以及基于上述算法的众多改进算法。
PID控制在压电陶瓷驱动的FTS的伺服控制中具有重要工程实践意义，PID控制的简便易实现及对系统模型的准确性没有过度要求，使得PID控制成为FTS伺服控制中的主要控制手段[29]。这是因为压电陶瓷的模型易受温度、使用时间和环境等因素影响。实践中压电陶瓷的模型参数不是固定的，使得一些复杂控制算法的简便性、时效性和经济性较差。但是PID控制不能抑制非线性误差，故前馈控制和PID的复合控制算法在工程上大量应用。此处的前馈控制主要是补偿系统的非线性误差，前馈控制能够大幅增大系统的响应时间，在中高频工况下具有较大优势，而且前馈控制可以补偿掉系统的非线性部分。对迟滞模型求解，然后再求得迟滞逆模型，进而对FTS系统进行基于逆模型的前馈补偿。
本文主要对压电陶瓷迟滞非线性误差进行实验分析，将迟滞非线性误差分为频率无关迟滞现象和频率相关迟滞现象。接着对BW和PI的频率无关迟滞模型进行修正和对比，确定了采用PI模型描述本文的频率无关迟滞现象。然后设计基于Hammerstein模型的频率相关迟滞模型。提出了压电陶瓷驱动器迟滞非线性误差的建模方法，并分析了其有效性和准确性，给FTS伺服控制提供了一种实用的前馈控制器。
为了提高FTS系统的控制精度，则须对FTS的驱动元件(压电陶瓷)的工作特性具有准确的认识。通过对压电陶瓷输入不同信号的开环性能测试，分析了其输出位移的非线性现象，其结果如图1～3所示。
对压电陶瓷输入频率为1 Hz幅值为5 V的正弦信号，将压电陶瓷的时间-位移曲线转化为相应的电压-位移曲线，如图1所示，此曲线为类似椭圆形，该现象表述了压电陶瓷的非线性。图1所示现象称为压电陶瓷的迟滞现象，也反映压电陶瓷的多值性。
图1 1 Hz正弦输入信号的压电陶瓷非线性效应
对压电陶瓷输入相同频率为1 Hz幅值分别为8 V、7 V、6 V、5 V、4 V、3 V、2 V和1 V的正弦信号，绘制出电压-位移曲线如图2所示，从图2中可得压电陶瓷的迟滞现象会随着幅值的增大而增大，而且还可以看出压电陶瓷运动的第一个周期和之后的周期有差异，即初载曲线和之后的曲线有区别，这也是压电陶瓷的蠕变和记忆性造成的，压电陶瓷使用之前预热一段时间可以避免初载曲线的问题，故本文将不再讨论该问题。
图2 不同幅值信号下的压电陶瓷非线性效应
对压电陶瓷输入相同幅值5 V，频率分别为0.5 Hz、1 Hz、5 Hz、10 Hz、100 Hz、200 Hz和400 Hz的正弦信号，绘制压电陶瓷的电压-位移曲线如图3所示，从图3可以看出压电陶瓷的迟滞现象会随着频率的增大而增大，从而表现出压电陶瓷的迟滞现象是和输入信号频率相关的，其中10 Hz以上的迟滞曲线和10 Hz以下的迟滞曲线具有较大差异，本文将10 Hz以下的迟滞现象认为和输入信号频率是无关的(以下简称频率无关)，将10 Hz以上的迟滞现象认为和输入信号的频率是相关的(以下简称频率相关)。
图3 不同频率信号下的压电陶瓷非线性效应
从以上分析可知，压电陶瓷的迟滞非线性表现出和幅值、频率相关的复杂性，这种非线性误差严重影响了压电陶瓷的伺服控制精度。并且这种非线性误差的最大值可以超过输出行程的15%，故不能通过传统、便捷和线性的控制算法(如PID算法)解决，只能通过前馈算法进行补偿。设计前馈算法的前提是可以准确获得迟滞非线性的数学模型，然后对模型求逆就可得到前馈控制器。
2.1 Bouc-Wen模型及其修正
Bouc R于1967年提出了一种描述位移与恢复力之间迟滞关系的数学模型，在此基础上由Y.K.Wen完善[30]，最终形成用非线性微分方程描述系统迟滞曲线的Bouc-Wen模型，其具有控制器简单的优点，该模型的标准表达式为[1]：
(1)
将式(1)转化为描述压电陶瓷迟滞曲线的公式如下：
(2)
式中，y(t)为压电陶瓷位移输出，u(t)为压电陶瓷的电压输入，其中n具有控制曲线平滑性的作用。
在压电陶瓷的前馈控制中，需要将标准BW模型编写入PowerPMAC运动控制器，故需将式(2)进行离散化处理以便于编程。其离散化公式为：
(3)
标准BW模型迟滞曲线是关于中心对称的，而压电陶瓷的迟滞曲线是复杂的、非对称的。本文在标准Bouc-Wen模型的非线性项 里增加若干非对称项，对标准BW模型修正使之具有非对称性，称之为非对称BW模型，其公式如下[31]：
(4)
式中，β1、β2、β3、β4、β5、β6为新引入的参数用来改变标准BW模型的对称性。对式(4)进行离散化处理以便于编程。其离散化公式为：
(5)
2.2 Prandtl-Ishlinskii模型及其修正
PI模型最早是用来描述塑性弹性变形[22]。PI模型具有解析形式、且有逆模型公式以及模型参数易于辨识等优点，但是无法表述频率相关的迟滞模型，其在压电陶瓷迟滞非线性控制方面具有重要应用。
双边play算子，如图4所示[12]，其输出定义为：
图4 双边play算子
(6)
fr(u,x)=max(u-r,min(u+r,x))
(7)
式中，当u(t)∈Cm[0,tE]，Cm[0,tE]为分段连续函数的集合,0=t0 在描述压电陶瓷的迟滞非线性时由于双边算子中的输入始终具有u(t)≥0的关系，故采用单边play算子可以更简洁的描述压电陶瓷的迟滞现象，如图5所示[13]，其输出定义为： (8) fr(u,x)=max(u-r,min(u,x)) (9) 式中的各参数及含义和双边play算子的一致，此处不再赘述。 图5 单边play算子 利用不同阈值r的单边play算子加权叠加而成的PI模型可以表示为： (10) 式中，y(t)为PI模型输出，w(t)为不同阈值的单边Play算子对应的密度函数，Fr[u](t)为单边算子。 PI模型需要编写入PowerPMAC运动控制器中，则需将式(8)～(10)进行离散化处理以便编程，其离散化公式为： (11) 式中，y(k)为PI模型输出值，u压电陶瓷输入值，ri是PI迟滞模型的阈值，wi是PI迟滞模型的密度函数，Npi为单边play算子的个数。 本文提到一种基于Hammerstein模型的频率相关迟滞模型来描述压电陶瓷的频率相关迟滞现象[32]。Hammerstein模型由一个非线性模块和一个线性模块组成，其模型结构如图6所示。 图6 Hammerstein模型示意图 Hammerstein模型中的非线性模块N(·)定义为PI迟滞模型，Hammerstein模型中的非线性模块G(z)定义为ARX模型，其模型公式为： (12) ARX模型中的参数需要进行辨识，具体辨识方法是：将压电陶瓷的电压输入信号u(t)通过PI迟滞模型，得到Hammerstein模型的中间变量yu(t)； 对于非对称Bouc-Wen模型，其中d、α、β1、β2、β3、β4、β5、β6、n八个参数需要用智能算法对其估计。对于修正的Prandtl-Ishlinski模型，其中Npi个密度函数wi需要用智能算法对其估计。本文采用粒子群算法(PSO，particle swarm optimization)作为估计上述参数的智能算法。 粒子群算法是由Eberhart和Kenndy在1995年提出[33]。PSO本质上是一种随机搜索算法，一种智能优化技术。该算法具有概念简明、参数设置少、实现容易、精度高和收敛快等优点，广泛应用于优化问题和搜索问题中。相比其他传统优化算法，PSO收敛到全局最优解的概率较大，它适合在动态、多目标优化的环境中寻优，具有更快的计算速度和全局搜索能力，对待优化的函数没有严苛要求(如可微分、时间连续等)，对多变量、高度非线性、不连续及不可微分的情况优势更加突出。 PSO初始化为一群随机粒子，然后通过迭代找到最优解。在每次迭代中，粒子通过跟踪两种极值完成迭代更新：一种是个体极值，也就是粒子本身所找到的最优解； 为了避免陷入局部最优解，本文采用随机权重的PSO。其原理是将标准PSO中的随机权重w设定为随机数，当粒子在接近最优解时，随机产生的w可能为较小值，由此加快收敛速度，还能避免w线性递减不能收敛到全局最优解的问题。 随机权重的PSO的计算步骤如下： 1)设置模型参数和最大迭代次数M，随机设置每个粒子的位置和飞行惯性(速度)。 2)计算每个粒子的适应度值。 3)获取当前的个体极值和全局极值。 4)更新粒子的飞行惯性(速度)和位置。 5)更新惯性权重。 6)判断是否达到算法停止条件(误差足够小或达到最大迭代次数)，如果达到条件则退出算法并输出结果，否则返回到第2)步继续迭代。 采用图2所示的迟滞曲线，设置随机权重粒子群算法模型参数如表1所示，其中采样时间为0.1 ms。 表1 随机权重PSO算法模型参数 表中，Npso为粒子群的粒子数，c1为个体学习因子，c2为群体学习因子，wmax为惯性权重最大值，wmin为惯性权重最小值，σ为随机权重方差，M为最大迭代次数。 定义随机权重粒子群算法中的适应度函数，其公式为： (13) 定义以下函数来评价拟合效果： (14) (15) NMAX为归一化最大拟合误差，NRMS为辨识精度(归一化均方根误差)。 对标准BW中的d、α、β、γ、n五个参数进行辨识，其中n限制为正整数，随机权重PSO的模型参数如表1。辨识结果如图7所示。 图7 标准Bouc-Wen模型辨识结果 标准BW模型辨识结果为，Frms=0.522，NMAX=2.754%，NRMS=1.848%。从中可以看出辨识效果较差，对FTS所要求的精度而言，该误差是不能接受的。 为了克服标椎BW模型的缺点，采用非对称BW模型来辨识迟滞曲线。利用随机权重粒子群算法对非对称BW模型中的d、α、β1、β2、β3、β4、β5、β6、n九个参数进行辨识，其中n限制为正整数，随机权重粒子群算法模型参数和表1一致，辨识结果如图8所示。 图8 非对称Bouc-Wen模型辨识结果 非对称BW模型辨识结果为，Frms=0.282，NMAX=1.739%，NRMS=0.997%，辨识精度相对于标准BW模型有了较大提升，但是拟合曲线的形状不理想。 为了进一步提高辨识精度，采用修正的PI模型拟合迟滞曲线。利用随机权重权重算法对式(11)中的Npi个密度函数wi进行辨识，为了兼顾模型精度和PI模型在控制系统中的处理速度，Npi的值取为10，随机权重粒子群算法模型参数和表1一致。辨识结果如图9所示。 图9 Prandtl-Ishlinskii模型模型辨识结果 PI模型辨识结果为，Frms=0.111，NMAX=0.650%，NRMS=0.392%，辨识的变量参数如表2中的Case2所示。 更改随机权重粒子群算法中的粒子群数量和最大迭代次数，重新对式(11)辨识，结果如表2所示。对比表中的Case1、 Case2和 Case3，可知最大迭代次数的增加并不一定会使辨识精度增加，其一般取值在[300,600]之间； 为了验证随机权重粒子群算法的可靠性，将辨识算法更换为遗传算法(GA, genetic algorithm)[30]。具体利用Matlab软件中的遗传算法工具箱对式(11)中的10个密度函数wi进行辨识，对这10个参数的区间全部设置在[-1×10-16,1×10-16]之间，设置种群大小为50个。 PI模型由GA算法辨识的结果为，Frms=0.111，NMAX=0.653%，NRMS=0.393%，辨识的变量参数如表3所示，相对比图12的辨识结果和辨识精度并没有提升。对比随机权重粒子群算法和MATLAB遗传算法工具箱对式(11)的辨识情况可以发现，在辨识精度和运算时间方面随机权重粒子群算法具有明显优势，证明了随机权重粒子群算法的可靠性和实用性。 综上所述，随机权重粒子群算法的最大迭代次数一般取值在[300,600]之间、粒子群数量一般取值在[30,80]之间，这两个参量对辨识精度的影响有限，并且将随机权重粒子群算法与Matlab遗传算法工具箱的辨识情况进行比较，以上证明了随机权重粒子群算法可靠性和实用性； 表2 修正PSO辨识Prandtl-Ishlinskii模型参数的结果 表3 GA辨识Prandtl-Ishlinskii模型参数的结果 由图5所示分析结果，可知压电陶瓷的迟滞现象是和频率相关的。对压电陶瓷进行1 Hz至100 Hz，幅值为5 V的扫频实验，绘制如图14所示的迟滞曲线。在10 Hz以下时迟滞模型的建立按照频率无关的迟滞模型处理，当压电陶瓷频率大于10 Hz时，频率的影响将不能忽略。 用表2中Case2所述的PI模型去拟合图10的迟滞曲线，得到的拟合效果如图11所示，其拟合误差的最大值为4.692 μm，误差的均方根值为1.827 μm。从中可以看出PI模型的拟合曲线几乎全部重合为一个迟滞环，这是因为PI模型为频率无关的迟滞模型，故不能用来描述频率相关迟滞模型。 图10 1～100 Hz开环扫频的迟滞曲线 图11 1～100 Hz开环扫频的迟滞曲线 本文采用一种基于Hammerstein模型的频率相关迟滞模型来描述压电陶瓷的频率相关迟滞现象，其中 定义为PI迟滞模型，其公式为(11)，公式中的参数为表2中Case2的数据。利用Matlab的系统辨识工具箱对ARX模型中的参数进行辨识，可得辨识精度为97.73%的ARX模型，其表达式为： (16) PI迟滞模型和式(16)所示的ARX模型共同组成Hammerstein模型。 用Hammerstein模型拟合图11所示的迟滞曲线，其拟合误差的最大值为0.531 μm，误差的均方根值为0.218 μm。相比较于只用PI模型拟合，误差的最大值降低了88.683%，误差的均方根值降低了88.068%，说明Hammerstein模型在描述频率相关的迟滞模型时具有较好效果。Hammerstein模型为描述压电陶瓷频率相关迟滞现象提供了一种简单实用的方法。 本文为了补偿压电陶瓷迟滞非线性误差，实验分析了该非线性误差的特性，将其分为频率无关现象和频率相关现象。利用随机权重粒子群算法分别对标准BW模型、非对称BW模型、修正PI模型进行辨识，对应的辨识精度分别为1.848%、0.977%、0.392%，最终选择PI模型作为压电陶瓷的频率无关迟滞模型。设计了基于Hammerstein模型的频率相关迟滞模型，相比PI模型对频率相关迟滞曲线的拟合，其误差的最大值降低了88.683%，误差的均方根值降低了88.068%。本文提出了压电陶瓷驱动器迟滞非线性误差的建模方法，并分析了其有效性和准确性，给FTS伺服控制提供了一种实用的前馈控制器。
2.3 基于Hammerstein模型的频率相关迟滞模型
将yu(t)和扫频实验的结果y′(t)作为ARX模型的输入和输出，利用Matlab的系统辨识工具箱(system identification)对ARX模型中的参数进行辨识，可得ARX模型。PI迟滞模型和得到的ARX模型共同组成Hammerstein模型，该模型就是压电陶瓷的频率相关迟滞模型。3.1 随机粒子群算法的设计
另一种是全局极值，也就是整个粒子群当前找到的最优解。
3.2 频率无关迟滞曲线的辨识
对比表中的Case2、 Case4和 Case5，可知粒子群数量的增加并不一定使得辨识精度的增加，其一般取值在[30,80]之间。
采用随机权重粒子群算法分别对标准BW模型、非对称BW模型、修正PI模型进行辨识，对应的辨识精度分别为1.848%、0.977%、0.392%，本文最终选择PI模型作为压电陶瓷的迟滞模型。
3.3 频率相关迟滞曲线的辨识
