自学习模型

2024-07-11

自学习模型（共11篇）

自学习模型篇1

一、引言

随着大数据 (Big Data) 技术和数据密集型科学的发展, 数据已经渗透到各个行业和业务功能中, 并已成为生产的一个重要因素。美国国家科学基金会提出了教育大数据研究的目标[1]:更好地了解人们在智慧环境中学习;通过对创新技术的教学设计和教学工具的学习过程, 提高学习和测试方法;将大数据融入学习环境中。在教育大数据理论的机遇到来时候, 带来了更多困难和挑战, 最突出的挑战是如何获得广泛的教育数据来源。随着成熟的信息技术环境的大数据技术, 学生的学习和评价不仅要关注学生成绩的评价, 而且要更加重视学生的形成性评价。文献2中的“促进与信息技术整合的教学建议部分”, 认为对智能化教学环境的建设, 提供优质数字教育资源、鼓励评价软件工具的发展。

二、国内外研究现状

学习经历 (Learning Experience) 指的是在学习过程中发生的任何一种互动或其他经验, 涉及课程, 程序或其他教学参与者。既可以在传统的学习环境中发生, 例如学校、课堂, 也可以在非传统的学习环境, 更能在像学习者的直接请教老师、教授传统的交互式教学情况或非传统的像通过学习游戏和互动应用学习的互动教学情境中。

通过现实生活中的复杂的社会系统的数据的开展探测挖掘, 例如使用包括传感器、GPS定位、智能手机等移动终端收集数据, 识别用户的日常社会背景的活动, 然后将实时数据和历史数据关联, 可以做人际关系的物联网。学习经验得到的数据更接近于“相关性”的实时数据流和历史数据, 这些数据将“数据流”在高地产的形式, 同时综合环境信息协会还需要获取和记录的数据流, 实现环境及个人信息空间。传统意义上的学习经验数据采集通常忽略环境信息的处理, 个人信息的数据主要是从描述性的反馈、绩效信息等任务完成的教育测量方式。

三、教育大数据个性化自适应学习模型

为了实现个性化自适应学习评价数据模型, 我们综合评价内容和评价结果的两个评价过程, 在此基础上, 建立个性化评价模型和自适应评价模型。通过学习活动来确定对学习者的学习行为通过个人评价内容评价的主要信息点, 对评价的个性化的学习过程和结果的基础上的层次评价模型评价和个性化评价的结果, 确定其学习水平。

1. 通过学习和学习活动为支撑的学习行为评价模型, 包括信息的学习行为和学习内容分割聚类分类系统分析不同功能的学习行为的研究内容分析个体学习评价模型的行为, 从而建立一个个性化的学习评价模型。该模型的内容包括四个部分:学习评价、评价、评价、考核和评价, 以及课外资源的评价。该模型的评价过程涉及到学生的正式学习和日常学习活动, 评价内容包括学科知识评价和非学科知识评价, 评价方法涉及定量评价和定性分析。

2. 基于大数据支持的仪器面板、报表和可视化工具的分析大数据个性化自适应分析方法。

学习仪表板提供了数据和报告的可视化分析, 方便个人做出关于教学和学习的决策。学习仪表板包括四个用户视图:学习者的观点, 教育者的观点, 研究者的观点, 以及组织视图。不同的视角是不同的, 但都是相互关联的。根据不同利益相关者的需求, 提供不同的数据显示。

四、结论

综上所述, 教育大数据的环境下, 论文详细分析了网络在线学习环境中师生学习经历数据的获取、识别、分析及交流与生成等基本的计算的理论和方法;阐述了大数据分析技术在网络在线学习模式中学习绩效评价需要解决的关键科学问题;根据研究成果构建了在线学习中学习绩效评价系统。

参考文献

[1]Dean J, Ghemawat S.Map Reduce:Simplified Data Processing on Large Clusters.Proc of 6th OSDI.San Francisco:USENIX Association, 2004:137-150.

[2]Dean J, Ghemawat S.Experiences with Map Reduce:an abstraction for large scale computation Proc 15th Inter-Conf on PACT.Washington DC, 2006:1-2.

[3]Dean J, Ghemawat S.Map Reduce:a flexible data processing tool.Communications of the ACM, 2010, 53:72-77.

[4]D.J.De Witt, R.H.Gerber, G.Graefe, M.L.Heytens, K.B.Kumar, and M.Muralikrishna.GAMMA-A High Performance Dataflow Database Machine.In VLDB’86, pages 228–237, 1986.

自学习模型篇2

空间自回归模型的局部影响分析和运用

由于空间数据的相依特性,使得通常的.删除点诊断异常值的方法不适合采用.为了寻找数据中的异常点和影响点,采用局部影响分析技术,通过引入扰动的方法来发现影响点,最后通过实例说明局部影响分析技术能够有效地发现模型中可能的影响点,并且能够揭示更多的细节信息.

作者：罗玉波吴喜之 LUO Yu-bo WU Xi-zhi 作者单位：中国人民大学统计学院,北京,100872 刊名：统计与信息论坛 CSSCI英文刊名：STATISTICS & INFORMATION FORUM 年，卷(期)： 23(9) 分类号：O212.1 关键词：局部影响空间自回归模型异常点影响点

自学习模型篇3

摘要：针对传统木材纹理分类的准确率低且难度大的问题，依据LBP（局部二值）算子和AD-ABOOST（自适应增强）算法理论，提出了LBP-ADABOOST模型对木材纹理进行识别分类，通过均匀旋转不变特性与原始LBP算子相融合，提取纹理的特征值，结合自适应增强算法，从而训练得到每类纹理所对应的分类器模型参数，构造分类器，实现对木材纹理准确高效分类，实验结果表明相比于BP神经网络，SVM支持向量机等分类算法，该模型的实验结果误差率为4%左右，准确率高，实用性强，

关键词：木材纹理分类；LBP算子；ADABOOST算法；分类器

DOI：IO.15938/j.jhust.2015.02.011

中图分类号：TP391.4

文献标志码：A

文章编号：1007-2683（2015）02-0057-06

0 引言

木材纹理分类是木材优化利用过程的重要部分，木材纹理结构精细复杂、无规律的天然属性，使得纹理分类一直是木材学的前沿课题.针对不同的研究，国内外学者提出了不同的特征提取及纹理分类算法，非负矩阵分解，灰度共生矩阵法，马尔可夫随机场，尺度不变特征变换法等存特征提取方面取得了一定的突破，而常用的分类算法有BP神经网络，SVM支持向量机，决策树，极限学习机等.大多纹理分类模型基于以上算法的结合，且取得了一定的成果，此外，我国对于木材纹理分类的研究起步较晚，初期主要是对国外的经验总结，优化传统的分类算法，现阶段我国的纹理分类技术发展较好，但也存在一定缺陷，主要由于实际应用的训练样本在个体之间存在着差异，导致分类算法对分类结果的差异性较大.如何提高木材纹理分类的准确性至关重要，也将是本文的重点.

近些年，一种简单高效的纹理特征分析方法——局部二值模式（logical binary pattern，LBP）成为了众多学者研究的对象，在描述、提取局部纹理特征方面取得了很好的效果.此外，Adaboost自适应增强算法是通过训练样本特征得到弱分类器，对弱分类器的线性组合得到最终的强分类器，进行分类学习，其在机器学习和数据挖掘方面应用较广.基于以上表述，本文采用LBP与ADABOOST的融合尝试进行木材纹理分类，且这种分类算法在木材纹理分类方面应用极少.经实验，本文提出的基于LBP-ADABOOST模型的木材纹理分类算法达到预期目的，并且正确率明显高于BP神经网络，SVM支持向量机，不失为木材纹理分类提供一种有效方法.

1 特征提取

局部二值模式（logical binary pattern，LBP）最早是由Ojala等在1996年提出，是一种描述图像局部纹理特征的算子，原始的LBP算子定义为在3×3的窗体内，以窗体中心像素为阈值，与相邻的8个像素的灰度值比较，若周围像素值大于中心像素值，则该像素点的位置被标记为1，反之为0.这样，3×3领域内的8个点可产生8bit的无符号数，即得到该窗体的LBP值，并用这个值来反映该区域的纹理信息，如图l所示：

随着LBP算子在图像中的应用，扩展定义一个半径为R（R>O）的圆形邻域，B（B>O）个邻域像素点均匀分布在圆周上面.设定该邻域中心像素值是C，则C可用该邻域内的B+1个像素的函数来定义，即其中：gc为中心像素值；g0，g1，gb-1为B个邻域像素值.中心像素的坐标是（Xc，Yc），则其邻域坐标（Xi，Yi）为：

当坐标（Xi，Yi）不在中心，通过双线性内插法，将邻域的像素值减去中心的像素值计算得到局部的纹理特征像素值C：

假定在实际情况中中心像素值与邻域像素值的差值g0-gc，gb-1-gc独立于gc，则（3）可表示为：

由于c（g。）代表的是中心点的像素值，与图像的局部纹理特征无必然联系，故可以忽略不计，得到：

上式（5）中的像素差值描述了每个纹理模式，若纹理不受像素值单调变化，只须考虑差值符号即可：

由此可以得到一个B为二进制数，乘以相应的权重2i求和得到LBP特征值：

原始LBP算子是灰度不变的，但不是旋转不变的，图像经过旋转后可得到不同的LBP值.为了便于描述图像信息，针对图像旋转之后还可以得到相同的LBP值，Maenpaa等人提出具有旋转不变的LBP算子（rotation invariant LBP），不断旋转圆形邻域得到LBP值，取其最小的作为邻域的LBP值，表述如下式：其中ROI（x，i）是旋转函数，表示将x的二进制数按位循环右移i次，

但随着邻域采点集数的增加，二进制模式的种类急剧增加，达到2B个，这种数量无论是对纹理特征的提取，纹理的分类都是不利的.针对二进制模式的降维问题，Ojala提出一种均匀模式（uniform pat-tern）.均匀模式即根据编码模式出现频率的高低，在圆形二进制编码中，至多有两个0到1或1到0的变化，表示如下：

上式的结果满足U≤2的模式时称为均匀模式，用LBP表示，通过改进，二进制的模式减少，且不会丢失信息，使得原来2B种减少到B（B-1）+2种，降低了特征矩阵的维数.

基于以上两种性质的优越性，可以将LBP的旋转不变性与均匀模式结合得到更好的效果，称为旋转不变均匀模式用符号LBPriu2表示，定义如下：其中：U（LBPb，r）的计算方法是（10）的表述，该模式下不仪保留了图像的纹理特征，而且邻域二进制编码的种类降到了B+2种，大幅减少了特征总量，本文采用了旋转不变均匀模式LBP算子提取特征值，并对比了几种模式下不同分布情况，如图2所示.可以看出LBP算子在均匀旋转不变模式下的特征维数相对较少，同时保留了图像信息，效果最佳，

2 算法原理

自学习模型篇4

话务自学习[1],就是智能学习特定话务的历史变化规律,并根据历史变化规律预测未来话务量的变化趋势,从而对话务数据设置告警的上、下限阈值,达到在第一时间发现话务量异常并告警的目的。

为达到话务自学习的目的,需要对未来话务量进行预测,并根据预测值进行告警。目前各种常用的预测方法有简单函数的拟合预测,如惯性预测[2]、Kalman滤波[3]等、ARMA模型[4,5]、BP模型[6] 、SVM预测模型[7,8]等。ARMA模型作为一种成熟的时间序列分析预测方法,能很好地描述话务量序列,适合在话务量预测中广泛使用。

1AR模型定义和特点

根据话务量数据的特征,可以把话务量作为一种时间序列来进行建模预测。时间序列预测是一种比较成熟的预测理论和方法,尤其是在20世纪70年代以后,时间序列分析理论的日益成熟,为时间序列预测和分析的广泛应用打下了坚实的基础。常用时间序列预测模型有自回归模型(AR)、滑动平均模型(MA)、自回归-滑动平均模型(ARMA)、累积式自回归-滑动平均模型(ARIMA)。

对于p阶的AR模型,其一般表达式为:

$x (n) + \sum_{i = 1}^{p} φ_{i} x (n - i) = u (n) (1)$

式(1)中,x(n)为时间序列,φi是模型的参数,p是模型的阶数,u(n)是均值为0的高斯白噪声,满足N(0,σ2)分布。

AR模型可表示为一白噪声u(n)通过图中所示线性系统得到时间序列x(n)。

2AR模型的建模方法

AR模型的建模流程如图2所示。

2.1数据的预处理

AR模型要求建模的数据平稳,对不平稳的时间序列需要将话务数据进行均值化处理和差分平稳化处理。

时间序列平稳的定义为:如果一个随机过程的均值和方差在时间过程上都是常数,并且在任何两时期的协方差值仅依赖于该两时期间的距离或滞后,而不依赖于计算这个协方差的实际时间,就称它为平稳的。

2.2AR模型定阶

确定AR(p)模型当中p的值称为模型定阶,不同的p值对模型效果的度量有决定性的影响。模型阶数的选取对谱估计特性有很大的影响,阶数过低,会使谱估计分辨力低,得不到理想的结果;阶数过高,不仅会加大计算量,还会在谱估计上产生虚假的细节。常用的AR模型定阶方法有:AIC检验法,最小二乘法,极大似然法和MSE最小值检验法。

现采用MSE最小值检验法,其具体步骤为:

(1) 设定AR模型的最高阶数,设该阶数为n;

(2) 分别求出1阶至n阶的各阶AR模型的参数;

(3) 利用各个识别出的AR模型对某一业务数据进行拟合,再由拟合值和该业务的真实值求出各阶AR模型的MSE值;

(4) 迭代寻找MSE值最小的所对应的AR模型阶数p,即为最优AR模型阶数。

2.3AR模型参数估计

AR模型参数估计的方法总体有直接法和间接法两类:直接法包括最小二乘法、解Yule-Walker方程法、UIrych-Clayton法等;间接法包括LUD法、BSMF法、Burg法等。上述方法中,采用解Yule-Walker具有简单、参数估计无偏、精度高等优点,现在AR模型参数估计方法采用解Yule-Walker方程的方法。

2.4AR模型的稳定性检验

AR模型的稳定性检验可由下式表达

$| φ_{1} | + | φ_{2} | + \dots + | φ_{p} | < 1 (2)$

即模型参数的绝对值之和小于1时,认为模型是稳定的。

3AR模型应用实例分析

在实际工作中,可通过选取一段时间的话务量的历史数据,选择合适的建模算法进行建模,根据模型的预测结果及事先设定的告警阈值,在话务量出现异常时告警,以方便运维人员及早地对异常话务进行分析,查找原因,保证后续话务的正常。

在实验中选取某电信公司的7～8月的某业务的话务数据作为仿真数据。以7月份的话务数据为样本数据进行建模,使用建模结果对8月份的话务数据进行预测,通过MAPE值来判断仿真结果的准确度。

3.1话务数据的预处理

首先选取了7月份某业务的每日话务数据,共31条数据形成时间序列,记为x(n)。该时间序列如图3所示。

通过对话务数据序列观察可知,该业务月份数据显示出呈现一定的周期性规律,在7月份以内(31天)近似以7天(一周)为一周期重复流量。为满足建模序列平稳化需求,初步对原始序列进行差分处理,即Xt-ᐁ7(Xt)。经差分处理后损失7天样本,样本容量缩小为24,序列图初步呈现一种不相关的随机变量序列,可判断其期望和方差为一常数。差分后的序列如图4所示。

