动态过程数学模型参数估计的最小二乘方法LeastWord文件下载.docx
- 文档编号:7531928
- 上传时间:2023-05-08
- 格式:DOCX
- 页数:24
- 大小:624.96KB
动态过程数学模型参数估计的最小二乘方法LeastWord文件下载.docx
《动态过程数学模型参数估计的最小二乘方法LeastWord文件下载.docx》由会员分享,可在线阅读,更多相关《动态过程数学模型参数估计的最小二乘方法LeastWord文件下载.docx(24页珍藏版)》请在冰点文库上搜索。
J=yTy+TT-yT-TTy
=yTy+TT-2TTy式(2-1-4)
假设:
(T)(n+1)(n+1)满秩,由
利用线性代数的以下两个矩阵对向量求偏导数的公式:
和
有:
所以:
解出参数估计向量:
Ls=(T)-1Ty式(2-1-5)
P=(T)-1则参数估计向量Ls=PTy
参数估计向量Ls被视为以下“正则方程”的解:
(T)=Ty式(2-1-6)
注:
为了便于区别,我们用红体字符表示估计量或计算值,而用黑体表示为参数真值或实际测量值。
三、关于参数最小二乘估计Ls性质的讨论
以上求解参数最小二乘估计Ls时并为对{k}的统计特性做任何规定,这是最小二乘估计的优点。
当{k}为平稳零均值白噪声时,则Ls有如下良好的估计性质:
a)参数最小二乘估计Ls是y的线性估计
Ls=PTy是y的线性表出;
b)参数最小二乘估计Ls是无偏估计,即ELs=(参数真值)
[证明]:
ELs=E[PTy]=PTE(y)=PTE(+)=
PT+E()=+0=
c)最小二乘估计Ls的估计误差协方差阵是2P(n+1)(n+1)
即:
E[(Ls-)(Ls-)T]=2P
E[(Ls-)(Ls-)T]=E[PT(y-
)(y-)TP]=E[PTTP]=PTE(T)P=
PT2INNP=2P
d)若{k}为正态分布零均值白噪声时,则Ls是线性无偏最小方差估计(证明从略)。
如若{k}是有色噪声,则Ls不具有上述性质,即为有偏估计。
四、最小二乘估计Ls的的几何意义和计算问题
1。
最小二乘估计的几何意义
最小二乘估计的模型输出值为yk=kTLsk=1,2,…N
输出实际测量值与模型输出值之差叫残差:
k=yk–yk
模型输出向量为y=Ls,而残差向量为:
=y–y=y–Ls
Tk=Ty–T(T)-1Ty=Ty–Ty=0
即残差向量与由测量数据矩阵的各个向量:
1,2,…,N张成的超平面(估计空间)正交,而最小二乘模型输出向量y为实际输出向量y在估计空间上的正交投影,这就是最小二乘估计的几何意义。
2.关于最小二乘估计计算中的病态问题
估计参数向量Ls一般是求解正则方程:
得出。
可以利用消元法等一系列求解多元线性一次方程组的方法,计算得出,其有解的条件是(T)=P–1矩阵非奇异(行列式数值大于零)。
但有时在求解式(2-1-6)方程组是会出现矩阵接近于奇异(行列式数值接近于零),即所谓“病态”的情况。
由此导致参数估计的结果不稳定,不可信。
出现上述情况的原因可能是由于被辨识的过程受到的外加激励不够,采样间隔太密;
或者A/D转换的位数太短,计算舍入误差累计所致。
为解决最小二乘计算中可能出现的病态问题,提出了不少改进算法,例如:
Householder变换法、改进的平方根法和U—D分解算法。
后者是Bierman1977提出的改善P阵计算性质(对称性、正定性和稳定性)而又不增加计算量的算法。
正定P阵可以分解成一个上三角阵U(其对角线元素都为1)和一个对角阵D
P=UDUT
由此可解决最小二乘计算中可能出现的病态问题,具体可参阅关于计算方法的文献。
总之,我们在使用最小二乘的辨识方法时,应该注意避免出现和克服病态问题。
应用举例
在建立生产过程的静态模型时,特别是在机理不清之时常用多元线性回归方法,例如:
水泥凝固放热量与水泥成分的关系模型
y=a0+a1x1+a2x2+a3x3+a4x4+
y水泥凝固时的放热量(卡/克);
x1~x1水泥的几种成分。
五、非线性最小二乘法(NonlinearLeastSquare)
以“误差平方总和为最小”的估计准则,估计非线性模型参数的方法。
假设非线性静态系统模型为
y=f(x,)
非线性模型f的形式是已知的,参数未知。
经过N次实验,取得N组数据(x1,y1)(xN,yN)。
准则:
需要用优化算法求解,常用的算法有两类——搜索法和迭代法。
前者如单纯型法;
后者如梯度法、高斯法、牛顿—拉夫森法、变尺度法等等。
该类方法也还可应用于动态模型和时间连续模型的参数估计。
单纯型法是先给出参数空间的几个猜想点,构成正多面体,计算各点的目标函数值,比较各值后舍去最差的点,按照反射、开拓、收缩等步骤确定新的估计点,直到预定的精度要求后停止搜索。
迭代法是先猜想一个估计的初值,确定一个向量为可接受的方向和步长,进行迭代计算k+1–k=μ•v…。
具体内容可阅有关计算方法的参考书籍。
近年来发展出一系列基于生物进化论的优化新算法,如遗传算法(GeneticAlgorithms)(基于改进遗传算法的系统辨识方法.北方工业大学学报.第10卷,第1期,1998年3月,P62-67.)和免疫算法(ImmuneAlgorithms),使得优化算法的性能得到了很大的改善。
应用举例:
化学反应速度模型
2—2动态过程参数估计的线性最小二乘法
一、模型类:
考虑CAR模型
式(2-2-1)
其中{y(k)}和{u(k)}为可测的输出和输入,{(k)}为不可测的随机干扰。
上式还可表示成:
A(z-1)y(k)=B(z-1)u(k)+(k)式(2-2-2)
A(z-1)=1+a1z-1+a2z-2+。
+anz-n
B(z-1)=b1z-1+b2z-2++bnz-n
还可表示为
式(2-2-3)
当进行了k=1-n,2-n,..,0,1,2,…,N共计(N+n)次采样,得到N个方程:
y
(1)=-a1y(0)-…-any(1-n)+b1u(0)+…+bnu(1-n)+
(1)
y
(2)=-a1y
(1)-…-any(2-n)+b1u
(1)+…+bnu(2-n)+
(2)
………………………………………………
y(N)=-a1y(N-1)-…-any(N-n)+b1u(N-1)+…+bnu(N-n)+(N)
用矩阵表示成
yN=N+N式(2-2-4)
yN=[y
(1),y
(2),…,y(N)]T
N=[
(1),
(2),…,(N)]T
二、参数最小二乘估计Ls的导出
估计准则为
式(2-2-5)
由
解出:
N=(NTN)-1NTyN式(2-2-6)
上式可视为以下正则方程的解
(NTN)N=NTyN式(2-2-7)
称为最小二乘的“一次完成算法”,是离线算法,有解的条件是(NTN)2n2n满秩。
用消元法或平方根法解线性方程组,得出N。
三、用配方法导出Ls
=(yN-N)T(yN-N)+yNTN(NTN)-1NTyN
-yNTN(NTN)-1NTyN=
{-(NTN)-1NTyN}TNTN{-(NTN)-1NTyN}
+yNTyN-yNTN(NTN)-1NTyN
上式的后两项中均不含,能使得J=nim的条件是:
=(NTN)-1NTyN即前式(2-2-6)
用两种方法推证出相同结论。
四、线性动态参数最小二乘估计Ls的性质
在静态模型:
yk=kT+k式(2-1-2)
中的k=[1,x1,,xN]T为确定性量,取值与yk统计性质无关;
在动态模型:
y(k)=kT+(k)式(2-2-3)
的最小二乘估计Ls虽然形式上与静态的相同,但是式的k中包含y(k-1)、y(k-2)、…,导致有关估计的统计性质的证明要困难得多,不能简单地套用静态模型多元回归的结果。
动态参数最小二乘估计Ls的估计性质的主要结果是:
当N时:
EN=0(渐进无偏估计)
a.s.
N(强一致性收敛)
如若{(k)}为有色噪声,Ls是有偏估计。
五、数字仿真
y(k)-1.5y(k-1)+0.7y(k-2)=u(k-1)+0.5u(k-2)+(k)
一般取N=100~200。
仿真结果如下表:
参数
a1
a2
b1
b2
真值
-1.5
0.7
1
0.5
估计
-1.5091
0.7074
1.0022
0.4656
2-3递推最小二乘方法(RLS法)
一、递推算法的导出
对式(3-2-1)的CAR过程{y(k)}和{u(k)}进行了k=1-n、…、N,共(n+N)次观测,组成了yN和N,可得出参数估计N,现在再进行一次新的采样,又得出N+1和新的估计N+1。
如何由估计向量N经过递推直接得到由新的估计向量N+1,而不必反复做一次完成LS法的计算?
先分析yN+1与yN以及由估计向量N+1与N的关系,用分块矩阵:
式(2-3-1)
式(2-3-2)
由式(3-2-6)
有由估计向量
N+1=(N+1TN+1)-1N+1TyN+1=PN+1N+1TyN+1
其中
=(NTN+N+1N+1T)-1=(PN-1+N+1N+1T)-1
根据以下矩阵求逆公式:
若A是nn满秩矩阵,B和C是nm阵,且(A+BCT)nn和(I+CTA-1B)mm都满秩(I为m维单位阵),则有以下矩阵恒等式成立:
(A+BCT)–1=A–1-A–1B(I+CTA-1B)–1CTA-1
现令:
P–1=A和N+1=B=Cm=1,考虑到前面得到的PN+1=(PN+N+1N+1T)-1得出
式(2-3-3)
其中:
即相当于矩阵恒等式中的(I+CTA-1B)–1
将式(2-3-3)、式(2-3-2)和式(2-3-1)代入由估计向量
N+1=PN+1N+1TyN+1并令:
式(2-3-4)
将上式打开:
N+1=PNNTyN+PNN+1y(N+1)–
-KN+1N+1TPNNTyN–KN+1N+1TPNN+1y(N+1)
N+1=N+(PNN+1-KN+1N+1TPNN+1)y(N+1)–KN+1N+1TN
可以证明:
(PNN+1-KN+1N+1TPNN+1)=KN+1式(2-3-4)
N+1=N+KN+1(y(N+1)–N+1TN)式(2-3-5)
二、递推算式的物理含义
N+1=N+KN+1(y(N+1)–N+1TN)式(2-3-5)
式(2-3-4)
式(2-3-3)
N+1T=[-y(N),-y(N-n+1),u(N),,u(N-n+1)]式(2-3-6)
模型:
y(N+1)=N+1T+(N+1)
N时刻对N+1时刻的预报
y(N+1N)=N+1TN式(2-3-7)
预报误差(被称为新息),用绿色表示
y(N+1N)=y(N+1)-y(N+1N)=y(N+1)-N+1TN
式(2-3-8)
则式(2-3-5)可表达成
N+1=N+KN+1y(N+1N)式(2-3-5)
物理意义:
新的参数估计N+1是对上次老的估计N进行修正而得出的,修正的依据是利用在N对新的输出y(N+1)预报的预报误差。
KN+1是修正系数向量,它需要递推计算得出,在递推计算KN+1时要用到估计误差的协方差阵PN,而后者也是递推得出的。
三、初值选择和计算框图
递推计算需要初值0和P0。
初值的选择有两种方法:
其一是,可以用初始的3n组数据用一次完成算法解出:
P0=(T)-1和0=P0Ty,但是实际上并不用上述办法,而是用以下更为简便的方法。
即令:
0=0和P0=I2n2n,其中=103~106(即假设为足够大的正数)。
在递推2n步后,估计结果与前面介绍的精确初值相接近。
递推计算框图如下图所示。
四、LS法和RLS法的讨论
(1)LS和RLS法数学等价
a)均由使准则J=[y(k)-kT]2=min得出
b)不要求对{(k)}的统计特性有任何验前知识。
c)如果{(k)}为零均值白噪声,则可得渐进无偏估计,即当n时,E=0,且。
d)若{(k)}为有色噪声,是有偏估计,但是因该算法简单,所以应用广泛。
e)均可推广到多输入多输出系统。
开始
置初值
0N
由式(2-2-6)
构成N+1
进行第N+1次采样
由式(2-2-4)NN+1
计算KN+1
由式(2-2-5)计算N+1
由式(2-2-3)计算PN+1
否
终点判断
是
输出模型
(2)LS法和RLS法的比较
a)LS法是一次完成算法,适于离线辩识,要记忆全部测量数据,程序长;
b)RLS法是递推算法,适于在线辩识和时变过程,只需要记忆n+1步数据,程序简单;
c)RLS法用粗糙初值时,如若N较小时,估计精度不如LS法。
2–4时变过程的参数估计(时变递推最小二乘法)
参数的递推最小二乘估计与一次完成的最小二乘估计是数学等价的,它们都仅适用于估计时间定常过程的参数,而不适用于估计时变过程的参数。
时变过程的特点是过程的参数可能随着时间变化而改变。
因此,它的数学模型参数具有“时间性”,在利用动态过程的输入—输出数据来辨识模型参数时,“老”的数据往往只能反映“老”的过程参数;
而改变后的“新”参数,要靠用新的和比较新的实验数据来估计。
因此时变过程参数估计的特点是,不同时段的实验数据的作用是有区别的。
时变过程的参数估计有多种不同的算法,限于课时,本课程只讲授最常用的“带遗忘因子的递推最小二乘估计算法”,以下简称为“遗忘因子法”。
一、遗忘因子算法的思路
算法的主要思路是“厚今薄古”,即对新老数据给予不同的对待,逐渐遗忘老数据的影响。
具体做法是每采得一个新的y(N+1)时,将以前的所有数据乘以小于1的加权因子(0<
<
1),
乘以加权因子后的估计准则为:
即对不同时刻的残差予以不同加权的平方和,老数据的作用按照指数衰减,被遗忘。
二、遗忘因子递推最小二乘算法的递推算式
利用矩阵求逆公式得出:
=2,遗忘因子(0<
1);
取值范围(0.95~0.995),值愈小,“遗忘”愈快。
带遗忘因子的RLS法由以下几个递推算式组成:
式(2-4-1)
式(2-4-2)
三、遗忘因子的作用
对于遗忘因子的不同值,得到了不同的遗忘效果:
值较小时的估计跟踪参数时变的能力强,但是噪声干扰影响造成的估计波动亦大,可通过以下的计算机仿真说明值的影响。
模拟如下2阶过程:
y(k)+a1y(k-1)+a2y(k-2)=b1u(k-1)+b2u(k-2)+(k)
a2=-0.7,b1=1,b2=0.5,而
a2=0.6(N<
1000时)
a2=1.2(N≥2000时)
给出值分别为1;
0.995;
0.99;
0.95时的a1和b1:
四、算法的改进
1)变算法
模型的一步预报误差为
(k)=y(k)-kTk-1
将2(k)与事先给定的m1值相比较(m1>
0)
当2(k)>
m1时选择较小的值
当2(k)≤m1时选择较大的值
一般选择m1=(4~6)Var{(k)}。
2)变P算法
将2(k)与事先给定的m2值相比较(m2>
m2时将P阵重新置为:
106I
比较:
以2000步的预报误差平方总和Q做比较
Q的理想值为2000
RLS法
变算法
=0.999/0.9
m1=4
变P算法m2=25
=0.995
=1
=0.99
2889
6270
2502
2122
2137
3)附加R阵法
只要在P阵的递推公式上附加一个小小的正定附加方阵R,即可大大提高的跟踪能力。
为事先给定的正值
例1、y(k+1)=a1y(k)+b1u(k)+(k+1)
a1=-0.20k40;
a1=-0.8k40
b1=0.40k40;
b1=1.6k40
例2.、参数按照正弦规律快速周期性变化的一阶过程
a(k)=0.5-0.4Sin(k20)
b(k)=1+0.8Sin(k20)
2–5最小二乘法的应用
例1.热交换器在线辨识(R.Isermann“Automatica“)
辩识以蒸汽流量调节阀为输入,热水出口温度为输出的动态模型,t=0.3s、N=31、n=3,用RLS法。
结论:
对采样的数据进行适当的滤波处理,对辩识的效果大有改善。
例2.直流它激电动机参数在线辨识“自动化学报”
直流电机模型的形式是已知的,主要参数如:
机电时间常数Tm、电磁时间常数Ta、磁通量和负载转距Mc是需要辩识的,其中:
假设Mc是常值,其方向总是与电机旋转方向相反,即:
方案一、利用u和的采样值辩识Tm、Ta、Ce和Mc
方案二、利用i和的采样值辩识J和Mc
辨识结果如下表所示
N
8
10
20
30
40
50
J(km.m2)
0.00365
0.00465
0.00448
0.00450
0.00451
0.00454
MC(N.m)
0.4380
0.4644
0.4630
0.4549
0.4552
0.4149
实验电机ZK-22FC,110V,5.7A,2800rpm,0.5KW,产品目录给出的转动惯量J=0.00425(km.m2)。
可见,当N=20以后,估计结果与实际很接近了。
2–6关于卡尔曼滤波(KalmanFiltering)与递推最小二乘参数估计
卡尔曼(RudolfEmilKalman)美籍数学家兼电气工程师,1930年生于布达佩斯,哥伦比亚大学博士,斯坦福大学教授。
主要学术贡献为系统的能控性与能观测性和卡尔曼滤波。
关于卡尔曼滤波的主要结果:
基于状态空间描述并受噪声干扰的信号进行滤波,提取有用信号的方法。
在此之前是著名的维纳滤波,但是维纳滤波是基于频率域的,用积分方程描述,不好计算,而且只适用于平稳过程,计算时需要利用全部数据,不便于在线计算。
卡尔曼滤波是基于时间域的,用状态方程描述,计算仅需用有限步的观测数据,特别适合于在线计算,可以推广到非平稳过程。
线性离散随机系统,用离散时间状态方程和观测方程描述为:
k=1,2,…
x(k)—状态向量y(k)—观测输出
w(k)和v(k)均为零均值正态白噪声干扰。
初值x(0)和E(x)=x0,以及方差Varx0=P0均已知,最小方差
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 动态 过程 数学模型 参数估计 最小 方法 Least