3.2模型的确定

通过自相关和偏自相关函数图形,来判断是否适合建立AR模型,如图5所示。

通过自相关和偏自相关图形,可以初步判经差分处理后的序列应该是平稳的,自相关系数ρk在负指数的控制下有逐步收敛到0的趋势,属拖尾情形,而偏自相关系数φkk在k=16以后显著趋于0,属截尾情形。根据平稳AR模型的两个重要特性:自相关系数拖尾与偏自相关系数截尾,由此判断序列适合建立AR模型。

通过MSE最小值检验法来确定AR模型的阶数,易见MSE指标在模型7阶时MSE值趋于平稳。遵循建模精简性原则,考虑建立AR(7)模型。如图6所示。

在AR模型的阶数确定后,即可通过解Yule-Walker方程求出模型的参数φi,从而确定了和该业务相对应的模型。

3.3预测结果

使用模型对8月份的话单数据进行预测,如图7所示。图中蓝色为真实数据差分后的结果曲线,红色为预测结果的曲线,从图中可以看出两条曲线的轨迹有很高的吻合度,说明由AR(7)模型预测出的结果是相当准确的。计算得到MAPE=2.6137261×10-2。

4 总结

通过对AR模型的研究,结合话务自学习对话务模型的要求,提出了一种基于AR模型的话务自学习模型。通过对话务数据进行平稳化差分处理,模型的定阶、模型参数的确定等步骤,得到话务数据的模型,并根据话务模型对未来的话务数据进行预测。仿真结果表明,建立的话务自学习的AR模型和实际的话务数据有很高的一致性,且仿真速度快,适合在实际的工作中应用。

摘要：通过对AR模型的研究,提出了一种基于AR模型的话务自学习模型。通过对话务数据进行平稳化差分处理,得到话务数据的模型,并对未来的话务数据进行预测。仿真结果表明,建立的话务自学习的AR模型和实际的话务数据有很高的一致性,且仿真速度快,适合在实际的工作中应用。

关键词：话务自学习,AR模型,话务预测

参考文献

[1]中国电信计费模型(CTG—BM3.0—06):《网管与监控》分册V1.0.ppt,2009年10月:59—62

[2]陶伟宜.基于回归方法的移动通信话务量预测.邮电设计技术,2004;(5):29—32

[3] Harvey A C.Forecasting,structural time series models and the kal-man filter.Cambridge,UK:Cambridge University Press,1991

[4] Box G E P,Jenkins G M,Reinsel G C.Time series analysis Forecas-ting and control.3rd ed New Jersey:Prentice-Hall,1994

[5]于艳华,王军,宋俊德.应用乘积季节ARIMA模型的话务量预测及结果分析.计算机工程与应用,2009;45(20):99—102

[6]牟颖,王俊峰,谢传柳,等.大型呼叫中心话务量预测.计算机工程与设计,2010;31(21):4686—4689

[7]陈蓉,宋俊德.基于SVM分块回归分析的话务量预测模型.计算机应用,2008;(9):2230—2232

自学习模型篇5

摘要：由于传统BP算法存在收敛速度慢，容易陷入局部极小值等弊端，目前的BP优化算法又使得控制过程变得复杂，继而基于BP神经网络的模型参考自适应控制过程也存在实时性差，收敛性慢，精度不高等不足。现针对改进的BP算法和非线性系统的可逆性，分析设计了一种基于激励函数自寻优的BP网络模型参考自适应控制，并通过Matlab仿真结果表明，在满足控制精度的情况下控制系统中的辨识器和控制器效果都很理想。因此，对工程应用有很大的实际参考利用价值。

关键词：BP算法；神经网络；模型参考自适应控制；激励函数；Matlab仿真

中图分类号：TP183

文献标识码：A

DOI：10.3969/j.issn.1003-6970.2015.07.024

0 引言

在现代实际工业生产中，被控对象存在各种不确定性和时变性，因而使得工业控制过程变得繁琐复杂，针对线性时变系统或非线性系统的控制，人们不断的研究其解决方法，Nare ndra等人提出了神经网络控制和模型参考自适应控制相结合的神经网络模型参考自适应控制（Neural Network Model Reference AdaptiveControlˉˉˉNNMRAC）方法。近来神经网络的研究已成为智能控制研究的热点，因其自身具有自学习的特点，可以有效地解决不确定和复杂的非线性控制系统控制问题。因此将神经网络与模型参考自适应控制相结合，组成基于神经网络的模型参考自适应控制系统，进而使其在复杂非线性过程控制中具有不可替代的优势。目前神经网络模型参考自适应控制系统中应用最广泛的神经网络是BP神经网络。

BP神经网络（Back Propagation Network）是一种多层前向型神经网络也被称为反向传播网络，在BP网络中信号是前向传播的，而误差是反向传播。一般三层BP网络结构就可以使其对有限个不连续点的函数进行逼近，也可以逼近任意非线性映射关系。然而，传统BP神经网络算法存在很多缺点，各种优化改进的BP算法也层出不穷。就目前的BP优化算法，常常忽略算法本身存在的自适应、自学习等特点，改进算法如蚁群算法和粒子群算法大都是直接对BP神经网络中的参数进行训练。文中采用基于激励函数自寻优的方法改进BP神经网络模型参考自适应系统的控制方法，改进后的BP神经网络模型参考自适应控制方法收敛速度快、精度高，系统控制过程中被控对象的输出能够很好地跟随参考模型的输出，具有很好的控制效果，在实际工程中也可以得到很好的验证和应用。

1 BP神经网络结构及其算法改进

BP神经网络已经被证明具有很强的学习能力，能够逼近任意连续有界的非线性函数。一般的BP神经网络包括输入层、隐含层、输出层，其中隐含层可以有多个。其中隐含层和输出层的激励函数通常直接采用Sigmoid函数，其函数表达式为：

上式中β称为Sigmoid函数的斜率参数，不同的β取值，引起曲线的弯曲程度不同，β越大，f（x）图形越陡峭。当斜率参数β接近无穷大时，Sigmoid函数将转化成简单的阶跃函数。但与阶跃函数不同，Sigmoid函数对应于0～1之间的一个连续取值区域，但阶跃函数只对应0和1两个取值。

一般的三层前馈神经网络拓扑结构如图1所示。

输入向量为X_i=（Xl，X2…，X_n）^T，i=1，2，…，n，n表示输入神经元的个数，隐含层的输出向量为y_j= （Y1，y2，…，y_m），j=l，2，…，m，输出层的输出向量为Ok= （Ol，O2….，o_l），k=l，2，…，l。每层之间的权值用w表示，W_ij为隐含层和输入层之间的权值，W_jk为输出层和隐含层之间的权值。

这里对上述BP算法的改进，也就是通过改进激励函数f（x），进而优化神经网络，最终使得基于神经网络的模型参考自适应控制在不增加复杂性及确保精度的情况下，系统性能进一步得到提高改善。由于BP神经网络产生局部极小值的一个重要原因就是误差函数是一个以Sigmoid函数为自变量的非线性函数，而Sigmoid函数存在饱和区，所以改进和优化激励函数对于BP算法的应用是至关重要的。通过实验发现，在函数表达式中增加一个控制参数η，可以控制激励函数的压缩程度。改进的激励函数形式如下：

上式描述的f（x）的定义域为（-∞，+∞），值域为（0，1），函数也是单调的，满足激励函数的条件。

以往出现的改进BP算法学习过程中，η和β的赋值都是经验值，本文使得η是一个可以自适应的参数，就是通过判断网络不断学习过后的权值能否减小网络误差来自动的调整η的值，其调整方法为：

其中，al，E为网络误差，f指迭代次数。在误差信号反向传播时，自适应参数η是随着误差信号不断进行修正的。

此算法可以提高BP网络的收敛速度，同时也避免了陷入局部极小值。文中使其结合模型模型参考自适应控制明显提高了系统的整体控制效果，进一步验证了算法改进的实用性，与传统BP算法相比，改进后的算法在实际运用中更具有意义。

2 神经网络模型参考自适应控制系统结构

典型的神经网络模型参考自适应控制系统结构如图2所示。

图2中NNC（ Neural Network Controller）为神经网络控制器，NNI（ Neural Network Identifier）为神经网络辨识器，r为参考输入，u为NNC的输出，Y_m和y分别为参考模型和被控对象的输出，e_c是参考模型输出和被控对象输出之差，e_i是被控对象输出和辨识器输出之差，NNC的权值修正目标是使e_c达到系统设定值（理想值为零），NNI的目标也是使e_i尽可能最小（理想值为零），且为NNC传递梯度信息。

神经网络辨识器NNI的训练误差表示为，其中，y（k）当前k时刻被控对象的输出数据，为下一时刻的预测输出数据。则辨识器的调整规则就是使误差E_i尽可能小，E_i表示为：

神经网络模型参考自适应控制系统的控制目标在于使被控对象的输出y与参考模型的输出Y_m渐近的匹配，即

其中，s为一个给定的小正数。

神经网络控制器NNC的训练则由误差e_c=y_m-y来训练，训练准则如上式（8），控制系统中神经网络辨识器和控制器的学习算法就采用改进后的BP算法。

在神经网络模型参考自适应控制系统的控制策略设计中，改进的BP算法能够在满足系统控制规律符合要求的情况下，使得神经网络模型参考自适应控制效果更好。虽然改进的BP算法是激励函数自寻优的自适应方法，不能够使神经网络辨识器NNI进行离线训练，但是快速的BP算法仍然可以使网络具有很好的实时性。首先在线训练辨识器，待参数训练好以后，再进行控制器NNC的训练，最终可以保证被控对象的输出y很好的跟踪参考模型的输出Y。

3 仿真实例研究

3.1 改进的BP算法验证

本文采用BP神经网络进行预测控制来验证改进算法的有效性。利用简单的一组样本训练集和样本目标集进行神经网络的训练，再给定一组输入样本数据，观测输出层输出数据和误差。分析样本数据设计BP神经网络结构为3个输入、2个输出、隐含层的神经元数目为8。网络学习次数为100次，目标误差设置为0.001。

使用MATLAB软件进行网络训练，传统BP算法的网络训练过程收敛情况如图3所示，经过56步循环达到了网络误差要求的精度。改进的BP算法网络收敛情况如图4所示，仅需要10步就达到了误差精度要求，其中a和b取值分别为0.9和1.5。输出误差和网络实际输出数据在表1中展示，直观的看出，改进后的BP网络可以得到更有效的输出。

根据得到的误差收敛曲线比较看出，改进后的BP算法所用训练步数少即需要的训练时间少，说明收敛速度明显加快。

测试输出结果如表1所示。

从表中可以直观清晰看出改进后的BP算法实际输出误差明显减小，提高了算法精度。

3.2 改进的神经网络模型参考自适应控制仿真实例

结合参考文献中提到的污水处理的例子进行改进算法的验证。污水处理系统结构图如图5所示，

在污水处理系统控制结果是否达标，主要是通过需氧量（OD）、溶解氧（DO）等几个重要参数来衡量。本例中为了提高污水处理效果，系统控制目标设置为使误差e_c控制在±0.05mg/L以内，污水处理控制系统中采

用离散的参考模型：

y_m（k）= 0.375 y_m（k-1）+0.623r（k）

其中，控制输入r（k）=2为系统给定的阶跃信号。

污水处理系统的实验仿真中，BP网络辨识器设定4个输入变量和1个输出变量，隐含层包含10个隐节点，对于BP网络控制器取3个输入层节点，隐含层的节点数为6。根据BP神经网络控制器和辨识器的改进算法，采用MATLAB进行仿真，取采样周期t_s=0.OOls，这里a取0.8，6取1.5，仿真结果如图6所示。

在图6 （a）中第一条线为控制输入r，中间的第二条曲线代表参考模型的输出Y_m，最下边的曲线代表污水被控对象的输出y，图6 （b）中的曲线代表误差e_c（系统实际输出与参考模型输出之差）的变化。从图中可以分析看出，改进后模型参考自适应控制方法在该控制系统中的控制效果很好，氧的溶解浓度（DO）保持在2mg/L左右，参考输出和实际输出最终相吻合，误差e_c控制在±0.05mg/L以内，因此仿真结果满足控制系统的控制要求。

4 结语

自学习模型篇6

近年来,随着冷轧产品带钢厚度需求越来越薄,板形质量要求越来越高,轧机基本设定模型的精度需求也随之增高。轧制力作为轧制过程控制基本设定模型的核心参数,其预报精度不仅影响板形质量、板厚精度,也对整个冷轧过程的工艺控制精度起着至关重要的作用[1]。而轧制过程是受多种因素影响的非线性、耦合、实时性的复杂过程,因此按照传统数理模型建立的轧制力预报模型具有很大的局限性[2]。文献[3]把小脑模型神经网络CMAC用于未知的动态控制系统中,通过严格的公式推导证明了CMAC用于闭环控制系统具有可靠的稳定性和良好的逼近能力,但是传统CMAC神经网络的学习过程中误差平均分配腐蚀了历史学习结果,限制了它的学习速度和收敛精度[4]。为了进一步提高轧制力预报模型的精度,笔者利用权值更新次数的倒数与单个样本本次激活的地址更新次数倒数和的比,作为网络权值更新的信度,建立基于信度分配的小脑模型CA-CAMC网络与轧制力自学习相结合的轧制力预报模型。仿真实验证明该模型具有更高的预报精度和较快的收敛速度,能够更好地满足现代工业越来越高的控制精度需求。

1传统轧制力及自学习模型

1 . 1轧制力模型

带钢冷轧时金属不仅会发生塑性变形,而且在入口及出口处还存在弹性变形,因此目前计算冷轧过程控制轧制力时广泛采用考虑了两种变形的工程计算公式———布兰德福特希尔公式:

其中,

上述式中: Fpc为轧制力; b为带钢宽度; Fk为动态变形抗力( 由静态变形抗力及延伸率计算) ; k为张力影响因子; Dp为摩擦力影响因子;为带钢与轧辊的接触长度,R'为工作辊弹性压扁半径,H和h分别为轧机入口、出口带钢厚度; Ftb为后张力; Ftf为前张力; r为总的压下率,r = ( H - h) /H; μ 为摩擦因数。

1 . 2轧制力自学习模型

轧制力模型的自学习功能就是利用大量的在线数据及反馈数据不断在线更新模型的修正系数,以适应系统状态的变化,逐渐减小传统数据模型计算的理论轧制力值与现场实际轧制力的误差,从而提高以后带钢轧制力的设定精度[5]。由于现场生产时带钢规格变化多且存在各种参数波动,如果选取一个常数作为修正系数,则无法实现设定轧制力时刻快速逼近实际轧制力,因此需要根据历史轧制数据和实际测量数据计算轧制力自学习系数,并保存到特定数据表中,以适应现场的多变性。轧制力自学习系数计算如下:

上述式中:zp为轧制力自学习系数,初始值为1;z-1p为前一钢卷的轧制力自学习系数;β为自学习系数的平滑常数;z-1pa为前一钢卷的实际轧制力自学习系数;Fp为自学习预报轧制力;F-1pa为前一钢卷的实际轧制力;F-1p为前一钢卷的预报轧制力。

2CA-CMAC网络轧制力预报模型

2 . 1CA-CMAC预报模型

影响轧制力计算精度的主要因素有带钢宽度、厚度、张力、变形抗力及工作辊半径等,其中变形抗力主要由带钢屈服强度决定[6]。因此本文把带钢宽度b、入口带钢厚度H、出口带钢厚度h、前张力Ftf、后张力Ftb、带钢屈服强度yp、工作辊半径R作为CA-CMAC轧制力自学习预报模型输入X的部分参数,另外为了保持传统轧制力自学习模型的完整性,把其预报值Fp作为CACMAC轧制力自学习预报模型输入X的另一个参数,轧制力预报值Fpca作为模型输出,建立基于信度分配的小脑模型CA-CMAC网络轧制力预报模型,模型结构如图1 所示。

图1 由线性映射X→A和非线性映射A→W两个映射组成,A为状态激活联想空间,W为物理地址空间。在映射X→A中,输入空间X中的每个输入向量首先被离散化到状态空间S,然后每个状态激活联想空间A中的C个存储单元经编码形成C个地址。对于每个输入,A中只有C个单元为1,而其余的均为0,因此A是一个稀疏矩阵。另外随着输入的增加A所需要的存储空间也会急剧增大,在实际应用中很难实现,因此我们可以使用哈希编码的方式对A的空间进行压缩,将A映射到一个小得多的物理地址空间W中。在这个映射过程中,输出数据随机分布存储在这C个地址单元格中。假设压缩后物理存储空间共有Nm个存储单元,则第i个输入对应的输出Fpcai由C个被激活的权值求和得到,其表达式如式( 7) 所示:

式中: aj( x) 为输入向量x激活的联想空间的第j个地址的状态标志,若第j个地址被激活,则aj( x) 其值为1,否则为0; wj为第j个存储单元的权值。

S1~S8—8个输入向量对应的状态空间;a1~aN—联想空间地址是否被激活的状态标志,其值为0或1;Δw(t)—为第t次学习时权值更新变化量;Fpa—实际轧制力。

2 . 2CA-CMAC模型的学习算法

CA-CMAC模型的学习过程就是根据期望输出和实际输出之间的误差不断调整权值的过程。传统CMAC网络学习过程中把误差平均分配到C个存储单元作为权值更新变化量,并没有一种具体的方法决定哪一个存储单元应对本次误差负有较大的责任,能够记录且应用的信息只有在学习过程中每个单元调整的次数。文献[7]中提出存储单元更新的次数越多,存储的信息越精确,权值调整量越少,因此笔者把存储单元权值更新次数的倒数与本次激活地址更新次数倒数和的比作为误差分配的信度,对CA-CMAC模型的权值进行更新,如式( 8) 所示:

式中: η 为学习常数; ei为第i个输入的期望输出与模型预报输出的误差; g( j) 为第j个单元地址在t - 1 次循环时更新的次数。

CA-CMAC轧制力自学习预报模型的算法描述如下。

( 1) 初始化网络参数,包括量化空间A和每个输入激活的存储单元C;

( 2) 采用mapminmax( ) 函数归一化输入X和实际轧制力Fpa;

( 3) 根据地址映射函数dj= mod ( aj( x) + j,C) + 1 计算X激活的单元地址; 利用hash ( ) 函数压缩联想空间的地址到物理单元地址;

( 4) 被激活物理地址的更新次数加1,计算本次激活的单元地址更新次数倒数的和;

( 5) 计算第i个期望输出也就是第i个实际轧制力Fpai与第i个预报输出Fpcai的误差ei;

( 6) 根据式( 8) 对网络模型的权值进行更新;

( 7) 判断误差或者迭代次数是否满足结束条件,若满足则终止循环,否则转到第4 步。

3模型测试及结果

为了验证本文建立的基于CA-CMAC轧制力自学习预报模型的精度,利用某厂1420 酸轧实际生产数据对模型进行训练与测试。随机选取141 个钢卷的生产数据,其中120 卷数据为训练集,21 卷为测试数据集。该厂使用5 机架6 辊轧机,分别对1 ~ 5 机架的轧制力预报模型进行训练和测试,误差测试结果如表1 所示。

从表1 中可看出,与传统模型的预报误差相比,基于CA-CMAC网络的轧制力自学习预报模型预报误差较小,且较稳定,即使是轧制力偏小的第5 机架的误差也有一定程度的减小,且模型训练的时间只有10 s,因此能够满足现场生产快速性、连续性的需求。

10 k N

图2 是10 卷钢卷的1 机架实际轧制力与传统模型和CA-CMAC模型预报轧制力对比图。可明显看出CA-CMAC模型的预报值更逼近实际轧制力,经过CA-CMAC网络学习补偿后轧制力设定模型的设定精度得到了提高,更能够满足现代工业高精度的需求。

4结论

本文将传统轧制力自学习模型的输出和轧制力主要影响因素作为CA-CMAC网络的输入,轧制力作为输出,建立了一种基于信用度分配的CA-CMAC网络轧制力自学习预报模型。实验证明,本文提出的预报模型可以大幅度提高预报精度和预报稳定性,更能满足现代工业越来越高的控制精度需求;模型训练时收敛速度快,也能更好地满足现场连续生产及时性需求。另外,该模型具有可移植性,对建立轧制控制系统中其他参数的预报控制模型具有一定的借鉴作用。

参考文献

[1]孙文权,何安瑞,邵建,等.高精度冷轧自动控制系统研究及实践[J].冶金自动化,2015,39(3):44.

[2]高雷,郭立伟,陈丹.冷连轧轧制力模型中变形抗力和摩擦系数的分析[J].轧钢,2013,30(4):12.

[3]COMMURI S,LEWIS F L.Control of unknown nonlinear dynamical systems using CMAC neural networks:structure,stability,and passivity[C]//Intelligent Control Proceedings of the IEEE International Symposium on.Monterey:IEEE Control Systems Society,1995:123.

[4]苏刚,陈增强,袁著祉.小脑模型关节控制器(CMAC)理论及应用[J].仪器仪表学报,2003,24(增刊1):269.

[5]刘挺.攀钢1 220 mm冷连轧机轧制力模型参数自学习分析[J].轧钢,2014,31(5):45.

[6]王飞,刘顺心.冷轧变形抗力模型研究[J].冶金自动化,2014,38(1):60.

自学习模型篇7

目前针对装载机的轨迹跟踪控制问题已经取得一定程度的进展。柳波,等基于模糊自整定PID控制实现了装载机自动铲掘时的轨迹控制问题[1]。龚捷,等在考虑作业阻力补偿的情况下,基于计算力矩法设计装载机铲掘作业的轨迹跟踪控制器,实现铲斗轨迹跟踪控制[2]。

自适应迭代学习控制整合自适应控制和迭代学习控制的优点,基于PD反馈控制器增益自适应调节,迭代学习控制器修正控制力矩,该算法具有快速的渐进收敛性能[3]。拉格朗日法与牛顿-欧拉法是机构动力学建模的两种常见方法,封闭形式的拉格朗日方程是控制器设计和综合的有效工具,而递推形式的牛顿-欧拉方程则是实时控制、求解逆动力学问题的有效方法[4]。采用牛顿-欧拉法推导了装载机工作装置的二自由度动力学模型,在matlab环境下设计自适应迭代学习控制器并进行稳定性证明,对装载机工作装置的轨迹跟踪问题进行仿真控制,以期达到高精度追踪。

1 基于牛顿-欧拉法的动力学方程

将装载机工作装置机器人化处理,视为平面二连杆操作臂机构,建立如图1所示的D-H坐标系。设定动臂与车架铰接点为O0,铲斗与动臂铰接点为O1,铲斗斗尖位置为O2,分别以O0、O1、O2为原点建立坐标系O0x0y0z0、O1x1y1z1、O2x2y2z2。

图1中,θ1、θ2分别为动臂和铲斗相对于车架和动臂的转角;l1=lO0O1、l2=lO1O2分别为动臂和铲斗长度;m1、m2分别为铲斗和动臂的质量;假定动臂和铲斗的质心处于连杆O0O1、O1O2的中心位置。

1.1 牛顿-欧拉动力学方程

根据旋转关节的运动推导关节力矩的过程可细分为外推法和内推法两个步骤。外推法采用牛顿欧拉方程,按连杆编号从小到大的顺序逐次计算连杆的速度和加速度;内推法则按连杆编号从大到小的顺序计算连杆间的力、力矩[5]。运动学方程如下。

外推法求解过程:i=0,1,…,n。

内推法求解过程:i=n,n-1,…,1。

式中:iωi、 $^{i} ω_{\dot{i}}$ 、ivi、 $^{i} b_{\dot{i}}$ 为连杆i的角速度、角加速度、线速度、线加速度相对于坐标系{i}的表达式;i+1iR和 $_{i + 1}^{i}$ R为坐标系{i}和{i+1}之间的姿态变换矩阵,且i+1iR $_{i + 1}^{i}$ R=I;i $\hat{Ζ}$ i为Zi轴方向上的单位矢量;iPi+1为节点i到节点i+1的矢量;iPCi为节点i到连杆i质心位置处的矢量;CiIi为连杆i质心处的惯性张量;mi为连杆i的质量;为连杆i质心上的惯性力和力矩相对于坐标系{i}的表达式;为连杆i-1作用在连杆i上的力和力矩相对于坐标系{i}的表达式;τi为关节i上的旋转力矩。

1.2 装载机工作装置的动力学计算

以图1中装载机工作装置的二连杆机器人操作臂为研究对象,推导其牛顿-欧拉动力学方程。假设装载机铲斗斗尖处不受力的作用f3=0,n3=0;车架固定不旋转有 $ω_{0} = 0, \overset{\cdot}{ω_{0}} = 0$ ;计及重力因素[6]的影响,在式(1)中令 $v_{\overset{\cdot}{i + 1}}_{0} = G - G y_{0} ‚ g = 9.8 m / s^{- 2}$ 。相邻坐标系之间的姿态变换矩阵为(其中ci+1=cos(θi+1)),(si+1=sin(θi+1))。

$\begin{array}{l} _{i + 1}^{i} R = [\begin{matrix} c_{i + 1} & - s_{i + 1} & 0 \\ s_{i + 1} & c_{i + 1} & 0 \\ 0 & 0 & 1 \end{matrix}], \\ _{i}^{i + 1} R = [\begin{matrix} c_{i + 1} & s_{i + 1} & 0 \\ - s_{i + 1} & c_{i + 1} & 0 \\ 0 & 0 & 1 \end{matrix}], i = 0, 1 (3) \end{array}$

连杆在坐标系{i}中的惯性张量为:

根据平行移轴定理[7]得到以刚体质心{C}为原点的坐标系与{i}坐标系间的转换式

${\begin{cases} ^{i} Ι_{x x} =^{C_{i}} Ι_{x x} + m (y_{c}^{2} + z_{c}^{2}) \\ ^{i} Ι_{y y} =^{C_{i}} Ι_{y y} + m (x_{c}^{2} + z_{c}^{2}),^{i} Ι_{z z} = \\ ^{C_{i}} Ι_{z z} + m (x_{c}^{2} + y_{c}^{2}) \\ ^{i} Ι_{x y} =^{C_{i}} Ι_{x y} - m x_{c} y_{c},^{i} Ι_{x z} = \\ ^{C_{i}} Ι_{x z} - m x_{c} z_{c},^{i} Ι_{y z} =^{C_{i}} Ι_{y z} - m y_{c} z_{c} \end{cases} (6)$

式(6)中pc=[xcyczc]T表示连杆质心在坐标系{i}中的位置矢量。

将式(3)—式(6)代入式(1)和式(2)计算装载机工作装置的动力学方程如下(式子太长未列出)。

$\begin{array}{l} ^{C_{1}} Ι_{1} = [\begin{matrix} 0 & 0 & 0 \\ 0 & \frac{m_{1} l_{1}^{2}}{12} & 0 \\ 0 & 0 & \frac{m_{1} l_{1}^{2}}{12} \end{matrix}], \\ ^{C_{2}} Ι_{2} = [\begin{matrix} 0 & 0 & 0 \\ 0 & \frac{m_{2} l_{2}^{2}}{12} & 0 \\ 0 & 0 & \frac{m_{2} l_{2}^{2}}{12} \end{matrix}] (7) \end{array}$

设 $\vec{i}$ 、 $\vec{j}$ 、 $\vec{k}$ 分别为坐标系O0-x0y0z0中的单位向量

${\begin{cases} ^{1} ω_{1} = \dot{θ}_{1} \vec{k} \\ ^{1} \dot{ω}_{1} = \ddot{θ}_{1} \vec{k} \\ ^{1} \dot{v}_{1} = g s_{1} \vec{i} + g c_{1} \vec{j} \\ ^{1} \dot{v}_{C_{1}} = (- \frac{l_{1}}{2} \ddot{θ}_{1} + g s_{1}) \vec{i} + (\frac{l_{1}}{2} \ddot{θ}_{1} + g c_{1}) \vec{j} \\ ^{2} ω_{2} = (\dot{θ}_{1} + \dot{θ}_{2}) \vec{k} \\ ^{1} \dot{ω}_{2} = (\ddot{θ}_{1} + \ddot{θ}_{2}) \vec{k} \\ ^{2} \dot{v}_{2} = (l_{1} \ddot{θ}_{1} s_{2} - l_{1} \dot{θ}_{1}^{2} c_{2} + g s_{12}) \vec{i} + (l_{1} \ddot{θ}_{1} c_{2} + \\ l_{1} \dot{θ}_{1}^{2} s_{2} + g c_{12}) \vec{j} \\ ^{2} \dot{v} C_{2} = (l_{1} \ddot{θ}_{1} s_{2} - l_{1} \dot{θ}_{1}^{2} c_{2} + g s_{12} - \frac{l_{2}}{2} (\dot{θ}_{1} + \dot{θ}_{2})^{2}) \vec{i} + \\ (l_{1} \ddot{θ}_{1} c_{2} + l_{1} \dot{θ}_{1}^{2} s_{2} + g c_{12} + \frac{l_{2}}{2} (\ddot{θ}_{1} + \ddot{θ}_{2})) \vec{j} \\ ^{1} Ν_{1} = \frac{m_{1} l_{1}^{2}}{12} \ddot{θ}_{1} \vec{k} \\ ^{1} F_{1} = (- \frac{m_{1} l_{1}}{2} \ddot{θ}_{1} + m_{1} g s_{1}) \vec{i} + (\frac{m_{1} l_{1}}{2} \ddot{θ}_{1} + m_{1} g c_{1}) \vec{j} \\ ^{2} Ν_{2} = \frac{m_{2} l_{2}^{2}}{12} (\ddot{θ}_{1} + \ddot{θ}_{2}) \vec{k} \\ ^{2} F_{2} = m_{2} (l_{1} \ddot{θ}_{1} s_{2} - l_{1} \dot{θ}_{1}^{2} c_{2} + g s_{12} - \frac{l_{2}}{2} (\dot{θ}_{1} + \dot{θ}_{2})^{2}) \vec{i} + \\ m_{2} (l_{1} \ddot{θ}_{1} c_{2} + l_{1} \dot{θ}_{1}^{2} s_{2} + g c_{12} + \frac{l_{2}}{2} (\ddot{θ}_{1} + \ddot{θ}_{2})) \vec{j} \end{cases} (8)$

${\begin{cases} τ_{1} = (\frac{m_{1} l_{1}^{2}}{3} + \frac{m_{2} l_{2}^{2}}{3} + m_{2} l_{1}^{2} + m_{2} l_{1} l_{2} c_{2}) \ddot{θ_{1}} + \\ (\frac{m_{2} l_{2}^{2}}{3} + \frac{m_{2} l_{1} l_{2} c_{2}}{2}) θ_{2}^{\ddot{2}} - \frac{m_{2} l_{1} l_{2} s_{2}}{2} \overset{\cdot}{θ_{2}}^{2} - \\ m_{2} l_{1} l_{2} s_{2} \overset{\cdot}{θ_{1}} \overset{\cdot}{θ_{2}} + \frac{m_{2} g l_{2}}{2} c_{12} + (\frac{m_{1}}{2} + m_{2}) g l_{1} c_{1} \\ τ_{2} = (\frac{m_{2} l_{2}^{2}}{3} + \frac{m_{2} l_{1} l_{2} c_{2}}{2}) \ddot{θ_{1}} + \frac{m_{2} l_{2}^{2}}{3} \ddot{θ_{2}} + \\ \frac{m_{2} l_{2} s_{2}}{2} \overset{\cdot}{θ_{1}}^{2} + \frac{m_{2} l_{2} g c_{12}}{2} \end{cases} (9)$

2 工作装置轨迹的自适应迭代学习控制

2.1 问题描述

考虑实际中存在不确定性、摩擦及外部干扰的影响,式(9)描述的装载机工作装置动力学方程可写为

$Μ (q_{i} (t)) \ddot{q_{i}} (t) + C (q_{i} (t), \ddot{q_{i}} (t)) \ddot{q_{i}} (t) + G (q_{i} (t)) = τ_{i} (t) + d_{i} (10)$

式(10)中:qi(t)=θi(t),i=1,2,t∈[0,T]为时间变量,i∈Z+为迭代次数, $q_{i} \in R^{n}, \overset{\cdot}{q_{i}} \in R^{n}, \ddot{q_{i}} \in R^{n}$ 表示第i次迭代的关节角度、角速度、角加速度量。M(qi)∈Rn×n为惯性矩阵, $C (q_{i}, \overset{\cdot}{q_{i}}) \overset{\cdot}{q_{i}} \in R^{n}$ 表示科氏力和离心力量,G(qi)∈Rn为重力项,τi(t)∈Rn为作用在关节位置处的控制力或力矩,di∈Rn为未建模动态和扰动。

2.2 收敛性分析

控制目的是∀t∈[0,T]及∀i∈Z+,设计有界的自适应迭代学习控制器τi(t)来保证 $q_{i} \in R^{n}, \overset{\cdot}{q_{i}} \in R^{n}, \overset{\cdot \cdot}{q_{i}} \in R^{n}$ 有界,且当i→∞时, $q_{i} (t), \overset{\cdot}{q_{i}} (t)$ 收敛到理想的参考位置 $q_{d} (t), \overset{\cdot}{q_{d}} (t)$ 。为实现控制目的,给出以下基本假设:

(A1) qd(t)为可实现的关节参考轨迹角度。

(A2) 对于∀t∈[0,T]及 $\forall i \in Ζ_{+} ‚ q_{i} (t), \overset{\cdot}{q_{i}} (t), \overset{\cdot \cdot}{q_{i}} (t) ‚ d_{i} (t)$ 有界。

(A3) ∀i∈Z+,初始条件满足:

$\overset{\cdot}{q_{d}} (0) - \overset{\cdot}{q_{i}} (0) = q_{d} (0) - q_{i} (0) = 0$ 。

且满足如下4个基本特征:

(B1) M(qi)∈Rn×n为有界正定对称矩阵;

$(B_{2}) \dot{Μ} (q_{i}) - 2 C (q_{i}, \overset{\cdot}{q_{i}}) \in R^{n \times n}$ 为对称矩阵,且 $ω^{Τ} (\dot{Μ} (q_{i}) - 2 C (q_{i}, \overset{\cdot}{q_{i}})) ω = 0, \forall ω \in R^{n}$ ;

$(B_{3}) G (q_{i}) + C (q_{i}, \overset{\cdot}{q_{i}}) \overset{\cdot}{q_{d}} (t) = ψ (q_{i}, \overset{\cdot}{q_{i}}) ξ^{Τ} (t), ψ (q_{i}, \overset{\cdot}{q_{i}}) \in R^{n \times (m - 1)}$ 为已知矩阵,ξ(t)∈Rm-1为未知向量;

$(B_{4}) \forall q_{i}, \overset{\cdot}{q_{i}} ‚ t \in [0, Τ] ‚ i \in Ζ_{+} ‚ | | C (q_{i}, \overset{\cdot}{q_{i}}) | | \leq k_{c} | | \overset{\cdot}{q_{i}} | | ‚ | | G (q_{i}) | | \leq k_{g}, k_{c}$ 、kg均为正实数。

定理定义关节位置跟踪误差和速度误差分别为: $e_{i} (t) = q_{d} (t) - q_{i} (t), \overset{\cdot}{e_{i}} (t) = \dot{q}_{d} (t) - \dot{q}_{i} (t)$ ,采用控制律

${\begin{cases} τ_{i} (t) = Κ_{Ρ} e_{i} (t) + Κ_{D} \overset{\cdot}{e_{i}} (t) + φ (q_{i}, \overset{\cdot}{q_{i}}, \overset{\cdot}{e_{i}} (t)) \overset{⌢}{θ_{i}} (t) (11 a) \\ \overset{⌢}{θ_{i}} (t) = \overset{⌢}{θ}_{i - 1} (t) + Γ φ^{Τ} (q_{i}, \overset{\cdot}{q_{i}}, \overset{\cdot}{e_{i}} (t)) \overset{\cdot}{e_{i}} (t) (11 b) \end{cases}$

式(11)中: $\overset{⌢}{θ}_{i - 1} (0) = 0 ‚ φ (q_{i}, \overset{\cdot}{q_{i}}, \overset{\cdot}{e_{i}}) \in R^{n \times n}$ ,且 $φ (q_{i}, \overset{\cdot}{q_{i}}, \overset{\cdot}{e_{i}}) = [ψ (q_{i}, \dot{q}_{i}) sgn (\overset{\cdot}{e_{i}})]$ ,矩阵KP、KD、Γ均为对称正定矩阵,则ei(t)、 $\overset{\cdot}{e_{i}} (t)$ 有界,且

$\lim_{i \to \infty} e_{i} (t) = \lim_{i \to \infty} \overset{\cdot}{e_{i}} (t) = 0, \forall t \in [0, Τ]$ 。

证明第i次迭代时,构造Lyapunov函数

$W_{i} (t) = V_{i} (t) + \frac{1}{2} \int_{0}^{t} \overset{⌢}{θ_{i}} (t) Γ^{- 1} \overset{⌢}{θ_{i}} (t) d τ (12)$

式(12)中,设定 $\overset{⌢}{θ}_{i} (t) = θ (t) - \overset{⌢}{θ}_{i} (t) ‚ θ (t) = [ξ^{Τ} (t) β]^{Τ}, \hat{θ}_{i} (t) = [\hat{ξ}_{i}^{Τ} (t) \hat{β}_{i} (t)]^{Τ}$ 为θ(t)的估计值,由假设(A1)、和基本特征(B1)知,

$| | Μ (q_{i}) \overset{\cdot \cdot}{q}_{d} - d_{i} | | \leq β ‚ β > 0$

$\begin{array}{l} V_{i} (e_{\dot{i}} (t), e_{i} (t)) = \frac{1}{2} e_{i}^{\dot{Τ}} (t) Μ (q_{i}) e_{\dot{i}} (t) + \\ \frac{1}{2} e_{i} (t)^{Τ} Κ_{p} e_{i} (t) (13) \end{array}$

由于 $\bar{θ_{i}} = - θ + \overset{⌢}{θ}_{i} - \overset{⌢}{θ}_{i + 1} + θ = - \overset{⌢}{θ}_{i} + \overset{⌢}{θ}_{i + 1}$ ,

即 $\overset{⌢}{θ}_{i + 1} = \bar{θ}_{i} + \overset{⌢}{θ}_{i}$ ,则

$\begin{array}{l} \tilde{θ}_{i} (t) Γ^{- 1} \tilde{θ}_{i} (t) - \tilde{θ}_{i - 1} (t) Γ^{- 1} \tilde{θ}_{i - 1} (t) = \\ - 2 \bar{θ}_{i}^{Τ} (t) Γ^{- 1} \tilde{θ}_{i} (t) - \bar{θ}_{i}^{Τ} (t) Γ^{- 1} \bar{θ}_{i} (t) (14) \end{array}$

2.2.1 Wi(t)的非递增性证明

$Δ W_{i} = W_{i} - W_{i - 1} = V_{i} - V_{i - 1} - \frac{1}{2} \int_{0}^{t} (\bar{θ}_{i}^{Τ} Γ^{- 1} \bar{θ}_{i} +$

$2 \bar{θ}_{i - 1}^{Τ} (t) Γ^{- 1} \bar{θ}_{i - 1} (t)) d τ$ (15)

式(15)中 $\bar{θ}_{i} = \overset{⌢}{θ}_{i} - \overset{⌢}{θ}_{i - 1}$ 。由于 $\int_{0}^{t} \dot{V} (t) d τ = V_{i} (t) - V_{i} (0)$ ,即 $V_{i} (t) = V_{i} (0) + \int_{0}^{t} \dot{V} (t) d τ$ ,

又由于

${\begin{cases} V_{\dot{i}} (t) = e_{i}^{\dot{Τ}} (t) Μ (q_{i}) e_{\overset{\cdot \cdot}{i}} (t) + \frac{1}{2} e_{i}^{\dot{Τ}} (t) \times \\ \dot{Μ} (q_{i}) e_{\dot{i}} (t) + e_{i}^{\dot{Τ}} (t) Κ_{p} e_{i} \\ V_{i} (e_{\dot{i}} (t), e_{i} (t)) = \\ V_{i} (e_{\dot{i}} (0), e_{i} (0)) + \\ \int_{0}^{t} [e_{i}^{\dot{Τ}} Μ (q_{i}) e_{\overset{\cdot \cdot}{i}} + \frac{1}{2} e_{i}^{\dot{Τ}} \dot{Μ} (q_{i}) e_{\dot{i}} + e_{i}^{\dot{Τ}} Κ_{p} e_{\dot{i}}] d τ \end{cases} (16)$

由式(10)和特性(B2)和特性(B3)得

${\begin{cases} e_{i}^{\dot{Τ}} Μ e_{\overset{\cdot \cdot}{i}} = e_{i}^{\dot{Τ}} Μ (q_{\overset{\cdot \cdot}{d}} - q_{\overset{\cdot \cdot}{i}}) = e_{i}^{\dot{Τ}} Μ q_{\overset{\cdot \cdot}{d}} - \\ e_{i}^{\dot{Τ}} (- C q_{\dot{i}} - G + τ_{i} + d_{i}) \\ \frac{1}{2} e_{i}^{\dot{Τ}} \dot{Μ} (q_{i}) e_{\dot{i}} = e_{i}^{\dot{Τ}} C e_{\dot{i}} = e_{i}^{\dot{Τ}} C (q_{\dot{d}} - q_{\dot{i}}) = \\ e_{i}^{\dot{Τ}} C q_{\dot{d}} - e_{i}^{\dot{Τ}} C q_{\dot{i}} \end{cases} (17)$

由已知得

$\begin{array}{l} e_{i}^{\dot{Τ}} (Μ (q_{i}) q_{\overset{\cdot \cdot}{d}} - d_{i}) \leq \\ β | | e_{i}^{\dot{Τ}} | | = e_{i}^{\dot{Τ}} β sgn (e_{\dot{i}}) ψ (q_{i}, q_{\dot{i}}) ξ^{Τ} + β sgn (e_{\dot{i}}) = \\ [ψ (q_{i}, q_{\dot{i}}) sgn (e_{\dot{i}})] [ξ^{Τ} β]^{Τ} = ϕ (q_{i}, q_{\dot{i}}, e_{\dot{i}}) θ (18) \end{array}$

则

$\begin{array}{l} V_{i} (e_{\dot{i}} (t), e_{i} (t)) = \\ V_{i} (e_{\dot{i}} (0), e_{i} (0)) + \int_{0}^{t} e_{i}^{\dot{Τ}} (Μ (q_{i}) q_{\overset{\cdot \cdot}{d}} - d_{i} + \\ C (q_{i}, q_{\dot{i}}) q_{\dot{d}} + G (q_{i}) + Κ_{Ρ} e_{i} - τ_{i}) d τ \leq V_{i} (e_{\dot{i}} (0), e_{i} (0)) + \int_{0}^{t} e_{i}^{\dot{Τ}} (ψ (q_{i}, q_{\dot{i}}) ξ^{Τ} + Κ_{Ρ} e_{i} + β sgn (e_{\dot{i}}) - \\ τ_{i}) d τ \leq V_{i} (e_{\dot{i}} (0), e_{i} (0)) + \\ \int_{0}^{t} e_{i}^{\dot{Τ}} (ϕ (q_{i}, q_{\dot{i}}, e_{i}) θ + Κ_{Ρ} e_{i} - τ_{i}) d τ (19) \end{array}$

将控制律式(11a)代入式(19),得

$\begin{array}{l} V_{i} (e_{\dot{i}} (t), e_{i} (t)) \leq V_{i} (e_{\dot{i}} (0), e_{i} (0)) + \\ \int_{0}^{t} e_{i}^{\dot{Τ}} (ϕ (q_{i}, q_{\dot{i}}, e_{i}) \tilde{θ}_{i} - Κ_{D} e_{\dot{i}}) d τ (20) \end{array}$

据式(11b)得: $\bar{θ}_{i}^{Τ} (t) = (Γ ϕ^{Τ} e_{\dot{i}})^{Τ} = e_{i}^{\dot{Τ}} ϕ Γ$ ,有

由假设(A2)得 $V_{i} (\overset{\cdot}{e_{i}} (0), e_{i} (0)) = 0$ ,将式(21)代入式(15),有

$\begin{array}{l} Δ W_{i} = - V_{i - 1} + V_{i} - \frac{1}{2} \int_{0}^{t} (\bar{θ}_{i}^{Τ} Γ^{- 1} \bar{θ}_{i} + 2 \bar{θ}_{i}^{Τ} Γ^{- 1} \tilde{θ}_{i}) d τ \leq \\ - V_{i - 1} + \int_{0}^{t} e_{i}^{\dot{Τ}} (ϕ \tilde{θ}_{i} - Κ_{D} e_{\dot{i}}) d τ - \frac{1}{2} \int_{0}^{t} (e_{i}^{\dot{Τ}} ϕ Γ ϕ^{Τ} e_{\dot{i}} + 2 e_{i}^{\dot{Τ}} ϕ \tilde{θ}_{i}) d τ \leq - V_{i - 1} - \frac{1}{2} \int_{0}^{t} e_{i}^{\dot{Τ}} (ϕ Γ ϕ^{Τ} + \\ 2 Κ_{D}) e_{\dot{i}} d τ \leq 0 (22) \end{array}$

因Vi-1、Γ、KD均为正定阵,ΔWi≤0,因此Wi为非递增数列,得结论如下:若W0有界,那么Wi必定有界。

2.2.2 W0(t)的连续有界性证明

由式(12)和式(20)得

$\begin{array}{l} W_{\dot{0}} (t) \leq e_{0}^{\dot{Τ}} (ϕ (q_{0}, q_{\dot{0}}, q_{\overset{\cdot \cdot}{0}}) \tilde{θ}_{0} - Κ_{D} e_{\dot{0}}) + \\ \frac{1}{2} \tilde{θ}_{0}^{Τ} Γ^{- 1} \tilde{θ}_{0} (23) \end{array}$

由于 $\overset{⌢}{θ}_{- 1} (t) = 0$ ,根据控制律式(11b)得

$\overset{⌢}{θ}_{0} (t) = Γ φ^{Τ} (q_{0}, \overset{\cdot}{q_{0}}, \overset{\cdot}{e_{0}}) \overset{\cdot}{e_{0}} (t)$ ,有

$\begin{array}{l} W_{\dot{0}} (t) \leq - e_{0}^{\dot{Τ}} Κ_{D} e_{\dot{0}} + (\hat{θ}_{0}^{Τ} + \frac{1}{2} \tilde{θ}_{0}^{Τ}) Γ^{- 1} \tilde{θ}_{0} = \\ - e_{0}^{\dot{Τ}} Κ_{D} e_{\dot{0}} - \frac{1}{2} \tilde{θ}_{0}^{Τ} Γ^{- 1} \tilde{θ}_{0} + θ^{Τ} Γ^{- 1} \tilde{θ}_{0} (24) \end{array}$

因为 $θ^{Τ} Γ^{- 1} \overset{⌢}{θ}_{0} \leq Κ | | Γ^{- 1} \overset{⌢}{θ}_{0} | |^{2} + \frac{1}{4 Κ} | | θ | |^{2}, Κ > 0$ ,可得

$\overset{\cdot}{W_{0}} (t) \leq - β_{1} | | \overset{\cdot}{e_{0}} | |^{2} - β_{2} | | \overset{⌢}{θ}_{0} | |^{2} + \frac{1}{4 Κ} | | θ | |^{2} (25)$

$\begin{array}{l} 式 (25) 中 ∶ β_{1} = λ_{\min} (Κ_{D}), \\ β_{2} = \frac{1}{2} λ_{\min} (Γ^{- 1}) - Κ λ_{\max}^{2} (Γ^{- 1}), Κ \leq \frac{λ_{\min} (Γ^{- 1})}{2 λ_{\max}^{2} (Γ^{- 1})} ‚ 因此有 \end{array}$

$\overset{\cdot}{W_{0}} (t) \leq \frac{1}{4 Κ} | | θ | |^{2} ‚ \forall t \in [0, Τ] (26)$

因为θ(t)连续有界,推出W0(t)亦连续有界。

2.2.3 Wi(t)的连续有界性证明

Wi(t)可表示为 $W_{i} (t) = W_{0} + \sum_{j = 1}^{i} Δ W_{j}$ 。

则由式(15)可得

$\begin{array}{l} W_{i} \leq W_{0} - \sum_{j = 1}^{i} V_{j - 1} \leq W_{0} - \sum_{j = 1}^{i} V_{j - 1} \leq W_{0} - \\ \frac{1}{2} \sum_{j = 1}^{i} e_{j - 1}^{Τ} Κ_{Ρ} e_{j - 1} - \frac{1}{2} \sum_{j = 1}^{i} e_{j - 1}^{\dot{Τ}} Κ_{Ρ} e_{\overset{\cdot}{j - 1}} (27) \end{array}$

$\begin{array}{l} 那么 (\sum_{j = 1}^{i} e_{j - 1}^{Τ} Κ_{Ρ} e_{j - 1} + \sum_{j = 1}^{i} e_{j - 1}^{\dot{Τ}} Κ_{Ρ} e_{\overset{\cdot}{j - 1}}) \leq \\ 2 (W_{0} - W_{i}) \leq 2 W_{0}, 即 W_{i} (t) 有界。可得结论 ∶ \end{array}$

$\lim_{i \to \infty} e_{i} (t) = \lim_{i \to \infty} \overset{\cdot}{e_{i}} (t) = 0, \forall t \in [0, Τ]$ ,定理得证。

3 仿真研究

对于式(10)所示的装载机工作装置的动力学方程及式(11)所设计的控制器,仿真研究如下。设置动臂、铲斗质量分别为m1=1 793.080 7 kg、m2=2 706.447 8 kg;动臂、铲斗长度l1=2.668 2 m、l2=1.262 1 m。对式(10)描述的系统,可写为

$[\begin{matrix} m_{11} & m_{12} \\ m_{21} & m_{22} \end{matrix}] [\begin{array}{l} \overset{\cdot \cdot}{θ_{1}} \\ \overset{\cdot \cdot}{θ_{2}} \end{array}] + [\begin{matrix} c_{11} & c_{12} \\ c_{21} & c_{22} \end{matrix}] [\begin{array}{l} \overset{\cdot}{θ_{1}} \\ \overset{\cdot}{θ_{2}} \end{array}] + [\begin{matrix} G_{1} \\ G_{2} \end{matrix}] = [\begin{matrix} τ_{1} \\ τ_{2} \end{matrix}] (28)$

式(28)中, $m_{11} = \frac{m_{1} l_{1}^{2}}{3} + \frac{m_{2} l_{2}^{2}}{3} + m_{2} l_{1}^{2} + m_{2} l_{1} l_{2} c_{2}$ ; $m_{12} = \frac{m_{2} l_{2}^{2}}{3} + \frac{m_{2} l_{1} l_{2} c_{2}}{2}$ ; $m_{21} = \frac{m_{2} l_{2}^{2}}{3} + \frac{m_{2} l_{1} l_{2} c_{2}}{2};$

在matlab环境下编制m函数[8],创建装载机工作装置的自适应迭代学习simulink主程序框图如图2所示。其中,ctrl为控制器子程序,adapt为自适应律子程序,plant为被控对象子程序,input为指令输入子程序。

设定式(11)的PD控制器参数KP=KD=diag[100,100],自适应律参数Γ=diag[150,150,150,150,150],迭代次数为5,仿真时间设定为1 s。干扰项di(t)=[dmsint dmsint],dm为幅值为1的随机信号。对装载机工作装置动臂、铲斗分别施加q1d=sin(2πt),q2d=cos(2πt)的关节位置指令信号,取被控对象初始状态为x(0)=[0 2π 1 0]T,得到动臂、铲斗的跟踪控制仿真结果如图3—图6所示。

图3显示5次迭代时动臂和铲斗的关节位置追踪情况,期望值与迭代值之间存在误差,且随着迭代次数增加,误差越来越小;图4显示5次迭代后动臂和铲斗有极好的速度追踪效果;

图5、图6分别显示了动臂和铲斗的关节位置跟踪误差和速度跟踪误差图,位置和速度跟踪误差随迭代次数逐渐减小。

4 结论

在运用牛顿-欧拉法对装载机工作装置进行动力学分析的基础上,针对其不确定性和外部干扰性, 采用自适应迭代学习控制策略对其进行轨迹跟踪仿真实验。仿真结果表明,该控制器能保证系统在一定时间内稳定地减小跟踪误差,追踪性能良好,充分验证了该算法的有效性和可行性。

摘要：在理想情况下用牛顿-欧拉法建立装载机工作装置的二自由度动力学模型,考虑实际情况中存在的建模误差和扰动干扰,设计自适应迭代学习控制器对其进行轨迹跟踪控制研究,基于Lyapunov函数证明了跟踪误差的稳定性和收敛性,在matlab环境下对所设计的控制器进行仿真研究。仿真结果表明,所设计的控制器对装载机工作装置有极好追踪效果,验证了该算法的可行性与有效性。

关键词：装载机,工作装置,牛顿-欧拉方程,自适应迭代学习控制,轨迹跟踪

参考文献

[1]柳波,黄杰,等.滑移装载机自动铲掘轨迹控制研究.武汉理工大学学报,2010;32(23):107—115

[2]龚捷,鲍金峰,等.基于计算力矩法的装载机工作装置轨迹控制.机械工程学报,2010;46(13):141—146

[3]刘松国.六自由度串联机器人运动优化与轨迹跟踪控制研究.浙江大学学位论文,2009

[4]刘美华,常文森,等.机械手的显式递推牛顿-欧拉动力学及其分散自适应控制的并行方案.国防科技大学学报,1988;10(4):45—55

[5]田志祥,吴洪涛.闭环双臂空间机器人的动力学建模与仿真.华南理工大学学报(自然科学版),2011;39(8):42—47

[6]贠超,等译.机器人学导论.北京:机械工业出版社,2006

[7]张继元.机构分析与综合的解.北京:人民交通出版社,2007

自相似网络流量模型研究篇8

1994年Leland对Bellcore局域网的测试与分析成果问世以后,大量的业务流(如WAN、LAN、VBR及ISDN等)监测和分析相继表明,计算机网络上的各种业务,均呈现了统计自相似性(长相关性),即网络流量的时间序列存在着突发性。而作为计算机网络基础理论研究的前沿热点问题之一,网络流量统计分析、网络流量建模及网络性能评价也是现代通信网络规划与设计的基础。而与其相适应,基于自相似业务流的数学建模和排队分析已经成为当前网络性能评价和优化、流量控制和网络构建过程中不可或缺的方案实现要素,并且对网络规划、网络控制以及高质量的网络服务等方面的优势设置也有着重要的理论意义与现实应用价值。各种具有突发特性的业务源呈现出的自相似特性显著影响到网络的传输性能和流量控制策略,例如对时延、丢包率、吞吐量等网络性能指标的直接影响,正使得网络的设计、控制、分析和管理变得复杂。因而,只有对自相似流量下的网络性能进行正确的分析与评价,才能降低流量自相似性所带来的不利影响,使网络性能得到优化。另外,为了能够给丰富的新型业务提供强力支持,在对网络节点设备系统进行设计时,基于网络流量特性的有效性能评价将为整个系统性能指标的设计提供科学合理的计算方法。对自相似流量下的网络性能开展探索研究,则显得至关重要。而这也是本文的研究目标。下面将展开详细的分析与论述。

1 自相似基本理论

分形和自相似 (Self-Similar)[1]的概念最早形成于上世纪中期,源起于美籍法国数学家Mandelbrot对诸如海岸线长度,流体中的湍流、对流等非线性问题的研究。具有自相似规律的不规则事物称为分形(Fractals)。依据分形的自相似特性,分形主要有三类:由迭代函数系统定义出的精确自相似分形;由递推关系式定义出的半自相似分形;由不同尺度下保持统计测度的特性定义的统计自相似分形。三种分形约束依次递减。统计自相似最弱,是对自然分形对象进行定义分析的最基本约束。本文讨论的自相似就是统计自相似,又称为随机自相似。分形的度量称为分形维数D,主要有豪斯多夫维数(Hausdorff Dimension)、计盒维数(Box Dimension)和分配维数等定义方式,描述了分形空间特征。

作为分形的基本特性,自相似指的是复杂系统的整体与部分,一部分与其他部分之间在精细结构或性质上所具有的相似性。自相似具有伸缩对称性,即线性或非线性变换下的不变性,对分形对象进行放缩或者剪切等操作,只能改变其外部表现形式,而表征自相似特性的参数即分形维度则不会有任何变化。可以是在几何结构与形态、过程、信息、功能、属性和成分等表现形式上,可以是在时间、空间和数量等测度上,也可以是随机的、统计的、复杂的,但绝不仅仅是简单的按比例缩放后的重合。自相似随机过程是平稳过程。自相似性的数学表示为:

f (λr)= λαf(r)或者f(r)～rα (1)

其中,λ称为标度因子(缩放因子),α称为标度指数(分形维数)。函数f(r)是面积、体积、质量、流等属性的一种可测测度。

时间序列在统计意义上的自相似(Self-Simalry,SS),可定义为:{X(t)}满足式(2),则称{X(t)}是统计自相似的过程。

X(t)～|λ|-HX(λt),t∈(-∞,+∞) (2)

其中,～表示统计意义上的同概率分布,λ是缩放比例,H 是赫斯特指数/系数(Hurst Exponent/ Coefficient)。

随机过程{X(t), t≥0}具有长相关性[2,3,4](Long Range Dependency,LRD) 或称长记忆性(Long Memory),如果满足:{X(t), t≥0}的自相关函数R(n)=E[X(t)X(t+n)]存在,且存在常数C和α,0<α<1,使得满足R(n)～Cnα,或者有 $\lim_{n \to \infty} \frac{R (n)}{C n^{- α}} = 1$ ,同时有级数 $\sum_{n = 1}^{\infty} R (n)$ 发散,即 $\sum_{n = 1}^{\infty} R (n) = \infty$ 。常常也可用 $\sum_{n = 1}^{\infty} E [X (1) (X (n + 1) - X (n))] = \infty$ 进行判定。

有些自相似过程具有长相关性,LRD反映了自相似过程中的持续现象,意味着未来的统计信息蕴含在过去和现在的信息之中,这种信息可以通过预测和估计实现和获得,实际网络流量业务的到达就是长相关的。在要求不是很严格的情况下(默认1/2<H<1),自相似和长相关可以用来描述同一个随机时间序列,只是两者的侧重点不同,SS侧重序列的数学相似性,LRD侧重序列的统计相关。LRD用来描述时空序列自相关函数的幂定律(Power Law)衰减特性,特性表现是慢于指数式衰减。其中,幂定律表征标度不变性(Scale Invariance),如给定函数f(x)=axk,对自变量的标度缩放常数因子λ,只是造成应变量按幂次比例λk缩放,即有f(λx)=a(λx)k=λkf(x)∝f(x)。与之相对,短相关(Short Range Dependency,SRD)或称无记忆性(Memoryless),描述的就是自相关函数的指数式衰减,意味着未来的统计信息只蕴含现在的信息之中,而和过去无关,实则就是马尔科夫性。一般情况下,采用自回归分数求和滑动平均过程离散式模型和分形布朗运动连续式模型来描述随机过程的长相关性。

赫斯特系数[5] 用于描述长相关时间序列的自相关性(Autocorrelation),用于表征时间序列是回归、平均还是聚集等相关趋势,而在分形理论中则表征分形的随机程度(Randomness)并直接依赖于分形维度,有H=2-D。H的取值如下:

(1)0< H <1/2 表示负相关,表征时间序列中的一个高值之后是一个低值的可能性很高,低值之后又可能变为高值,这种高低值的交替的趋势很可能持续一段时间;

(2)1/2< H <1 为正相关,表征时间序列中的一个高值之后是另外一个高值的可能性很高,这种维持高值的趋势将很可能持续一段时间;

(3)H=1/2 表征没有相关性的时间序列,但是在很小的时间间隔中可以是正相关或者负相关,自相关数绝对值服从指数式衰减,不同于正负相关的幂定律衰减。

正相关时,H 越大,自相似程度越高。H越大,分形维度越小,曲线越粗糙;反之,H越小,分形维度越大,曲线越光滑。如图1所示。

在信号系统分析中,常常给出自相似的另外一种定义,是基于时间序列的。该定义在时间序列分析和预测中会经常用到。自相似的这一定义表述如下。

给定一个广义平稳随机过程(时间序列),X={Xn,n =1,2,…},对其进行非重叠的顺序分块,分块长度是m,得到 $X_{n}^{(m)} = \sum_{i = n m - m + 1}^{n m} X_{i}, n = 1, 2, 3, \dots$ ,称为X的聚类(过程)。如果存在聚类使得,当m→∞,D[X $_{n}^{(m)}$ ]=m2H-2D[Xn],1/2<H<1,则称X是自相似系数为H的严格二阶自相似过程。m称为聚类粒度,表征X $_{n}^{(m)}$ 对X的分辨率,m越大,分辨率越高,X $_{n}^{(m)}$ 和X越接近。如果随机时间序列聚类的统计特性不随聚类粒度的变化而有所变化,则称该序列是自相似的。

由定义可以看出,随机时间序列的聚类方差按照幂定律衰减,且衰减速率是-1<2H-2<0,将其称为聚类具有慢衰减方差(Slow Decaying Variance)。和泊松过程相比,该衰减较慢,泊松过程的聚类方差有D[N $_{n}^{(m)}$ ]=m-1D[Nn]。现在已经知道方差可以描述时间序列的波动性,也就是可以粗略描述业务流的突发性。因此,可以断定自相似业务流的突发性比泊松过程突发性大。这就解释了为什么自相似过程即使在很大的时间单位下也会保持震荡,而泊松过程却趋于平滑,也从另一个角度展现了统计自相似不随时间尺度的变化而变化的独有特性。

2 网络流量模型

分析网络流量特征是探索网络的第一步,有助于根据网络特征设计合理的拥塞控制策略,也有利于完善进一步的网络仿真实验。

根据分析粒度不同,网络流量分析可分为以下几类[6]:

(1)位级(bit),主要分析网络流量的数据特征,如链路传输速率,吞吐量等;

(2)包级(packet),主要分析IP 分组的到达、延迟、乱序和丢包等;

(3)流级(traffic),主要分析流的到达过程、到达间隔及其统计特征;

(4)应用级(application),主要基于网络应用提供的服务使用情况,搜集和分析网民行为和网络经济收益,如网络视频收视、网络广告投放和网络搜索等信息,主要用于经济领域。

在科研领域,流级流量分析研究为主导。网络流量测量技术主要有两种:

(1)被动嗅探式,不注入数据包、精度高、耗费资源,诸如TcpDump、WinDump和WinPCap等;

(2)基于SNMP参数测量,主动发送数据包、精度低、且可扩展性高。

网络流量模型(Traffic Model)是网络流量统计特征的研究工具。通常,采用时间序列表示某个特定时间或时间间隔内到达的数据包数量,流量则为该时间段内的平均值。用数学语言来描述即可,在时间点上的一个观测值序列表示流量,这个观测值序列的时间尺度可以是毫秒、秒、分钟、小时等,流量的单位一般是字节每秒(B/s)。网络流量模型一直以来就倍受学者关注,各类相关模型都已很多。按照网络流量相关性特点可以分为短相关流量模型和长相关流量模型两大类。短相关流量模型的自相关函数随时间间隔的增加呈指数衰减。长相关流量模型的自相关函数则随时间间隔的增加呈幂定律收敛,而且比指数衰减要慢。

近年来,研究发现局域网和广域网的流量都呈现统计自相似特征。这种特性就决定了自相似业务流模型已经成为模拟实际网络流量的主要手段,有关这方面的研究较多,目前主要有两类[7]:物理模型和统计模型。物理模型的建立基于自相似过程的物理意义,对于自相似的成因和特点有较强的表现力,使用频度较高的物理模型是流叠加法的ON/OFF模型[8]、TES模型[9]、α-Stable业务流模型[10];统计模型的建立则基于长相关性随机过程,形式灵活,精度较高,常用的统计模型有FARIMA模型[11]、分形高斯噪声模型(Fractional Gaussian Noise FGN)[12]、基于小波变换或马拉特(Mallat)模型、基于混沌映射的确定性模型、散粒噪声模型等,甚至是多重分形[13]。这类模型能很好地说明网络通信量中出现的长相关和重尾等现象,但对于瞬时性能的评估却非常困难。除前面介绍的流量模型外,目前已在研究的模型还有,漏桶模型、瀑布模型[14]、季节性神经网络流量模型等。其中,漏桶模型比较适合于特定网络应用流量的分析。例如,VBR的多媒体模型[15]。但是,该模型难以反映网络流量实际变化中的突发特性,尤其是混合流的特性。

经过研究发现,有一随机过程BH(t),如果满足以下条件:

BH(0)=0,数学期望E[BH(t)]|T=0,协方差函数为Cov(t,s)=E[BH(t)BH(s)]=1/2(|t|2H+|s|2H-|t-s|2H),0<s<t<T,其中,T为考察周期。随机过程BH(t)为分形布朗运动(Fractional/Fractal Brownian Motion,FBM)。参数H(0,1),是赫斯特指数,用于描述分形布朗运动的分形特性,H越大,该运动轨迹越平滑。当0< H <1/2时,为负相关;1/2< H <1时,为正相关;H=1/2 时,退化为无自相似性的维纳过程。正相关时,表现出长相关性,且H越大,自相似程度越高,故H也称为自相似系数。

FBM具有自相似性,即BH (αt) ～ α2HBH (t),由于FBM的协方差函数是齐次的2H阶。

FBM具有平稳增量,即BH (t)-BH(s) ～ BH (t-s)。FBM是一个连续时间的高斯过程,且是唯一的一个具有自相似性的高斯过程。故增量过程X(t) =BH(t+1)-BH(t),t≥1可称为分形高斯噪声(Fractional Gaussian Noise,FGN)。其对应期望E[X(t)]=0,自相关函数R(n)=E[X(t) X(t+n)]=2-1[(n+1)2H-2n2H +|n-1|2H],n≥0,满足当n→∞,H≠1/2时,R(n)～H(2H-1)n2H-2,而当H=1/2时,则转化为白噪声(White Noise)。进一步可以得到,当1/2 <H<1时,分形高斯噪声X(t)正相关且具有长相关性。

通常情况下,FBM没有独立增量,故其不是Lévy和鞅的。并且,普通FBM也不是平稳过程,但标准FBM却是平稳过程。

FBM具有长相关性,这是因为,当1/2< H <1时,有 $\sum_{n = 1}^{\infty} E [B_{Η} (1) B_{Η} (n + 1) - B_{Η} (n)] = \infty$ 。

FBM样本路径几乎处处均不可微(differentiable),然而基本上所有轨迹(Trajectory)都是少于H阶的任意霍尔德连续(Hölder continuous)。对于这样的轨迹, 存在一个常数c,使得任意ε>0,都有| BH (t)-BH(s)|≤c|t-s|H-ε。BH(t)曲线的Hausdorff维数和Box维数均为2-H。

最为经典的FBM模型,当属由Norros提出的分形布朗运动模型[16]:

$A_{Η} (t) = m t + \sqrt{m a} Ζ (t), t \in (- \infty ‚ + \infty) (3)$

其中,AH(t)表示t时间内到达的业务流,下文统一记为“A(t)”;Z(t)是标准的参数为H的分形布朗运动;分形布朗运动有3个参数,分别是:m>0表示平均发送速率;a>0表示偏差系数,H∈[0.5,1)为Z(t)的Hurst系数。Z(t)均值E[Z(t)]=0、方差D[Z(t)]=| t |2H,当H=0.5时,A(t)为无自相似性的布朗运动业务流。

3 网络自相似流量

因特网业务流量由响应业务流量(即TCP长流,Long Lived TCP Flows)和非响应业务流量(即UDP流和TCP短流,Short Lived TCP Flows)构成。有研究表明,Internet流量主要由TCP短流(如Web流量,HTTP应用)构成[17];TCP短流使得平均队列深度都表现出了类指数的特性。由此可知,应将短流的队列行为构设进入网络设计时的决策范畴[18]。2008年,IResearch咨询公司通过Cisco System 提供的数据统计结果表明, P2P和在线视频服务流量已经超过了网页浏览、电子邮件这些传统的网络服务, 正在占用越来越多的网络带宽。以2008年为例, 在线视频( 包括PC 和电视终端) 和P2P流媒体产生的数据流量已经占据全球互联网总流量的75. 9%,达到4 034 PB, 而当初预计到2012 年,这一数字将会继续上升至81.2%。另外,对等网(Peer-to-Peer,P2P)技术应用呈现了空前繁荣,如文件共享、协同计算、流媒体、IP-TV、VoIP语音视频通信及在线游戏等应用的陆续出现。基于实时性考虑,这些新兴应用协议多选择UDP作为其底层的传输协议,使得UDP流量呈上升趋势[19]。有研究分析可得,大概有30%-70%的网络流量产生于这种技术。而该技术又无疑带来日益增多的非响应业务流量(主要是UDP流)。因此则势必给网络流量的稳定性带来影响,而且也必将使路由器的排队性能面临严峻挑战。同时,文献[20]中基于对Internet 城域出口链路流量的准确测量,又一次提到了网络流量的短相关性不再明确,由此对于利用控制理论和排队理论开展AQM算法的性能分析制造了一定难度。

1994年,Leland等人发现了Bellcore的局域网网络流量具有自相似特性[21],开启了自相似网络流量的研究进程。此后,Paxson[22] , Crovella[28]等人分别验证了泊松采样测量的失效,表明网络流量具有广泛的统计自相似特性(Self-similar),从而诠释了马尔科夫链和泊松过程等短相关模型已不具说明功效的网络现象。其中的广泛性是指网络流量的时间序列在不同的时间尺度(毫秒～小时)上都存在着突发性。无论网络的规模、拓扑、应用、编码、传输介质如何变化,这种突发性始终存在。从某种意义来说,网络流量的突发性就是自相似特性的具体表现。而网络流量的自相似性质的发现,则成为网络测量领域和网络行为辨识方面的一个里程碑似的重大突破。其后,多种自相似模型相继引入,更新和拓展了人们有关网络流量特征的认识。自Taqqu等人提出实际网络流量呈现复杂的多重自相似[23]的理论后,基于网络多重自相似特征的研究已经日渐涌现[24,25,26],多重分形(Multifractal)业务流可以看成是多个基于不同网络应用的单分形业务流的合成和叠加。基于时间尺度对网络流量自相似的分析过程已变得越来越复杂,文献[27]就提出一种基于行为尺度的自相似,相关实验表明,在P2P业务流中存在时间尺度上的多重自相似,虽然简单的自相似无法描述,却存在较完美的行为自相似,即基于网络应用角度的网络行为(不同数据包收发粒度)的自相似。目前,已将基于多重分形的网络流量建模和网络性能评价确定为自相似网络的研究方向。

4 网络流量自相似性成因

关于网络流量自相似的成因,当前的解释很多。大致可以归纳为两个:网络文件的重尾分布[8,23,28,29] 和TCP的拥塞控制机制[30,31]。

Willinger等人认为计算机网络文件大小服从重尾分布是导致网络流量自相似性的主要原因。首先,可以把网络中端到端的链接看做是一个ON/OFF源,其中,ON对应有数据包发送,为链接忙时间,OFF对应无数据包发送,为链接闲时间,在网络流量足够大,网络链接足够多时,可认为这些ON/OFF源是独立同分布的;又由于忙时间和闲时间的持续都可能很长,且难以忽略,即呈现“Noah”效应,从而可以用重尾分布描述这些时间间隔。其次,Web文件(或者视频流文件)大小的分布(包括用户请求的文件、实际传输的文件、服务器端存储的文件等)服从重尾分布,使得网络传输时间具有无限方差,即服从重尾分布,从而导致整个链路层上的自相似;用户的想时间(think time)的重尾分布进一步导致了网络空闲时间服从重尾分布,并且ON状态比OFF状态还要重尾。大量重尾分布的ON/OFF源的聚合就产生了自相似性。最后,可靠的传输机制和流量控制机制又保留了由文件大小重尾分布所引发的长相关性——注意无流量控制和不可靠的UDP并未使得生成的流量具有长相关性。况且对流量自相似的估计并不因网络拓扑结构变化而改变,或者说自相似网络流量经过网关路由的转发处理并不能削弱自相似性。一个自相似过程的分支仍然是自相似的,若干个自相似过程的聚合也仍然是自相似的。网络业务流一直存在自相似性,不仅不会随着业务的聚合和分支而削弱,反而自相似系数还会增大,并使得自相似统计特征变得愈加复杂,由此产生了多重自相似。

Veres等人认为,TCP的拥塞控制机制亦是流量自相似特性的可能成因。文件的可靠性传输,也就是重传机制(Retransmission),即使改变其参数,如缓存大小、重传预定的次数和超时时限,也不能改变重传负载的自相似特性;同时,也并不随着网络源、拓扑、业务流汇聚和到达间隔时间分布的变化而变化;此外,还存在着传输层流量控制机制和可靠传输对网络流量的整形。当时间尺度超过10倍的数据包传输时间,重传数据包流量的方差将在总的流量(新数据包、重传数据包和丢失的数据包)中占据绝对优势成分,这就意味着极大的突发性,从而在某种程度上使得单个的TCP流量符合渐进自相似(H>0.5)。虽然在瓶颈缓存处堆叠的TCP流量是短时相关的(H≈0.5),但TCP拥塞控制可使瓶颈缓存占用率最大来平滑流量,堆叠的流量得到平滑,并在TCP拥塞控制和具有重尾特性的上层协议的共同作用下,使得堆叠的网络流量仍然显现了长相关性。TCP拥塞控制中包括着混沌特性,诸如:非线性(Nonlinearity)、确定性(Determinism)、混乱中的有序(Order in disorder)、对初始状态的敏感性(蝴蝶效应)(Sensitivity to initial conditions or the “butterfly effect”)、不可预见性(Unpredictability)。系统在特定的参数下产生自相似时间序列,从周期性到产生明显的混沌现象:对初值敏感,流量在广泛范围内具有自相似性,甚至是多重分形。周期性是由于流量数据的周期采集而引起,或者是由人们上网的行为习惯所引起。混沌性完全从网络流量自相似的频谱角度而进行更为完全的解释。由此得出,TCP本身就是一个产生自相似结果的确定性过程。

5 自相似网络流量的仿真和生成

5.1 赫斯特指数估计

Hurst指数估算有很多方法,可分为时域和频域两类,较为常用的有重标极差(Rescaled Range,R/S)分析法[32,33]、方差时间图(Variance Time Graph,VT-G)法[34]、留数法(Variance of Residuals,Res)[35]和绝对值法(Absolute Moment,Abs)[36]等时域算法;以及Whittle法[37]、小波基(Wavelet-transform Based,WB)法[38]和周期图法(Pariodogram Graph,PG)[39]等频域算法。其中,各算法的时间复杂度分别是: R/S法、Res法和Whittle法的结果为O(N2),WB法和PG法为O(Nlog(N)),而VT-G法和Abs法则为O(N2)[36]。文献[40]对这几种常见的算法进行了分析比较,并给出了各算法的时间复杂度;同时,进一步指出VT-G法、Abs法等聚类方差法(Aggregated Variance)[41],计算速度较快,结果相近;PG法在实现时可借助快速傅利叶变换来提高算法的速度;Res法和R/S法的速度相对较慢,但结果精度则相对较高;小波法虽然速度较快,实现过程却相当复杂。当H>0. 80以后,对FGN序列的时域方法估计值将失准且偏低,但频域方法估计值却仍能保持较高精度。文献[42]中,通过一系列实验也指出,时域类的方法性能全部都要低于频域类算法。而在频域类算法中,Whittle法的精度最高。文献[43]还指出,在人工合成数据序列下,各估算算法均表现良好,但在真实数据流的情况下,算法准确度却集体下降,这可能和真实数据的多重分形特性有关;另外,在突发噪声干扰下,对时域类算法的精度影响较大,而对频域类则具有较好的鲁棒性,加入现有滤波器技术也不能有效改进算法的精确度;建议研究者们避免采用单一算法,导出相应结论。目前,国内提出的最新计算方法主要有滑窗时变方差之差法[44]、局部Whittle法[45]、基于混合FBM的二次变差矩估计法,而国外提出的最新计算方法则有最大似然估计法(Maximum Likelihood)[46]、基于小波的变方差(Time-Varying)法[47]和基于FBM的锥多元自适应回归曲线法(Conic Multivariate Adaptive Regression Splines,CMARS)[48]。综上所述,时域算法大部分是通过作图估计H,精确度普遍不高,且需要较大样本空间(一般大于10 000);频域算法则只需要很小的采样序列即可,而且精度较高。

限于篇幅,此处仅粗略描述两个算法,以供诸位研究同仁参详与考量所用。具体内容如下。

首先,介绍一下最简单的算法VT-G法。给定一个时间序列X={Xn,n =1,2,…}, VT-G法主要是利用公式D[X $_{n}^{m}$ ]～m-βD[Xn],m→∞计算聚类方差,以此描绘得到方差的log坐标系图线,其拟合直线斜率即为-β,0<β<1,进而求得H=1-β/2。有时候也可用β来描述自相似程度。

然后,是应用最为广泛的R/S法。给定一个时间序列X={Xn,n =1,2,…}计算重标极差序列如(4)所示,并计算重标极差期望E[R(n)/S(n)]。Wallis已证明当t趋于无穷大时,E[(R/S)t]～C*tH。依据幂定律拟合数据估算Hurst系数。实际操作时,可以在log-log象限中作出重标极差期望的时间曲线,其拟合直线的斜率即是H。

$\begin{array}{l} \frac{R (n)}{S (n)} = \frac{\max {0, [\sum_{i = 1}^{k} X_{i} - k m} - \min {0, [\sum_{i = 1}^{k} X_{i} - k m]}}{\sqrt{\frac{1}{n} \sum_{i = 1}^{n} (X_{i} - m)^{2}}}, \\ m = \frac{1}{n} \sum_{i = 1}^{n} X_{i} (4) \end{array}$

5.2 分形高斯噪声

在已有的流量模型中,分形布朗运动(FBM)模型是最简单、最易于求解的自相似业务流模型。由于分形高斯噪声(FGN)是FBM的增量过程,故常常使用FGN生成FBM流量[49]。通过利用分形高斯噪声合成近似分形布朗运动网络业务流的快速技术主要有[50]:随机中点置位法(Random Midpoint Displacement),连续随机添加法(Successive Random Additions)和浮动比例法(Floating Proportionality)。更多关于分形的生成算法可参见文献[51]。近年来,还有一些其他生成自相似网络业务流的方法,诸如基于一般化柯西过程(Generalized Cauchy Process)的模型[52]。

当今,采用一种基于循环嵌入法(Circulant Embedding Approach,CEA)[53]的快速傅里叶变换(Fast Fourier Transform,FFT)[54]方法来生成FBM业务流,则是一种最快的算法。Perrin在文献[55]中已经验证了CEA算法生成分形高斯噪声是一种最优方案。考察不同的H生成的FBM,可知H越大,分形维度越小,生成曲线越光滑。生成曲线效果如图1所示。

6 展望

目前,自相似网络流量研究主要有三个方面:分析网络流量的统计特征并实现建模,包括“可信的”网络流量生成模型[56,57]和“可靠的”网络流量预测模型的构建[58,59]。其中,小波基(Wavelet Based)分形理论[60,61]和多重分形(Multifractal)模型[62,63]成为难点和趋势;基于自相似网络模型评估自相似流量对不同网络的各种性能影响[64,65];网络自相似或者长相关的成因[66,67]。自相似网络诸多问题还末得到彻底解决,很多问题尚处于讨论阶段,研究成果分散繁杂,也没有形成较为一致、清晰、且完整的体系,亟需进一步投入,加大研究发展力度。到目前为止,也还没有取得一个统一的、公认的数学模型来描述自相似网络流量。

7 结束语

本文从分形和自相似理论基础引入,主要介绍了网络流量的自相似特征以及形成原因;为了更深刻地认识网络流量特征,进一步介绍了常用的网络流量模型;最后对自相似业务流量的合成在Matlab中进行了仿真实现。借助本文,相信能对网络自相似形成了一个较为系统的认识。

基于自回归的高斯混合模型篇9

关键词：高斯混合模型,自回归模型,EM算法,聚类

一、简介

聚类问题是当今统计学领域的热点研究问题之一, 统计学家已经发明了诸多成熟的算法。随着数据特别是生物中基因数据的快速增长, 人们逐渐发现数据中的一些内在的性质, 利用传统的聚类算法处理这些数据无法利用这些性质因而无法提供较高的准确率。为此, 研究人员更多的是针对一类特殊模型提出特殊的算法, 这实际上是聚类细化的问题, 例如, 作者提出了一个混合模型 (Mixture Model) 来对基因表达数据进行聚类分析, 该模型运用了在基于基因表达试验中的设计矩阵, 作者将算法运用到三种不同的数据上并取得了非常好的效果。同样在中, 作者提出了一个混合模型来对周期性的基因数据进行聚类研究。

二、模型介绍

混合模型 (Mixture Model) 是聚类算法中常用的模型, 具有较强的理论基础, 通过利用期望最大化算法, 可以迅速地求解。现在我们将混合模型的思想运用在基于自回归数据中。假设自回归数据可以分为c个类, 类符合参数为π= (π1, …, πc) 的多重分布 (Multinom ialDis tribution) , 其中类h是符合正态分布, 其期望为uh= (uh1, uh2, Luhm) 方差Sk=δk Rk, 其中

混合模型的参数为μ, ρ, δh, 其中μ代表了c个类的期望值向量, ρh代表了c个类的自回归系数 (ρ1…ρk…, ρc) 。

若数据点yi= (yi1, yi2, …, yim) , 用zhi表示yi的类指标, 这样很自然的有 (∑ch=1zhi=1) , 若zhj=1表示数据点yi是由类h来产生的, 那么它出现的概率为

数据点yi出现的概率为

若有n个数据点 (n个样本) y1, y2, …, yn, 则数据的对数似然函数为

现在我们的目标是要估计出参数μ, ρ, δ以极大化似然函数L (1) .首先考虑对式子 (1) 进行求导并设置导数为零, 但由于log函数中含有和的形式, 无法通过求导来得出结果。期望最大化算法对于处理这种情形是非常自然的, 在下一小节, 我们将介绍如何利用期望最大化算法 (EMAlgorithm) 解决上述优化问题。

如果参数μ, ρ, δh已经被正确的估计出来, 那么数据点yi来自于类h的后验概率为τhi=f (zhi=1|u, ρ, δh, yi) =

因此可以将数据点yi分配到τhi最大的那个类中, 这也是混合模型进行分类的标准。

三、实验结果

对自回归系数估计的实验

本部分将对算法来进行测试以证明有效性, 算法的一个非常重要的部分是对相关系数ρh的估计, 如果能保证这部分的正确性, 则算法的整体正确性会得到保证。为此我们首先进行如下的试验:对于ρ=0.1, 0.3, 0.7, 0.9, δ=0.3, 1, 1.5分别生成期望为0, 从而估计ρ和δ的值, 设置m=10, 每个实验进行10次, 结果在表中。

四、总结

本文提出了一个基于自回归模型的混合模型并利用期望最大化算法给出了这个算法的迭代过程, 具有较完善的理论基础。在实验部分, 通过测试多种自回归数据并对比常用的Kmeans和Mclust算法, 我们证实了算法在处理自回归数据上的高精度。算法应用在基于时间点的基因表达数据上, 我们未来的工作将集中于算法的应用并进行算法的模型选择问题研究。

参考文献

[1]Kim BR, Zhang L, Berg A, Fan J, and Wu R.A computational approach to the functional clustering of periodic gene-expression pro?les.Genetics, 180:821–834, 2008.

[2]C.Fraley and A.E.Raftery.Enhanced model-based clustering, density esti-mation, and dis-criminant analysissoftware:Mclust.J.Classif., 20 (2) :263286, 2003.

[3]Hironori Fujisawa.The maximum likelihood estimatorsin a multivariate nor-mal distribution with ar (1) covariance structure for monotone data.Annals of the Institute of Statistical Math-ematics, 48:423428, 1996.

自学习模型篇10

[关键词]购买力平价理论；实际汇率；门限自回归（TAR）模型；非线性门限单位根检验

[中图分类号]F831 [文献标识码] A [文章编号] 1673-0461（2012）10-0012-06

一、引言

购买力平价理论（Purchasing Power Parity，简称PPP）是简单且较具实用价值的解释汇率的理论，该理论的核心思想在于，如果市场完全有效，任何相同的商品在世界范围内都应该以相同的价格销售。20世纪70年代以来，针对购买力平价理论的实证分析大量涌现，尽管研究所采用的分析工具不尽相同且研究发现多种多样，从研究结论和方法的角度看来，两大特点还是可以从目前围绕该理论建立的实证研究中总结出来：①针对PPP在普遍意义上成立与否目前仍然尚无定论，国际上对PPP 理论假说的较一致意见是：就长期而言，购买力平价可能发生市场出清的均衡现象；但从短期来看，其效果并不明显。②大量研究假设数据产生的过程为线性，故采用单变量单位根检验（如ADF方法、PP方法等）对其数据序列进行平稳性检验。

以上特点的第二点中针对实际汇率这一单一时间序列作线性假设的基础本身较为薄弱。Hechscher最早于1916年就提出实际汇率的调整可能是非线性的观点，结合其对国际贸易中出现的交易成本的交易成本的分析，他认为尽管同样商品的价格由同一货币在不同国家之间进行衡量可能存在不一致，但只有当预期的收益大于进行商品贸易的两国之间的运输成本时，套利和由之产生的价格修正才会出现（Taylor，2006）。

随着计量分析技术的发展，2000年以后开始大量出现采用非线性的计量方法对实际汇率的平稳性进行研究。国外的相关研究中，较具代表性的包括Taylor（2001），Chortareas等，（2002）， Chortareas和Kapetanios（2003）， Kapetanios等（2003），Kilian和Taylor（2003）， Canjels等（2004），如上研究运用非线性方法均发现了长期购买力存在的证据，但其研究数据大多取自于发达经济体，该类研究客体在市场有效性、外部冲击的影响性等方面的不同决定了其结果对于发展中经济国家或地区具备较为有限的参考价值。而针对我国的研究中，具备较强参考意义的包括王志强等人（2004）结合界限检验的方法发现1994年汇率制度改革至2004年，人民币的长期汇率购买力平价得到部分经验证据支持；而刘金全等人（2006）利用不同的方法和样本得出了不同的结论，他们通过马尔科夫区制转换的方法对1980年1月至2004年8月的人民币汇率进行了分析，认为长期购买力平价并不成立；此外，张卫平（2007）则运用平滑转换自回归模型（STAR）对人民币实际汇率的非线性行为进行了实证分析，并对其均值回复速度进行了估计，研究结果特出人民币与美元存在长期购买力平价关系；最新的研究包括丁剑平等人（2010）运用门限单位根检验方法研究亚洲6个主要新兴市场国家货币兑美元实际汇率，结果表明研究对象的实际汇率均存在门限非线性。

就研究结论而言，如上中外研究均采用非线性方法发现了PPP存在与否的经验事证，其结果由于样本、方法的不同而存在较为显著的差异；但从研究方法上看，这些研究存在的共同问题是无法在非线性和非稳态之间进行区分，故仍采用线性研究中广泛使用的传统检验方法对非线性关系的平稳性进行检验。考虑到如上问题的普遍性，本文采用最新发展起来的由Caner和Hansen提出的计量分析和检验方法，旨在为针对购买力平价理论的实证研究注人新的活力。本文的主要贡献在于通过有针对性地采用专门针对非线性关系的单位根检验，对真实汇率这一时间序列数据的长期稳定性特质进行全面与透彻的分析，以更好地发现购买力平价理论在长期的表现。

此外，从更为深遠的研究意义上考虑，组成一个共同货币区的必要（非充分）条件是影响这些国家（地区）或经济体的外部冲击因素。而如果通过一般购买力平价理论检验发现，这些国家（地区）或经济体的实际汇率存在趋同，则可证明这些国家或地区是潜在的最优货币区。从这一视角来看，一般购买力平价理论与Mundell提出的最优货币区理论存在共同之处。在此基础上，本文选择人民币兑换两岸三地（即中国香港、中国澳门和中国台湾）的实际汇率作为研究对象，可为未来针对中国大陆是否适宜于与中国香港、中国澳门、中国台湾组成同一货币区的相关研究从汇率相关性角度提供一定的研究事证和基础。

本文其他部分的安排如下：第二部分给出本文的模型介绍、设立方法及其相应的检验方法；第三部分解释相关数据的选择依据和处理方法；第四部分简明介绍相关分析结果数据并给出实证分析结果；第五部分结合前述研究成果提出总结性评述。

二、研究方法

1.模型介绍

实证分析中常见的非线性机制转换模型分为三个大类：即马尔科夫区制转换模型（Markov Switching Regime model）、门限回归模型（Threshold Regression Model）和平滑转换自回归模型（ Smooth Transition Autoregressive Model）（封福育，2011）。本文的实证研究模型的基础属于门限自回归模型（Threshold Autoregression Model， TAR）的一种，最早由Tong（1978）提出，并在非线性时间序列的实证领域得到广泛应用。

在Tong（1978）TAR模型核心建模思想的基础上，Cancer和Hansen（2001）通过考察回归单位根的渐近零分布，进一步发现该分布是非标准的，且其标准与否与门限效应是否存在直接相关，并由此发展其非限制性两区制门限自回归模型（Unrestricted Two-Regime Threshold Autoregressive Model）及其相应的检测方法。对比其他门限自回归方法，Cancer和Hansen的方法的最主要优点在于将非线性问题（即门限效应）和非平稳问题（即单位根）问题分开考虑，并通过Wald检验和非线性单位根检验分别检测时间序列的非线性和平稳性（P.K.Narayon and S.Narayon，2007）。从实证检验的结果上看：Cancer 和Hansen （2001）通过对美国失业率的平稳性考察认为，在检测非线性时间序列的平稳性问题上，其非线性单位根检测的方法要优于传统的ADF、PP等平稳性检验方法；S.H.Sekioua （2006）在实证研究远期汇率溢价的非线性调整问题的过程中，佐证了Caner和Hansen的观点；此外，P.K.Narayan（2005），J.Madsen，V.Mishra和R.Smyth（2008）等人针对不同研究对象的比较性检验结果均支持这一观点。由此，我们认为通过采用Cancer和Hansen模型和其相应的非线性单位根检验方法对两岸三地实际汇率的非线性和平稳性进行研究将得到更具参考意义的研究过程和更为稳健的实验结果。

2.模型设定

结合Caner和Hansen（2001）的文献与本文的研究对象，我们得到如下两区制TAR（k）形式：

？驻St= ？兹'1Xt-1I{Z■？燮？姿}+？兹'2Xt-1I{Z■>？姿}+？着t，（1）

其中，Xt-1=（St-1，1，？驻St-1，…，？驻St-k）′，（2）

如上模型中，

（1）St代表被检验的我国主要贸易国家（地区）的实际汇率，（其中t=1，…，T，用直接法以美元作为单位）；

（2）？着t是满足独立分布的随机误差项；

（3）I{.}代表指标函数，（当Zt-1？燮？姿时，I{.}=1，相反，当Zt-1>？姿时，I{.}=0）；

（4）k代表自回归阶数，满足k？叟1；

（5）？姿是阈值参数，Zt-1=St-1-St-m是阈值变量，其中m？叟1（在公式1中，？姿作为重要的未知变量，可使St随阈值变量Zt-1=St-1-St-m在两种区制间进行转化）；

（6）向量？兹1和？兹2表达为：

？兹1■■ ？兹2■■ （3）

公式3中，？籽1和？籽2为标量，代表St-1的斜率系数，？茁1和？茁2是以标量形式存在的截距项，且二者的值同为1，其代表确定性部分的斜率，？琢1和？琢2为1×k向量，其代表（？驻St-1，…，？驻St-k）在两个区制中的斜率系数。

（7）阈值的估计方法

为得出最终估计值，本模型用最小二乘法对m进行多次估计。在此过程中，未知变量？姿？缀[？姿1，？姿2]、？姿1、？姿2服从以下条件：

Pr（Zt？燮？姿1）=π1，Pr（Zt？燮？姿2）=π2（4）

其中，0<π1<π2<1，且π1+π2=1。

向量？兹1和？兹2中每一个参数的估计值均可通过最小化如下公式得出：

Q（？姿，m）=■？着t（？姿，m），（其中，t=1，…，T），（5）

给定？姿和m并使用最小二乘法，得到残差的估计值■t （？姿，m）由此，阈值变量的估计值可通过下式得到：

■■（■，m）。（6）

3.阈值检验

利用标准的Wald统计量，可对阈值进行检验：

WT=WT（■）=sup？姿？缀 [？姿■，？姿■]WT（？姿），（7）

Wald检验的零假设为：

H0∶？兹1=？兹2，（8）

如果不拒绝零假设，即代表实际汇率序列没有门限效应。实际汇率服从线性变化。

如拒绝零假设，则代表实际汇率序列服从非线性门限模型。

同时，sup？姿？缀[？姿■，？姿■]WT（？姿）服从非标准分布，Caner和Hansen（2004）建议采用自抽样法（Bootstrap Method）以得到渐近临界值和P值。

4.门限单位根检验

检验实际汇率变化的平稳性可应用完全单位根检验法和部分单位根检验法。

（1）假设

如果两个区制包括一个单位根，零假设为：

H0∶？籽1=？籽2=0，（9）

备择假设为：

H1∶？籽1<且？籽2<0，（10）

如接受零假设H0，则代表两个区制均包含一个单位根，实际汇率序列在两个区制中均不遵循一个平稳的过程；如拒绝零假设H0，则代表实际汇率序列服从平稳过程，即长期购买力平价成立。

如某一区制下存在单位根，备择假设可被重新构造为：

H2∶？籽1<0且？籽2<0，或？籽1=0且？籽2<0，（11）

如接受备择假设H2，则存在一个区制，在该区制下存在一个单位根；而另一区制下存在平稳性和均值回归。

（2）检验

门限单位根检验通过构造两个Wald统计量进行，即：

S1T=t■■I■+t■■I■，（12）

S2T=t■■+t■■，（13）

其中，S1T是單侧检验统计量，与备择假设？籽1<0或？籽2<0相对立；而S2T是双侧检验统计量，与备择假设？籽1≠0或？籽2≠0相对立；t1与t2则分别是估计值■1和■2的比率。Caner & Hansen（2001）指出，单侧Wald检验统计量S1T比双侧Wald检验统计量S2T更强，所以我们在如下的实证分析中主要采用单侧检验单位根的存在性。

三、数据选择

本文运用MATLAB 7.11.0进行实证检验。利用中国经济数据库（CEIC）选取1995年1月至2011年12月中国大陆、中国香港、中国澳门和中国台湾的月度名义汇率和消费者价格指数数据，共计1，428个样本值进行实证分析。实际汇率的得出遵循如下满足相对购买力平价理论的公式，通过该公式可将名义汇率转化为实际汇率。

St=■，（14）

其中，St表示实际汇率；S■■表示名义汇率；P■■表示外国物价水平；P■■表示我国物价水平。如上公式14中各变量的数据整理方法如下：

（1）对于名义汇率S■■，本文采用人民币直接标价法，选取港币、澳门币、台币兑人民币的月末汇率数据。

（2）对于一国价格水平P■■和P■■，本文选取各国CPI作为该指标的代理变量。由于消费者价格指数（CPI）作为衡量一国通货膨胀水平的常用指标包含一篮子商品和服务的价格，因此，本文选取CPI作为衡量一国价格水平的指标具备理论与现实意义。此外，为消除样本选取区间对分析的影响，本文以1995年为基年对所有CPI指数进行了调整（即令CPI1995=100）。

四、估计结果与实证分析

1.估计结果

用最小二乘法对TAR模型的估计结果如下表1所示。三组时间序列中，每一序列滞后12期，包含191个观测值，自抽样次数为300。

2.Wald检验

在估计结果的基础上，我们通过如上所述的Wald检验检测人民币兑港币、人民币兑澳门元、人民币兑台币的实际汇率是否存在门限效应。Wald检验结果如表2所示。

首先，最优延迟期数m来自Zt-1=St-1-St-m，由于本实证检验的最大滞后期数设置为12，故m的取值范围在1到12之间。最优延迟期数是使残差的方差最小的延迟期间（在该期相应的Wald统计量WT取值最大）。港币序列、澳门元序列、台币序列在最优延迟期数m分别等于4、11、6时，Wald统计量取值最大，分别为46.42、47.23和31.32。

其次，就Wald检验的结果而言，人民币兑港币、澳门元和台币的实际汇率序列在不同自抽样显著性水平下的临界值与Wald统计量的比较结果存在差异，其中，人民币兑港币的实际汇率在显著性水平为1%、5%和10%下的临界值均大于Wald统计量，门限效应显著存在；人民币兑澳门元的实际汇率在的显著性水平为10%时的临界值（45.73）小于Wald统计量（47.23），零假设不能被拒绝，同样被认为存在门限效应；对台币而言，在所有被考察的显著性水平下临界值均大于Wald统计量，零假设被接受，即可认为人民币兑换台币的实际汇率不存在门限效应。

再次，本次实证中的其他统计参数包括：自抽样P值分别为0.03、0.09、0.50，单位根P指分别为0.02、0.07、0.47，而门限估计值■分别为-0.01、0.01、0.01。

最后，港币、澳门元、台币的实际汇率落在区制1中的数量和比率分别为87（45.55%）、161（84.29%）和148（77.49%），该结果的折线图可由图1至图3直观表示。

总结Wald检验的结果可以发现，人民币兑港币和人民币兑澳门元的实际汇率呈非线性趋势；而从长期来看，人民币兑换台币的实际汇率则更多的表现出线性特征。针对如上线性和非线性的平稳性检验，应采取不同的检验方法以得出稳健的结论。

3.针对非线性结果的单位根检验

对呈现非线性特征的人民币兑港元实际汇率和人民币兑澳门元的实际汇率应用单侧单位根检验，结果如表3所示。

结果发现，港币和澳门元的Wald统计量S1T分别为0.10和8.13，均小于统一最优延迟期数下，不同显著性水平下的临界值（港币-9.43（10%）、11.10（5%）、15.79（1%），澳门元-10.73（10%）、12.53（5%）、21.15（1%））。因此，接受原假设H0∶？籽1=？籽2=0，即人民币兑港币和人民币兑澳门元的实际汇率存在单位根，呈非平稳状态，长期购买力平价在中国大陆和中国香港以及中国大陆和中国澳门之间并不成立。

4.针对线性结果的单位根检验

针对Wald检验中呈现线性特征的人民币兑台币实际汇率序列，我们分别使用ADF和PP检验法对其进行单位根检验，具体结果如表4所示。

由如上检验结果所示，PP检验中，显著性水平在1%、5%和10%的临界值均大于统计量；而在ADF检验中，统计量均大于相同显著性水平下的临界值。PP和ADF检验均给出了被检测时间序列为单位根过程的结论。综合以上，可认为人民币兑台币实际汇率为非平稳线性过程，长期购买力平价在中国大陆和中国台湾之间同样不成立。

五、研究结论

本文基于1995年1月至2011年12月的月度汇率和CPI数据，运用门限自回归（TAR）模型，根据Wald检验的线性与非线性检验结果，采用相应的非线性单位根和线性单位根检验方法分别探讨了人民币兑港币、澳门元和台币实际汇率的非线性和平稳性问题，核心结论主要包括以下几点：

（1）从实证研究的直观结果来看，人民币兑港币和澳门元的实际汇率呈非线性关系，进一步，我们通过非线性单位根检验发现其过程不具备平稳性。而人民币兑台币的实际汇率的非线性关系被研究所拒绝，结合PP和ADF单位根检验我们发现其过程存在单位根，同样不具备平稳性。以上实证研究成果表明，从长期来看，中国大陆与中国香港、中國澳门之间的购买力平价关系并不成立。

（2）从对未来相关研究的启示上来说，本文得出的人民币实际汇率存在非线性的实证结果在诸多采用非线性方法对此问题进行的研究基础上进一步表明：一方面，针对该问题的研究如仍采用单一的、传统的线性研究方法或将取得不稳健的结论，对购买力平价的检验应更多地采用非线性的方法来进行；另一方面，考虑到ADF、PP等方法在检验非线性平稳性方面的局限性，未来基于实际汇率长期趋势的研究应更多地采用针对检验非线性关系的平稳性而设计的、且更为有效的方法。

（3）從本文研究结论的经济意义来看，两岸三地作为一个整体，其货币与人民币之间的真实汇率在长期内呈现非稳定性的状态，代表其货币与人民币之间的名义汇率与相应的相对价格水平之间在长期并不存在稳定的内在关联性，我们认为此结果代表了时间趋势并非相关名义宏观经济变量之间的长期内在关系中的一个关键因素，并由此可以在一定程度上合理推知：若时间趋势构成了某一特定宏观经济变量长期走势的一个显著因素，则此变量必定包含有某种真实的成分。进一步来说，由于时间序列变量的相对跨度并非检验购买力平价的决定性因素，故长期购买力平价本质上或许可以被视为是内生经济收敛的一种结果。

[参考文献]

[1] Caner， M. and Hansen， B..Threshold autoregression with a unit root[J]. Econometrica， 2001， 69：1555 -1596.

[2] Yu-shao，Liu.， Meng-nan Zhu and Chi-wei Su. Threshold effects in purchasing power parity for African countries[J].African Journal of Business Management， 2011，5（24）：10235-10242.

[3] Narayan， Paresh Kumar and Narayan， Seema. Are real exchange rates nonlinear with a unit root？ Evidence on PPP for Italy： a note[J]. Applied Economics， 39： 19， 2483-2488， First published on： 27 August 2007.

[4] Narayan， Paresh Kumar. Are the Australian and New Zealand stock prices nonlinear with a unit root？ [J]. Applied Economics， 2005， 37： 18， 2161-2166.

[5] J.Madsen， V.Mishra， R.Smyth. Are labour force participation rates non-stationary？ Evidence from 130 years for G7 Countries[J]. Australian Economic Papers 2008， Blackwell Publishing Ltd.

[6] Sofiane H. Sekioua. Nonlinear adjustment in the forward premium： evidence from a threshold unit root test[J]. International Review of Economics and Finance， 2006，15：164 -183.

[7] Heckscher， E. F.. Vaxelkurens grundval vid pappersmynfot[J]. Economisk Tidskrift， 1916，18： 309 -12.

[8] Chortareas， G. and Kapetanios， G..The yen real exchange rate may be stationary after all： evidence from nonlinear unit root tests[J]. Department of Economics Working Paper Series， University of London， 2003.

[9] Chortareas， G. E.， Kapetanios， G. and Shin， Y.. Nonlinear mean reversion in real exchange rates[J].Economics Letters， 2002，77， 411 -7.

[10] Kapetanios， G.， Shin， Y. and Snell， A..Testing for a unit root in the nonlinear STAR framework[J].Journal of Econometrics， 2003， 112：359 -79.

[11] Kilian， L. and Talyor， M. P..Why is it difficult to beat the random walk forecast of exchange rates[J]. Journal of International Economics， 2003， 60：85 -107.

[12] Narayan， P. K.. New Evidence on purchasing power parity from 17 OECD countries[J]. Applied Economics， 2005，37：1063 -71.

[13] Narayan， P. K.. Are bilateral real exchange rates stationary？ Evidence from Lagrange multiplier unit root tests for India[J]. Applied Economics， 2006，38： 63 -70.

[14] Taylor， M. P.， Peel， D. A. and Sarno， L..Nonlinear mean reversion in real exchange rates： towards a solution to the purchasing power parity puzzles[J]. International Economic Review， 2001，42：1015 -42.

[15] 王志強，齐佩金，孙刚.人民币汇率购买力平价的界限检验[J]. 数量经济技术经济研究， 2004（2） .

[16] 刘金全，云航，郑挺国.人民币汇率购买力平价假说的计量检验：基于Markov区制转换的Engel- Granger协整分析[J].管理世界， 2006（3）.

[17] 丁剑平，谌卫学.汇率非线性因素在部分亚洲货币汇率中的特征：检验购买力平价论的新方法[J]. 财经研究，2010（2）.

[18] 张卫平.购买力平价非线性检验方法的进展回顾及其对人民币实际汇率的应用[J].经济学，2007（4）.

[19] 封福育.中、日、韩长期购买力平价的非线性检验-基于1994-2009 年的经验分析[J]. 统计与信息论坛，2011（1）.

Research on the Purchasing Power Parity of China Mainland，

China HK SAR， China Macau and China Taiwan： Empirical Examination Based on Canner & Hansen TAR Model

Zheng Xiaoya， You Haibo

（School of Economics， Xiamen University，Xiemen 361005，China）

Abstract： Purchasing power parity theory is a basic and wide applied exchange rate determination theory. Previous studies focus on whether the purchasing power parity is justified by using a traditional testing method of inspecting time sequence stability under the framework of linear or non-linear model. This paper， under the framework of non-linear study， based on Based on Canner & Hansen TAR Model， tries to make an empirical study of the Purchasing Power Parity of China Mainland， China HK SAR， China Macau and China Taiwan， targeting at threshold effect and stationarity of non-linearity state by using Wald test and non-linearity unit root test. The main findings is that the real exchange rate of RMB against KH D and MOP is non-linear， while the real exchange rate of the RMB against NT$ does not have threshold effect. In addition， the result of stationarity test based on non-linearity shows that the long-term Purchasing Power Parity may not hold among China Mainland， China HK SAR， China Macau and China Taiwan.