《系统辨识与建模》
课时：40
参考书：
1 徐南荣、宋文忠、夏安邦 系统辨识，东南大学出版社，1991
2 方崇智、萧德云 过程辨识，清华大学出版社，1988
3 Ljung L,And Soderstrom T,Theory and Practics of Recursive Identification,MIT.Press,Cambridge
4 [美]P.艾克霍夫 系统辨识-参数和状态估计，潘科炎等译，科学出版社，1977
5 [美]夏天长 系统辨识—最小二乘，熊光楞、李芳芸译，清华大学出版社，1983
期刊：
Automatica
Proc,IFAC Identification and System Parameter Estimation
IEEE Trans,On Automatic Control
自动化学报控制理论与应用
学习要求
1 培养独立学习一门新课程的能力，为今后学习和研究打下基础（要求大家尽量少依赖听课，多自学）。
2 掌握基本的辨识理论和辨识技术
3 能独立设计辨识实验，并编程计算
4 学习一些现代建模技术
考核办法
给出一个数据文件，通过编程对其进行辨识，并写出报告
第一讲 概论
实体与模型
实体：客观存在的事物及其运动状态，有时也称之为“系统”
模型：实体的一种简化描述。模型保持实体的一部分特征，而将其它特征忽略或者变化。不同的简化方法得到不同的模型。
模型分类：
直觉模型：地图、建筑模型、照片、软件演示文档等
物理模型：风洞、水力学模型、传热学模型、电力系统动态模拟模型等。（缩小的复制品）
数学模型：描述实体中一些关系和特征的数据模型。例如：投入/产出模型、热源与室温的关系模型等。
数学模型
数学模型还可分为：
图表模型：如阶跃响应、脉冲响应、频率响应、温度与热电偶输出关系表解析模型：代数方程、微分方程、差分方程、状态方程程序模型：神经网络仿真程序语言模型：模糊关系模型
获得数学模型的方法有：
经验总结法：模糊关系模型、静态线性关系模型
机理分析法：解析模型
实验法,图表模型
数据拟合法：解析模型、程序模型
用数据拟合法获得解析模型的过程即为系统辨识。
系统辨识
L,A,Zadeh[1962],辨识就是在输入和输出数据的基础上，从一组给定的模型类中，确定一个与所观测系统等价的模型。
这一定义给出了系统辨识的三要素：数据、模型类和准则。
数据：由观测实体而得。不唯一，受观测时间、观测目的、观测手段等影响。
模型类：规定了模型的形式。不唯一，受辨识目的、辨识方法等影响。
准则：规定了模型与实体等价的评判标准。不唯一，受辨识目的、辨识方法等影响。
系统辨识的三要素是评判数据拟合方法优劣的必要条件，只有在相同的三要素下，才可区分数据拟合方法的优劣；而在不同的三要素下，这种结论也会改变。（图1）
辨识目的明确模型应用的最终目的是很重要的，因为它将决定如何观测数据、如何选择三要素以及采用什么数据拟合方法等。而最根本的是它将影响辨识结果。
辨识目的主要取决于模型的应用。辨识模型应用有以下几个方面：
辨识目的 选择 准则
模型类
辨识实体 数据观测 数据 数据拟合 模型
图 1 系统辨识三要素
1 验证理论模型；要求：零极点、结构（阶次及时延）、参数都准确；模型类同理论模型。
2 设计常规控制器；要求：动态响应特性、零极点、时延准确；便于分析的模型类。
3 设计数字控制器；要求：动态响应特性、时延准确；便于计算机运算的模型类。
4 设计仿真/训练系统；要求：动态响应特性准确；便于模拟实现的模型类。
5 预报预测；要求：动态响应特性、时延准确；便于计算机运算的模型类。
6 监视过程参数，实现故障诊断；要求：参数准确；能直观体现被监视过程参数的模型类。
系统的定量与定性分析；要求：静态关系准确；模型简单，便于人脑判断。
辨识的一般步骤
我们将结合一个实际例子来说明辨识的一般步骤。
例：
上图为长网造纸的流程简图。对于造纸企业来说，质量控制就是要控制好成品纸的定量与水份。而纸的定量与水份与纸浆浓度D、纸浆流量F、车速V及蒸汽压力P都有关系：
G=f(D,F,V,P)
W=g(D,F,V,P)
为了采用计算机对上述过程进行控制，需要建立数学模型。
这就是系统辨识的第一步：明确辨识目的。
为了实现这个目标，我们要进一步了解系统的一些特点。
这是第二步：收集先验知识。
经过现场调查，我们发现：
1 车速调整存在同步困难，而不同步会引起断纸，因此，通常将车速设为恒定；
2 流量的改变到定量的改变存在约60秒的延迟，而响应过程只有约2秒；（保持浓度不变）
3 浓度的改变到定量的改变存在约120秒的延迟，而响应过程约80秒；（保持流量不变）
4 蒸汽压力的改变到水份的改变存在约45秒的延迟，而响应过程约60秒；（保持浓度与流量不变）
第三步：设计辨识试验。
辨识试验的目的是使采集到的数据能反映系统的动态特性，因此要对系统进行分块，设计对分块后系统施加的激励信号，设计数据采集时的采样频率。对于本例，其中的一个分块为流量和蒸汽压力对定量、水份的影响；
试验时，保持车速和纸浆浓度不变；
对流量和蒸汽压力，分别施加伪随机序列扰动，幅度以不引起断纸为限；
设定采样频率为2，试验时间为1000秒。
采集信号为：定量、水份、纸浆流量和蒸汽压力第四步：现场准备。（接线图）
现场准备要做以下几件事：
向企业领导申请试验时段；
准备扰动信号发生器，并通过预发信号，检验扰动信号是否准确；
测验现场信号的干扰情况，必要时设计模拟信号滤波器；
准备模数转换设备，调好信号的零迁和放大参数；
现场接线，将生产设备、试验设备与计算机连接。
第五步：数据采集。
将采集到的数据存盘，并编写数据说明文件；
第六步：数据预处理。
对采集到的原始数据进行变送器非线性校正、数字滤波、标准化、重抽样等加工，使数据适合辨识工具的处理，同时也应满足模型要求。
以上步骤为数据观测过程。
第七步：选择模型类。
选择模型类的工作有两部分：其一是选择应用模型，通常应依据辨识目的来选择；其二是选择参考模型，参考模型是便于进行结构辨识和参数估计的模型第八步：结构辨识与参数估计。
应用辨识理论和方法编制程序，对第六步所得的数据进行拟合，得到参考模型的阶次和参数。
第九步：模型检验。
对所得到的参考模型按评判准则进行检验，如不达要求，则分析问题所在，并返回到前期各相应步骤。
第十步：模型转换。
将参考模型转换为应用模型。
第十一步：应用评价。
从应用角度评价模型，如不符合应用要求，应分析问题所在，并返回到相应步骤。
下图描述了辨识各步骤之间的关系。
思考题：你认为系统辨识还有用吗？
明确辨识目的 收集先验知识
设计辨识试验
现场准备
数据采集
数据预处理
选择模型类
结构辨识
参数估计
不通过
模型检验 分析问题
通过
模型转换
不成功
应用评价 分析问题
成功
最终模型
系统辨识一般流程第二讲 辨识三要素一、数据本节介绍辨识数据的特点及获得适宜辨识的数据的方法。
随机过程X(t):在每一个时间点（t0）上，都是一个随机变量，其概率密度函数p(x,t)随时间变化。
平稳随机过程：在所有时间点上，概率分布都相同的随机过程，其概率密度函数p(x)不随时间变化。
各态遍历平稳随机过程：从整个时间轴上看，每个随机事件都会发生的平稳随机过程。其谱密度函数与概率密度函数类似。时间平均等于集合平均。
数字特征特征
随机过程
平稳随机过程
各态遍历平稳随机过程
均值
(期望值)
均方值
方差
=
自相关函数
协方差函数
=
=
=
一些关系
=
互相关函数
=
=
互协方差函数
=
=
其它
说明：离散计算时假设采样时间间隔为To，则时延τ=*To。
相关函数的性质
1 ≥0
2
3
4 若x(t)是周期为T的信号，则其自相关函数也是周期为T的信号。即：
5 若x(t)=y(t)+z(t)，且y(t)与z(t)互不相关（），则
6若x(t)=y(t)+z，其中，z是一个常数，则
7若x(t)均值为零，且不含有周期性成分，则当τ很大时，x(t)与x(t+τ)必然是互相独立的（不相关），因此，,τ充分大。
8 若x(t)均值为零，则。这是因为在通常情况下，等于向下平移，因此，当=0时，两者相等。
9 对于线性系统y(k)=G(z)u(k)，有
10 虽然x(t)是个随机过程，但却不是随机过程，而是一个确定性的时间函数。
Parseval定理与功率谱
Parseval定理：确定性信号x(t)的总能量为：
确定性信号x(t)的平均功率：
确定性信号x(t)的平均谱密度：
随机性信号x(t)的平均谱密度：
维纳—肯塔金关系式：
随机过程x(t)的谱密度与自相关函数构成一组傅立叶变换对：
定义互谱密度为互相关函数的傅立叶变换：
应用维纳—肯塔金关系式，可以证明，对于频率响应为的线性系统，在随机输入下的输出谱密度和互谱密度分别为：
输出谱密度关系告诉我们：要充分激励系统，就要使输入信号的频谱“宽”于系统频谱。
白噪声
如果一个零均值、平稳随机过程的谱密度为常数，我们称之为白噪声（由白色光联想而得）。白噪声有以下特点：
1
2 ，频谱宽度无限。
∞，τ=0
0，τ≠0
3 ，其中，为Dirac函数，即=
且
4 无记忆性，即t时刻的数值与t时刻以前的过去值无关，也不影响t时刻以后的将来值。从另一意义上说，即不同时刻的随机信号互不相关。
白噪声的用途：
1 作为系统输入时，有，τ=0,1,2,…，即为系统的单位脉冲响应。
2 作为被辨识系统输入时，可以激发系统的所有模态，可对系统充分激励；
3 作为被辨识系统输入时，可防止数据病态，保证辨识精度。
4 在辨识过程中，以输出估计误差是否具有白色性来判断辨识方法的优劣，也可用来判断模型的结构和参数是否合适。
5 产生有色噪声。
白噪声的产生方法：
1 （0，1）均匀分布白噪声：{，i=1,2,3…}，初值可取为：
2 正态分布白噪声：，
其中ξ为服从（0，1）均匀分布的白噪声。
有色噪声有色噪声是指每一时刻的噪声和另一时刻的噪声相关，因而其谱密度也不再是常数。在工业生产实际中，白噪声在物理上是不存在的，常见的往往是有色噪声。
有色噪声的表示定理：设平稳噪声序列{e(k)}的谱密度是ω的实函数，则必定存在一个渐近稳定的线性环节，使得在输入为白噪声序列的情况下，环节的输出是谱密度为的平稳噪声序列{e(k)}。
白噪声 线性环节 有色噪声
(成形滤波器)
{w(k)} H(z-1 ) {e(k)}
M序列（二位式最大长度伪随机序列）
例：4阶M序列 1 1 1 1 0 1 0 1 1 0 0 1 0 0 0 1 1 1 1 。。。。。
1 1 1 1 –1 1 –1 1 1 –1 –1 1 –1 –1 –1 1 1 1 1 。。。。。
特点：
周期性，周期长度为N=2n-1 (n为阶次)，是n个移位寄存器所能表示的最多状态数
M序列中某种状态连续出现的段称为游程。一个周期中有2n-1个游程，游程长度为1 ——n，但出现的概率是随机的。长度为1的游程有2n-2个，长度为2的游程有2n-3个，长度为3的游程有2n-4个，以此类推，最后，长度为n的游程有1个。
一个周期中1的个数比0的个数多1（例中有8个1，7个0）
若以-1代替0，则序列的自相关函数RM(τ)为
RM(τ)=1 τ=i*N,i=0,1,2,…
RM(τ)=1/N τ=i*N+k,i=0,1,2,…,k=1,2,…N-1
当N充分大时，M序列M1：{M(k)}与它的J步移位序列M2：{M(k+J)}在τ=0，1，…J-1和τ=J+1，…N时是不相关的，即
M序列的频谱：设为一步M序列信号的持续时间，N是一个周期的持续时间。M序列的频谱为：
其中：
上式告诉我们：
1 M序列的频谱不是光滑的曲线，而是线条谱。
2 M序列的直流分量（）与N2成正比，因此，加大N，可减少M序列中的直流分量。
3 M序列的频带为，因此，减少，可增加带宽。
4 谱线密度与N成正比，
应用,
将M序列作为扰动信号有以下好处：
1 幅值可取a,-a，容易选择，且当N充分大时，均值约等于0
2 在一个周期内，自相关函数RM(τ)近似为δ函数，因此，以M序列为输入的线性系统，其互相关函数序列等于脉冲响应序列（N大于过渡过程）
3 对于多输入单输出系统，可将同一M序列的不同移位序列（例：{M(k)}、{M(k+J)}、{M(k+2J)}等）作为各输入信号的扰动信号，当输出对各输入的过渡过程小于J时，可认为输入信号之间是互相正交的。
逆重复M序列：
将M序列与方波序列相乘，得逆重复M序列，其特点为：
周期=2N，一个周期中1的个数与-1的个数相等，自相关函数RM(τ)为
RM(τ)=1 τ=i*N*2,i=0,1,2,…
RM(τ)=-1 τ=i*N*2+1,i=0,1,2,…
RM(τ)=0 τ=i*N+k,i=0,1,2,…,k=1,2,…N-1
用作辨识系统的扰动信号的优点在于：具有白噪声的优点，而其幅值可控。
辨识试验设计设计原则 在安全的前提下，尽可能地激励系统；保持输入输出关系；适当解耦
明确目的与要求 模型用途了解辨识对象
划清要辨识系统的边界，选好输入/输出，从边界外连入的其它信号尽量保持稳定，并作为被辨识系统的噪声。
整体/局部 ——确定哪些输入需要叠加扰动信号，哪些输入要保持稳定输入/输出/噪声——确定是多输入还是单输入（耦合关系）、确定过渡过程是否有明显差异（时间常数）、了解噪信比的大小（滤波）及噪声类型（白色、有色）
值域范围——确定信号采集时是否需要零迁、放大安全工况——可叠加扰动信号的类型与幅值选择工况 生产负荷、试验时间、系统区域隔离、地理区域隔离、安全措施扰动信号设计要点,扰动信号频带应宽于过程的工作频带；持续时间为3-5个扰动信号周期；幅值由安全工况确定，对可中断生产的系统，且试验不引起原材料浪费的，可实施单独试验。
扰动信号类型——M序列/白噪声序列扰动信号幅值——由安全工况确定
M序列设计——级数n、步宽Ts、正交化
设系统过渡过程时间为T，最高工作频率为fmax，（通常，fmax=(3~15)/T）则以1.5*T< (2n-1)*Ts和1/Ts>3* fmax为原则，选择n 和Ts。
实际应用过程：
1 确认系统过渡过程T;
2 选择最高工作频率fmax=(3~15)/T；（设为fmax=k/T）
3 根据1/Ts>3* fmax，得出Ts=T/(3k)；
4 根据1.5*T< (2n-1)*Ts=(2n-1)*T/(3k)或4.5k=2n-1，求出n；
另一种设计方法：
1 确认系统过渡过程T，选择n；
2 根据n Ts=(0.1~0.8)T，算出Ts=(0.1~0.8)T/n。
对于多输入单输出系统，要求(2n-1-1)*Ts>1.5*∑(输出对各输入的过渡过程时间）。
白噪声序列——1/Ts>3* fmax
数据采集设计采样间隔——满足香农定理，设采样间隔为T0，1/T0>2* fmax。对于有扰动信号的系统，T0< Ts，且Ts应是T0的整数倍。
实验时间——数据长度：3-5个扰动信号周期，预实验1-2个扰动信号周期物理信号处理——信号零点选择、滤波器选择、放大器选择
（现场准备扰动信号发生器——测试、加限幅器、选择叠加方式、接线信号采集装置——测试、时间间隔标定、接线）
试验设计报告 上述内容及人员组织，报领导批准
数据预处理传感器非线性校正——根据传感器生产厂提供的静态非线性曲线，编制程序对相应的采样数据进行校正。
滤波——对于输入输出加相同的滤波器，输入输出关系不变重抽样——对于过渡过程较长的子系统，可用重抽样的方法，重新选择采样间隔，以降低模型阶次标准化——零均值、单位方差目的：降低舍入误差，防止病态方程，提高计算精度方法,y(k)+ay(k-1)=bu(k-1)+cu(k-2)+d-------------------(1)
∑y(k)/N+ ∑ay(k-1)/N= ∑bu(k-1)/N+ ∑cu(k-2)/N+d
=∑y(k)/N= ∑y(k-1)/N
=∑u(k-1)/N= ∑u(k-2)/N,d=(1+a) -(b+c)
Y(k)= y(k)- ,U(k)=u(k)-
Y(k)+aY(k-1)=bU(k-1)+cU(k-2)----------------------------------(2)
σ2y=∑Y(k)2/N,σ2u=∑U(k)2/N
Yσ(k)=Y(k)/ σy,Uσ(k)=U(k)/ σu
Yσ(k)+a Yσ(k-1)=b Uσ(k-1)* σu /σy +c Uσ(k-2) * σu /σy
B= b * σu /σy,C= c* σu /σy
Yσ(k)+a Yσ(k-1)=B Uσ(k-1)+C Uσ(k-2)-----------------------(3)
b= B*σy/ σu,c= C*σy/ σu ----反演公式二、辨识准则误差定义参数误差,J=∑/n,，适合于算法研究
方程误差,J=∑[e(k)]2/N e(k)= y(k)+y(k-1)-u(k-1)- u(k-2)，J是参数的线性函数；只能作一步预测；适合模型为：AY=B1U1+B2U2
输出误差,J=∑[e(k)]2/N,e(k)= y(k)-(k),(k)= -(k-1)+u(k-1)+u(k-2)，J是参数的非线性函数；当输入序列已知时，可作多步预测；适合模型为：
u(k) 过程 y(k) u(k) 过程 y(k)
B/A B/A
-1 /
e(k)
白化误差加权准则,J=∑[e(k) ]TW[e(k) ] /N
模型结构对误差的影响模型阶次越高，J值越小，导致模型过估计采样间隔对模型结构及误差的影响采样间隔增大一倍，模型阶次大致减少一半；
要进行不同采样步长下模型的比较，误差公式应改为：
J=∑[e(k) ]TW[e(k) ]*T02 /N,T0为采样步长三、数学模型对象特性静态模型 输入输出之间的关系不随时间变化动态模型 未来的输入输出关系受过去及现在的状态影响，模型的有效性不随时间变化动态建模 模型的有效性随时间变化，需要经常改变
集中参数/分布参数 集中参数模型只反映变量与时间的关系，用常微分方程描述；分布参数模型反映变量与时空之间的关系，用偏微分方程描述；
序列（权函数）/方程 序列（权函数）反映系统输出在标准输入下的响应；方程反映系统输入输出的关系，给定输入可计算出系统输出；
h(τ),y(t)=∫h(τ)u(t-τ)dτ
y(n)(t)+a1y(n-1)(t)+…+any(t)=b1u(n-1)(t-τ)+…+bmu(t-τ)
=AX+BU Y=CX+DU
线性/非线性 线性关系是：若y1=f(u1),y2=f(u2),则 ay1+by2=f(au1+bu2)
{hn(τ1,…,τn)},y(t)= ∑∫…∫hn(τ1,…,τn)u1(t-τ1)…un(t-τn)dτ1…dτn
=G(X,U)=GLX+GR(X,U) Y=F(X,U)
时变/时不变 模型的结构不边，但参数随时间变化。例如锅炉中的热交换器，由于灰尘的附着而改变热交换的效率
y(n)(t)+a1(t)y(n-1)(t)+…+an(t)y(t)=b1(t)u(n-1)(t-τ)+…+bm(t)u(t-τ)
确定型/随机型 若系统含有随机扰动，则称为随机型的。实际中很少没有随机扰动的，只是强度不同而已。扰动强度的度量用噪信比N/Sy=σW/ σX,(y=x+w)
单变量/多变量 Y(s)=G1(s)U1(s)+ G2(s)U2(s)
离散型/连续型 y(k)= ∑h(i)u(k-i)
A(q-1)y(k)=B(q-1)u(k)+w(k)
Y(z-1)=G(z-1)U(z-1)
X(k+1)=Ax(k)+Bu(k) y(k)=Cx(k)+Du(k)
有时延/无时延Y(s)=e-sτG1(s)U1(s)+ e-sτG2(s)U2(s)
y(k)+a1y(k-1)+…+any(k-n)=b1u(k-τ)+…+bmu(k-m-τ)
常用模型阶跃响应（函数、序列、曲线）
差分方程/微分方程连续传递函数/离散传递函数连续状态方程/离散状态方程一般描述/典型环节惯性环节 k/(1+Ts)n
微分环节 kTs/(1+Ts)n
积分环节 1/Tas(1+T0s)n
惯性+微分环节 k(1+aTs)/(1+Ts)n
模型间的相互转换
差分方程 离散传递函数 离散状态方程
阶跃响应
微分方程 连续传递函数 连续状态方程
差分方程与离散传递函数,A(q-1)y(k)= B(q-1)u(k) Y(z-1)=U(z-1)
微分方程与连续传递函数的关系与此类同。
差分方程、离散传递函数到阶跃响应：以阶跃信号为输入，通过计算机仿真得到阶跃响应。
阶跃响应到连续传递函数：
阶跃响应曲线转换为典型环节传递函数
1设系统传递函数为W(s)，在阶跃输入下，系统响应的拉氏变换为F(s)=W(s)/s，若y(kT) k=0,1,2,…为系统的阶跃响应序列，则F(s)的离散计算公式为F(s)=T∑e-skT y(kT)；
2 根据{y(kT)}的曲线形状，可选择一个典型环节传递函数w(s)，并令G(s)=w(s)/s为其阶跃响应的拉氏变换；
3 Jt=∫|g(t)-y(t)|dt,等价于 Js=|G(s)-F(s)|，或者 J=∑[G(sj)- F(sj)]2，sj充分小
4 对于充分小的s,选择sj>0，j=1,2,…,L（L大于G(sj)中参数个数）,可计算出F(sj)，同时G(sj)（以及J）中只含w(s)的未知参数；
5 令 =0，求出参数，即得G(s)，然后w(s)=s G(s)
离散传递函数与连续传递函数：
1 离散信号的保持
u(t) y(t)
G(s)
u(z)=Z[u(t)]=Z[u(kT)] y(z)=Z[y(t)]=Z[y(kT)]
u(kT) x(t) y*(t)
L(s) G(s)
u(z)=Z[u(kT)] X(s) y*(z)=Z[y*(t)]
G(z)= y*(z)/ u(z)=Z[G(s)L(s)] G(s)=Z-1[G(z)]/L(s)
若要求y*(z)= y(z)，则应满足y*(kT)= y(kT)，进而要求x(t)≈u(t)，这通过适当选择L(s)来实现。
L1：x(t)=u*(t)=，x(s)= =u*(s)，L(s)=1(或T)
当u(t)=δ(t)时，y*(t)= y(t)，故称为脉冲响应不变等效。G(z-1)=Z[G(s)]，G(s)=Z-1[G(z-1)]
L2：（零阶保持器）x(t)=u(kT)，kT≤t＜(k+1) T，L(s)=(1-e-Ts)/s
当u(t)=1(t)时，y*(t)= y(t)，故称为阶跃响应不变等效。G(z-1)=(1-z-1)Z[G(s)/s]，G(s)= Z-1[G(z-1)]/[(1-e-Ts)/s]
L3：x(t)=u(kT)*((k+1)T-t)/T+u((k+1)T)*(t-kT)/T，kT≤t＜(k+1)T，L(s)= (1/T+s)[(1-e-Ts)/s]2
斜坡响应不变等效。
Z变换与逆 z=eTs,s=ln(z)/T：（可用查表法进行）
G(z)=Z[G(s)]=G(s)| s=ln(z)/T，G(s)=Z-1[G(z)]=G(z)| z=eTs
近似变换：向后差分 s=
向前差分 s=
双线性变换 s=，
离散传递函数与离散状态方程
Y(z-1)=U(z-1)= H(z-1)U(z-1) X(k+1)=Ax(k)+Gu(k) y(k)=Cx(k)+Du(k)
从离散状态方程到离散传递函数,H(z-1)=C(zI-A)-1G+D
从离散传递函数到离散状态方程,设A(z-1)=1+a1z-1+......+anz-n,
B(z-1)=b0+b1z-1+......+bnz-n
则,0 1 0........,0 0 g1=b1-a1b0
0 0 1........,0 0 g2=b2- a1 g1-a2b0
A=,.....................................,G=,................................................
0 0 0,.........,0 1 gn-1=bn-1- a1 gn-2-..........- an-2 g1-an-1b0
-an –an-1 –an-2.........-a2 –a1 gn=bn- a1 gn-1-......................- an-1 g1-anb0
C=[1 0,...............0],D=b0
连续状态方程与连续传递函数的关系与此类同。
连续状态方程与离散状态方程:
=AX+BU Y=CX+DU X(k+1)=Fx(k)+Gu(k) y(k)=Hx(k)+Eu(k)
F=eAT≈I+AT，G=A-1(F-I)B，H=C,E=D (阶跃响应不变)
A=lnF/T≈(F-I)/T,B=(F-I)-1AG,C=H,D=E (阶跃响应不变)
辨识中的参考模型
ARMAX模型 A(q-1)y(k)= ∑Bi(q-1)/Fi(q-1)ui(k)+ D (q-1)/C(q-1)e(k)
AR自回归 X外生变量 MA滑动平均线性回归模型 y(k)=-a1y(k-1)-…-any(k-n)+b1u(k-1)+…+bmu(k-m)+ ν(k)
θ=[a1,…,an,b1,…,bm]T
φ=[-y(k-1),…,-y(k-n),u(k-1),…,u(k-m)]T
y(k)= θTφ+ν
伪线性回归模型 当φ中含有不可观测的变量时，通常用其估计值来代替。因此，含有待估计变量的线性回归模型称为伪线性回归模型。
预报误差模型 y=f(Y,U,θ) +ν，v为零均值的新息，f(Y,U,θ)为y的最优估计。
例：
考虑模型,y(k)=-ay(k-1)+bu(k-1)+e(k)+ce(k-1)
若e(k)、e(k-1)...,e(0)都是已知或可估计的，则e(k)的估计为：
(k)= y(k)+ay(k-1)-bu(k-1)- c(k-1)
则k时刻的输出估计为：
(k)= -ay(k-1)+bu(k-1)+ c(k-1) *
=-ay(k-1)+bu(k-1)+ c[y(k-1)+ay(k-2)-bu(k-2)- c(k-2)]
且有：c(k-1)= -acy(k-2)+bcu(k-2)+ cc(k-2)
两式相加得：
(k)+ c(k-1)= (c-a)y(k-1)+bu(k-1) **
此表明：若(0)是最优预报，则(1)、...、(k)都是最优预报。由此也得到预报误差模型为,f()：(c-a)y(k-1)+bu(k-1)- c(k-1)
v：y(k)- (k)= e(k)+ce(k-1)-c[y(k-1)- (k-1)]
==e(k)
或者,y(k)= (c-a)y(k-1)+bu(k-1)- c(k-1)+ e(k)
对于一般的差分模型：A()y(k)=B()u(k)+C()w(k)
其预报误差模型为：y(k)= (k)+ w(k)
其中,w(k)=y(k)-u(k)
(k)=(1-) y(k) +u(k)
请同学们将离散传递函数模型改写成预报误差模型，并用伪线性回归模型的形式表现之。
第三讲 参数估计的批量法最小二乘算法考虑差分方程：A()y(k)=B()u(k)+ w(k)，其中w(k)为白噪声。假定模型的结构已知（n,m），将其写成线性回归模型：
y(k)=-a1y(k-1)-…-any(k-n)+b1u(k-1)+…+bmu(k-m)+ w(k)
θ=[a1,…,an,b1,…,bm]T
φ(k)=[-y(k-1),…,-y(k-n),u(k-1),…,u(k-m)]T
y(k)=φT(k)θ+w(k)
若我们的观测数据可写出N个这样的等式，
YN=Φθ+WN,ΦT=[φ(1），。。。，φ(N)]
ΦT YN=ΦTΦθ+ΦT WN
θ=(ΦTΦ)-1ΦT YN-(ΦTΦ)-1ΦT WN
θLS=(ΦTΦ)-1ΦTYN,条件1：E(ΦT WN)=0 （w(k)白色）
条件2,ΦTΦ可逆以上结果等同于求使：J=Σ(y(k)-φT(k)θ)2最小的θ，因此称为最小二乘算法。
最小二乘算法的MATLAB程序：
clear all; %最小二乘法for MISO
load t2; %读入数据矩阵
tt=t2;
clear t2;
ll=size(tt); %得到矩阵维数200x2
r=ll(2)-1;
m(1)=input('order of A(z)');
tao(1)=0;
for i=1:r
i
m(i+1)=input('order of B(z)');
tao(i+1)=input('timedelay of B(z)');
end
if m(1)<max(m),m(1)=max(m);end
n=m(1)+max(tao);
lll=ll(1)-n;
in1=r+1;
kn=0;
for k=1:in1
for i=1:lll
for j=1:m(k)
jtao=j+tao(k);
ff(i,j+kn)=tt(i+n-jtao,k);
if k==1,ff(i,j+(k-1)*n)=-ff(i,j+(k-1)*n);end
end;end;
kn=kn+m(k);
end;
for i=1:lll
yy(i)=tt(i+n,1);
end;
fa=ff'*ff;
fay=ff'*yy';
zz=inv(fa);
zta=zz*fay
m，tao
加权最小二乘：若对各次观测数据加不同的权，即求使Jα=Σαk(y(k)-φT(k)θ)2最小的θ,则得到参数的加权最小二乘估计：
θLSα=(ΦTΑΦ)-1ΦTΑYN
参数估计的一般性质无偏性
如果E()=0，（=-） 或E()=，则称估计为无偏的。
有效性如果无偏估计满足Cov()=M-1，则称估计为有效的。其中：
称为Fisher信息矩阵，其逆M-1称为Crammer-Rao 下界。在一般情况下，有Crammer-Rao不等式：
Cov()=E{(-)(-)T}≥M-1
对于z=Hθ+w，若噪声w为零均值、协方差阵Σw=E(WWT)的正态分布噪声，即：
W～N(0，Σw)，则输出信号Z～N(E(Hθ)，Σw)，于是，
，因此：
，于是，M=E(HTΣw-1H)，其Crammer-Rao不等式为：Cov()≥E(HTΣw-1H)-1。
有效估计也称为最小方差估计、马尔可夫估计。
一致性若估计为渐近无偏的（E()=0），且Var()=0，则称为一致估计。Var()=E{(-)T(-)}
最小二乘、加权最小二乘估计性质估计性质
最小二乘
加权最小二乘
无偏性
= (ΦTΦ)-1ΦT WN，当:
1 w(k)为零均值独立随机过程（白噪声）时；或者
2 E(ΦT WN)=0时，有：E()=0。无偏条件较严格。
= (ΦTΑΦ)-1ΦTΑWN
无偏条件同最小二乘估计。
有效性
Cov()=E{(ΦTΦ)-1ΦTWNWTNΦ(ΦTΦ)-1}
=(ΦTΦ)-1ΦTΣwΦ(ΦTΦ)-1≥M-1,若W是服从正态分布的独立随机过程，则Σw=I，因而Cov()=I(ΦTΦ)-1= M-1。
Cov()=(ΦTΑΦ)-1ΦTΑΣwΑΦ(ΦTΑΦ)-1，当Α=Σw-1时，Cov()=(ΦTΣw-1Φ)-1= M-1，此时不要求W是服从正态分布的独立随机过程
一致性
若W是服从正态分布的独立随机过程，则Cov()==0,其中收敛于一个有界非奇异矩阵。而Var()只是Cov()的对角线上元的和，因而也有：Var()=0。
同最小二乘估计
噪声方差估计
e= YN-=Φθ+WN-Φ=Φθ+WN-Φ{(ΦTΦ)-1ΦTYN }=WN-Φ(ΦTΦ)-1ΦT WN
=(I-Φ(ΦTΦ)-1ΦT) WN=B WN,因为BT=B,B2=B
所以，E(eTe)=E(WN TB WN)= E(WN T WN- WN TΦ(ΦTΦ)-1ΦT WN)
=N- (TrΦ(ΦTΦ)-1ΦT)= (N-dimθ)
=Σ(y(k)-)2/(N-dimθ)
在有色噪声环境下，最小二乘估计是有偏的。下面的一些算法对最小二乘进行改进。
广义最小二乘考虑差分方程：A()y(k)= B()u(k)+ w(k)，其中w(k)为白噪声。假定模型的结构已知（n,m）。
如果噪声模型C()已知，显然用C()对输入/输出数据进行滤波，则可得到满足最小二乘估计无偏条件的模型：A()(k)= B()(k)+ w(k)，其中：
(k) = C()y(k)，(k) = C()u(k)。在C()未知时，我们可考虑采用迭代估计的方法去求得。
广义最小二乘的计算步骤为：
1 令C0()=1，i=0下标表示迭代次数；ξ0=1000000
2 计算(k) = Ci()y(k)，(k) = Ci()u(k)；i=i+1;
用最小二乘估计A()(k)= B()(k)+ w(k)中的参数；
用估计模型()、()以及各时刻的观测数据，计算出残差(k):
(k)= ()y(k)- ()u(k)
计算ξi=Σ2(k)及ξ=ξi-1-ξi，如果ξ小于一定数，则结束辨识。否则转下一步。
对于噪声模型C()(k)= w(k)，用最小二乘估计出参数，得到更新的Ci()后，返回2。
以上算法的每一次循环中都要进行滤波和两次求逆。下面的算法将在计算工作量上有所改进。
偏倚校正算法
仍考虑差分方程：A()y(k)= B()u(k)+ w(k)，其中w(k)为白噪声。
令(k)= w(k)，则A()y(k)= B()u(k)+ (k)。分别写成回归模型：
Y=Φθ+，=C+W，组合起来有Y=[Φ，]+W，其最小二乘解为：
=，利用分块矩阵求逆公式，有：
=(ΦTΦ)-1ΦTY-(ΦTΦ)-1ΦT
=D-1TMY
M=I-Φ(ΦTΦ)-1ΦT D=TM
须要注意，的求取仍然是一个迭代过程。偏倚校正算法的计算步骤为：
1 令C()=1，用最小二乘法求=θLS=(ΦTΦ)-1ΦTY，并保留Φ、Γ=(ΦTΦ)-1ΦT以及M=I-ΦΓ。
2 计算=Y-Φ，并依据=C+W构造，计算D=TM。
3 计算=D-1TMY，并计算=Γ及=θLS-。
4 若参数已收敛，则结束辨识，否则转2。
以上算法的一次循环中没有滤波，且只有一次求逆。如果将第3步中的计算改为：
=(T)-1T，则还可省去D的计算。（这一改进由夏天长首先给出。）
辅助变量法分析最小二乘法中，在Y=Φθ+W的各项乘上ΦT,然后利用ΦT W的期望值为零，得到参数的无偏估计。受此启发，若在Y=Φθ+的各项乘上ΨT，使其满足以下两个条件：1.ΨT的期望值为零；2,ΨTΦ可逆，则也可得到参数的无偏估计。
下面讨论辅助变量的选取：
设模型为A()y(k)= B()u(k)+ (k)，若u(k)与(k)不相关：
a 选取辅助模型D()z(k)= F()u(k)，用z(k)、u(k)构造Ψ；
b 若系统的纯时延τ已知，则可用u(k-τ)、u(k)构造Ψ；
c 用u(k)、u(k)构造Ψ
d (k)= D()w(k),若D()的阶次n已知，则可用y(k-n)、u(k)构造Ψ；
e 先求出最小二乘解θLS=(ΦTΦ)-1ΦTY,然后依据()z(k)= ()u(k)计算出输出估计z(k)，再用z(k)、u(k)构造Ψ；
多步最小二乘考虑差分方程：A()y(k)= B()u(k)+ w(k)，其中w(k)为白噪声。模型可改写为C()A()y(k)= C()B()u(k)+ w(k)或：
D()y(k)=F()u(k)+w(k)，其中:
D()=C()A()，F()= C()B()---------**。
此模型可用最小二乘解出D()、F()。这是第一步。
第二步可有两种不同方法：
a 解同次幂方程组由**式，可得关系：A()F()=B()D()，两边分别展开，并按的同次幂相等规则，可列出na+nb+nc个方程：
F=Hθ，其中θ=[a1,...ana,b1,...bnb]T，F=[f1,...fnb+nc,0,...0]
0 0 0...........,0 1 0 0...........,0 0
-f1 0 0..........,0 -d1 1 0.............0 0
H= -f2 –f1 0...........,0 -d2 -d1 1.............0 0
,...................................................
,...............................................-d1 1
,.............,-f2 -f1,.......................-d2 –d1
,.............,-f3 -f2...............................
,....................................................
-fnb+nc................................................
0 -fnb+nc...............,-dna+nc.......................
0 -fnb+nc..........,0 -dna+nc.......................
....................................................
0 0..........,-fnb+nc 0........................-dna+nc
na列 nb列用最小二乘可求得θ的无偏估计，即A()、B()。此法也可用来求C()。
B 传递函数等价降阶
由**式，可有,F()/D()=B()/A()，此说明两个传递函数是等价的。对F()/D()施加激励信号{u(k)}，可得输出z(k)= F()/D()u(k)。用最小二乘法处理{ z(k)，u(k)}，选择合适的阶次，可得A()、B()的无偏估计。
相关最小二乘设模型为A()y(k)= B()u(k)+ (k)，若u(k)与(k)不相关，则用u(k-j)乘模型中的各项并求期望，得：A()Ruy(j)= B()Ru(j)，用最小二乘法可得A()、B()的无偏估计。
第四讲 辨识原理
随机逼近法设广义误差e(k)是参数估计值θ的函数，参数辨识问题可通过极小化e(k)的方差来实现。即求参数θ使下列准则函数最小：J(θ)=1/2E{e2(k)}。J(θ)的负梯度为：
=E{e(k)}。如果可求解=0，则可求得参数的估计。但当e(k)的分布未知时，实际上是不可求解的。
在计算数学中，求二次函数的极小值常采用迭代法。首先给出参数的一个估计值，以二次函数在该参数估计值处的负梯度为修正方向，选取适当的步长后，修正参数估计值，直到收敛。仿此，我们有：θ(k+1)= θ(k)+ρ(k)。如果在求时不求期望，则得到一个随机的迭代算法，称之为随机逼近法：
例：考虑线性回归模型：y(k)=φT(k)θ+e(k)，其中e(k)是零均值随机噪声。J(θ)=1/2E{[ y(k)-φT(k)θ]2},=E{φ(k)[ y(k)-φT(k)θ]}。如果假定e(k)是各态遍历的，则梯度为零可用E{φ(k)[ y(k)-φT(k)θ]}=φ(k)[ y(k)-φT(k)θ]=0来代替，由此得到了最小二乘法。
应用随机逼近法，可得：
θ(k+1)= θ(k)+ρ(k)φ(k)[y(k)-φT(k)θ(k)]。
ρ(k)的选取应保证迭代收敛，可选取满足如下条件的ρ(k),
ρ(k)>0且；；，例如ρ(k)=1/k；
若以准则函数的二阶导数（即海赛矩阵）之逆来参与选择修正方向，则称为牛顿法：
θ(k+1)= θ(k)+ρ(k)R(k)-1φ(k)[y(k)-φT(k)θ(k)]，其中：
R(k)==E[φ(k)φT(k)]=R(k-1)+ρ(k)[φ(k)φT(k)-R(k-1)]
牛顿法的优点是收敛速度快。
模型参考自适应辨识方法（简略）
考虑线性回归模型：y(k)=φT(k)θ+e(k)，其中e(k)是零均值随机噪声。以输出估计误差为反馈信号，以PI调节器的方式来修正参数：
θ=θI+θP,θI(k)= θI(k-1)+Pφ(k)ε(k)，θP(k)=Qφ(k)ε(k)
其中P为对称正定矩阵，Q满足P/2+Q>0
ε(k)为广义误差，最简单的取法为ε(k)= y(k)-φT(k)θ(k-1)
极大似然法输出z是一个随机变量，它的概率密度p(z|θ)取决于参数θ。当获得观测序列ZL={z(1),z(2),...,z(L)}T时，由该观测序列组成的联合概率密度p(ZL|θ)应当取得最大值。（当一个随机事件发生了，我们有理由相信，外部条件一定处于使这随机事件发生的概率最大时的状态。）那么，θ的极大似然估计就是使p(ZL|θ)|θML=max的参数估计值。
由于p(ZL|θ)中ZL已知，因而它只是参数θ的函数，故称它为θ的似然函数。有时也记作L(ZL|θ)。极大似然原理可用下列等价的表示方式：
，，，
求解极大似然估计的下一步是要给出p(ZL|θ)的具体描述。
独立观测设z(1),z(2),...,z(L)是一组在独立观测条件下获得的随机序列，即各观测值是互相独立的，则p(ZL|θ)可简化为：p(ZL|θ)= p(z(1)|θ)p(z(2)|θ)...p(z(L)|θ)
= p(z(k)|θ)，其对数似然函数为：L(ZL|θ)=ln p(z(k)|θ)。
设z(k)～N(m,σ)，即：p(z(k)|θ)=，则其负对数似然函数为：
- L(ZL|θ)=+Llnσ+（L/2）ln2π
当σ已知时，准则函数就是：J(θ)=
当σ未知时，可先由min(J(θ))求及Jmin，再由min(- L(ZL|θ))求。
，
二、非独立观测若z(1),z(2),...,z(L)是非独立观测条件下获得的随机序列，即观测值z(k)是在已有观测z(1),z(2),...,z(k)的基础上得到的，则p(ZL|θ)应按条件概率的乘法规则写成：
p(ZL|θ)= p(z(L)|ZL-1,θ)p(ZL-1 |θ)，以此类推，有：
p(ZL|θ)= p(z(k)|Zk-1,θ)。
设z(k)～N(mk,σk)，并取mk为z(k)的条件均值：(k)=E[z(k)| Zk-1]，即：p(z(k)| Zk-1,θ)=，则其负对数似然函数为：- L(ZL|θ)= +（L/2）ln2π。
当σk≡σ已知时，准则函数就是：J(θ)=
当σk≡σ未知时，可先由min(J(θ))求及Jmin，再由min(- L(ZL|θ))求。
例：考虑模型：z=+，其中w～N(0,σ2I),则～N(,σ2I)，
p(|θ)=，
J=lnL(θ)=
，，，，
，i=1,2,...,n
，i=1,2,...,n
，i=1,2,...,n
，i=1,2,...,n
，→
直接求解以上非线性方程组是困难的，如果已有参数的某个估计，并用其构成如下预滤波器：，，，则上面前两个方程组可写成：
，，这等价于辅助变量法。对于噪声模型，我们已有：
，定义：，则w=C/Dv=Cv-(D-1)w，再次构造预滤波器：v*=v/D，w*=w/D,
则，，这也等价于辅助变量法。
预报误差估计法对于模型结构已知的系统，如果我们获得参数的某个估计，则可对输出进行预报，那么用预报误差的大小来衡量参数估计的优劣也是合理的。定义输出估计误差：
e(k)=y(k)-(k)，将收敛于e(k)的协方差阵。通常采用的预报误差有：
J1(θ)=Tr(ΛD(θ)),J2(θ)=log det(D(θ))。可以证明，当e(k)～N(0,Σe)时，若J1(θ)中的Λ取为Λ=LΣe-1，则预报误差估计等价于极大似然估计；
例：考虑模型：z=+。利用已经获得的参数估计和噪声w的到k-1时刻为止的估计，我们有输出估计：，并且可得噪声估计：=，由此可见，以为输出估计误差时，预报误差法与极大似然法的准则函数是等价的。
若记e(k)= (k)，或写成Ce=Dw，那么对Ce的估计可写成Ce=(D-1)w，或者，
于是，我们有输出估计：，那么输出估计误差（也作为噪声估计）：
Bayse估计
参数是与观测值关联的随机变量。在给定的观测数据下，使得验后的条件概率密度p(θ|ZL)最大的参数最有可能是真实参数。根据全概率公式：
p(θ|ZL)= p(z(L)|ZL-1,θ)p(θ|ZL-1)/p(z(L)|ZL-1)
= p(θ|Z0)
其中p(θ|ZL-1)为参数的验前概率密度函数，p(z(L)|ZL-1,θ)为数据的条件概率密度函数，p(z(L)|ZL-1)是与参数无关的项。根据验后条件概率密度p(θ|ZL)，有两种方法可求得参数估计，一是极大验后参数估计法，另一是条件期望参数估计法，它们统称为Bayes方法。
1 极大验后参数估计法参数的极大验后估计必定满足：=0=+
若令上式右边第一项为零，所得的参数估计即为极大似然估计。可见，极大验后参数估计比极大似然估计多考虑了参数的验前概率知识，因而一般情况下极大验后参数估计优于极大似然估计。但是参数的验前概率密度不易获得，因此，极大似然估计比极大验后参数估计用得更普遍。
2 条件期望参数估计法直接以参数的条件期望作为参数估计值是合乎逻辑的。(L)=E{θ|ZL}。可以证明，条件期望参数估计等价于极小化参数估计误差的方差,minE{[θ-(L)]T[θ-(L)]|ZL}
因此也称为最小方差估计。
例,考虑线性回归模型：y(k)=φT(k)θ+e(k)，其中e(k)是零均值随机噪声，
e(k)～N(0,)。设参数θ在数据Z0条件下的验前概率密度为：
p(θ|Z0)=，其中N=dimθ。那么，参数θ在数据ZL-1条件下的验前概率密度为：
p(θ|ZL-1)=。由噪声e(k)～N(0,)可知，
p(z(L)|ZL-1,θ)=，于是有：
log p(θ|ZL)=const--
=const-
若令，则log p(θ|ZL)= const-，显然极大验后参数估计应该是，而从正态分布的形式来看，它也是参数的条件期望估计。
第五讲 参数估计的广义递推算法广义递推算法的一般结构假定对象模型可用伪线性回归方程描述：y(k)=φT(k)θ+w(k)。设是参数θ的一个估计，e(k,)是广义误差，它是参数估计误差的一种描述。定义二次准则函数：
V(t,)=，其中Λ是对误差矢量各分量的加权矩阵，α(k,t)是对k=1,...,t 这个时间串上所有e(k,)的加权序列，满足：
α(k,t)=λ(t)α(k,t-1),k<t
α(t,t)=1，λ(t)≤1
当λ(t)≡1时，α(k,t)≡1。当λ(t)<1时，称为遗忘因子。记在时刻t使V(t,)最小化的参数估计为：t，则应有：。又若已知在时刻t-1处的最优估计t-1，则将V(t,)在t-1处展开成台劳级数：
V(t,)=
V(t,t-1)+[-t-1]+[-t-1]T[-t-1]+o(|-t-1|)
由可得：
+[t-t-1]= o(|t-t-1|)，若d*d矩阵可逆，则：
t=t-1-{}-1+ o(|t-t-1|)
定义记号：=-，则：
=-=--
=λ(t)-
=λ(t)+ -
假定：
1）=0，t-1为时刻t-1处的最优估计
2）=0
3）≈=R(t-1,t-2)
t=t-1+ R-1(t,t-1)，R(t,t-1)= λ(t)R(t-1,t-2)+ *
若定义：P(t)= R-1(t,t-1)，则P-1(t)=λ(t)P-1(t-1) +，
于是P(t)=，t=t-1+P(t)
此为广义递推算法的基本结构。如果再考虑辅助变量算法，则广义递推算法的一般结构为：P(t)=，t=t-1+P(t)
广义误差的选择假定对象模型可用伪线性回归方程描述：y(t)=φT(t)θ+w(t)。为了构造广义误差，我们先来构造输出估计。先定义记号：
t=[，，，]T，，，，)的构造中为广义输出,为输入观测值，为白色噪声估计，为有色噪声估计。对于广义输出，可选的例子有：
=[-y(t-1),-y(t-2),...,-y(t-n)]T,即系统输出的观测值。
=[-(t-1|t-1),-(t-2|t-1),...,-(t-n|t-1)]T,即系统输出的验前估计值。
=[-(t-1|t),-(t-2|t),...,-(t-n|t)]T,即系统输出的验后估计值。
有色噪声估计可选为0或：
=[-(t-1),...,-(t-r)]T,其中(t)=y(t)-，，0，0) t-1,或
=[-(t-1),...,- (t-r)]T,其中(t)=y(t)-，，0，0) t；
白色噪声估计可选为0或验前输出估计误差或验后输出估计误差：
=[(t-1),...,(t-r)]T,其中(k)=y(k)-，，，)t,
=[(t-1),...,(t-r)]T,其中(k)=y(k)-，，，)t-1
广义误差e(k,)通常为验前输出估计误差，即e(k,)=(k)，尽管如此，由于、及的不同选取，可构造的算法多达27种。
递推算法收敛性分析主要有三种分析方法：
ODE 法：将递推算法转化为伴随微分方程，通过研究微分方程的稳定性来判断递推算法的收敛性。
鞅理论：引进随机Lyapunov函数，并用鞅理论分析这个函数的收敛性，进而得到参数估计的收敛性。
间接法：通过研究R(t)t的收敛性，间接得到参数估计的收敛性。
ODE 法的启发式讨论令t=R(t,t-1)，并假定λ(t)=1，则：
t=t-1+ ，=+(-) **
t-t-1=dθ，=dt，-=dR，f(θ)= ,g(θ)= ，则有：
f(θ)，g(θ)-
ODE法的应用例请参考Ljung的著作。
Martingale收敛定理先介绍一些概念：
σ-代数：设S由集合X的某些子集组成，（X本身为S的一个元素），若S对“差”及可列和封闭，则称S为σ-代数。
递增的σ-代数序列：设S={At,t∈T}为σ-代数序列，若At(，(t∈T；
At(As，t≤s，则称S为递增的σ-代数序列。
随机变量序列向σ-代数的转化：对所有的t，若随机变量序列Xt可以从σ-代数序列At观测得到，即E{ Xt | At }= Xt，则称随机变量序列Xt可转化为σ-代数序列At。换言之，包含在At中的信息足以完全确定Xt。
Martingale：若①{ At }为递增的σ-代数序列；②{ Xt }可转化为σ-代数序列At；
③E(|Xt |)<∞；④Xs=E(Xt| As,s<t)，则称序列{ Xt }为对于{ At }的Martingale。
定理：设{ Xt }为对于{ At }的向量Martingale，且对于所有的t，均有||E{ Xt Xt T}||<∞，则存在一具有有界协方差的向量随机变量X，使得Xt→X,(a.s.)，Xt=E(X| At),(t>0。
Martingale收敛定理：设{Tn}、{αn}、{βn}均为非负随机变量序列，它们可可转化为σ-代数的递增序列Fn，使得E[Tn|Fn-1]≤Tn-1-αn-1+βn-1，若 a.s.(几乎处处)，则
Tn→T (a.s.)，且 a.s,
对于算法：
t=t-1+ R-1(t,t-1)，R(t,t-1)=R(t-1,t-2)+ *
定义随机变量：η(t)=e(t)[1-ψTR-1(t)ψ]，若η(t)可表示为：
η(t)=H-1(z-1)[ψT ]+w(t)，=t-θ0
其中w(t)满足,E[w(t)|Ft-1]=0，E[w2(t)|Ft-1]=σ2，e(t)-w(t)∈Ft-1，
Ft-1是由w(0),...,w(t-1),ψ(0),...,ψ(t-1)构成的σ-代数。
则当H-1(z-1)-1/2严格正实时→0 a.s.
以上结论的证明需引用Martingale收敛定理，其中：
Tt=[T R(t)+2αtβt]/t，αt=η(t)- w(t)+ψT/2，βt=-ψT
典型算法
P(t)=，t=t-1+P(t)
针对Ay=Bu+Cw的方程误差法
——RELS递推增广最小二乘计算式：φ(-y,u,)=[-y(k-1),...,-y(k-n),u(k-1),...,u(k-m),(k-1),...,(k-r)]T;
(k)=y(k)-(k),(k)=φT(-y,u,),
=[a1,...,an,b1,...,bm,c1,...cr]T;
误差结构图,
e(k)=Cw(k)/A
u(k) 过程 y(k)
B/A
(z-1) (z-1)
( z-1)-1 (k)
分析：
(k)= y-u-(-1)(k)=(-)u+[e-(-1)(k)]
若e=w，且u与w不相关，则→w意味着：→，→
但不能保证→A，→B，→C
针对Ay=Bu+Cw的输出误差法计算式：φ(-,u,)=[-(k-1),...,-(k-n),u(k-1),...,u(k-m),(k-1),...,(k-r)]T;
(k)=y(k)-(k),(k)=φT(-,u,),
=[a1,...,an,b1,...,bm,c1,...cr]T;
误差结构图,
e(k)=Cw(k)/A
u(k) 过程 y(k)
B/A
(z-1) -1(z-1) (k)
( z-1)-1
分析：
(k)= u+e-[u+(-1)(k)]=( -)u+[e- (-1)(k)]
若e=w，且u与w不相关，则→w意味着：→,→(-1)+1
若e=w，且u与w不相关，则→w意味着：→，→
针对Ay=Bu+w/D的方程误差法计算式：φ(-y,u,-ξ)=[-y(k-1),...,-y(k-n),u(k-1),...,u(k-m),-ξ(k-1),...,-ξ(k-r)]T;
ξ(k)=y(k)-φT(-y,u,0),(k)=y(k)-(k),(k)=φT(-y,u,-ξ),
=[a1,...,an,b1,...,bm,d1,...dr]T;
误差结构图,
e(k)=w(k)/AD
u(k) 过程 y(k)
B/A
(z-1) (z-1)
( z-1) (k)
分析：
(k)= [y-u]= [(-)u+e]
若e=w，且u与w不相关，则→w意味着：→，→DA
但不能保证→A，→B，→D
针对y=Bu/A+Cw的输出误差法计算式：φ(-,u,)=[-(k-1),...,-(k-n),u(k-1),...,u(k-m),(k-1),...,(k-r)]T;
(k)=y(k)-(k),(k)=φT(-,u,),
=[a1,...,an,b1,...,bm,c1,...cr]T;
误差结构图,
e(k)=Cw(k)
u(k) 过程 y(k)
B/A
(z-1) (k)
( z-1)-1
分析：
(k)= u+e- [u+(-1)(k)]=( -)u+[e- (-1)(k)]
若e=Cw，且u与w不相关，则→w意味着：→,→C
针对y=Bu/A+w/D的输出误差法计算式：φ(-y,u,-ξ)=[-y(k-1),...,-y(k-n),u(k-1),...,u(k-m),-ξ(k-1),...,-ξ(k-r)]T;
ξ(k)=y(k)-φT(-y,u,0),(k)=y(k)-(k),(k)=φT(-y,u,-ξ),
=[a1,...,an,b1,...,bm,d1,...dr]T;
误差结构图,
e(k)=w(k)/D
u(k) 过程 y(k)
B/A
(z-1) ( z-1) (k)
分析：
(k)= [y-u]= [(-)u+e]
若e=w，且u与w不相关，则→w意味着：→，→D
递推算法的实现
UD分解算式P(t)=，t=t-1+P(t)对舍入误差敏感，易导致病态。令P(t)=U(t)D(t)UT(t),P(t-1)=U(t-1)D(t-1)UT(t-1)。定义：d=dim(θ)
f=UT(t-1)ψ，g= D(t-1)f，β=λ(t)I+ψTP(t-1)ψΛ=λ(t)I+fTgΛ=λI+[]Λ
P(t)=[U(t-1)D(t-1)UT(t-1)- U(t-1)D(t-1)UT(t-1)ψΛβ-1ψTU(t-1)D(t-1)UT(t-1)]
=U(t-1)[D(t-1)-gΛβ-1gT]UT(t-1)
D(t-1)-gΛβ-1gT===+
=[ D(t-1)-gΛβ-1gT]-=-+DdededT-gΛβ-1gT
定义：Md=-DdededT+gΛβ-1gT，注意到Md的最后一行和最后一列全为零，因而有：
=Dd-VTddΛβ-1Vdd，=1，=- VTdiΛβ-1Vdd，其中：
g=Vd==。定义βk=λI+[]Λ,βd=λI+[]Λ=β，注意到gi=Difi=VTid,所以：βk+1=βk+gk+1Tgk+1Λ=βk+ Vdk+1VTdk+1Λ，或者写成：
βk-1=βk+1-1+βk-1Vdk+1VTdk+1Λβk+1-1，各项右乘Vdk+1得：
βk-1Vdk+1=βk+1-1Vdk+1+βk-1VdiVTdiΛβk+1-1Vdk+1=βk+1-1Vdk+1+βk-1Vdk+1[Dk+1-]，于是有：βk+1-1=βk-1。定义：
Vd-1=，我们有：=-Vd-1Λβ-1Vdd+，因而,
Md=Vd-1Λβd-1VTd-1+Vd-1Λβd-1VddVTddβd-1ΛVTd-1=Vd-1Λ[βd-1+βd-1-1VddVTddβd-1Λ]VTd-1
= Vd-1Λβd-1-1VTd-1。至此，我们把化成与相同的形式，即=-Vd-1Λβd-1-1VTd-1，
=D(t-1)-gΛβ-1gT=-VdΛβd-1VTd，仿此，可求出所有的和：
=Diβi-1βi-1,=-giΛβj-1-1fTj。以下是UD分解的计算流程：
初始化U(0)、D(0)，t=0，
1 t=t+1；计算f=UT(t-1)ψ，g= D(t-1)f，β0=λ(t)I；
2 对于j=1,...,d=dim(θ)，计算下列3-5步；
3 计算：βj=βj-1+fTjgjΛ；
Djj(t)= Djj(t-1)βj-1βj-1；
Vj=gj；
Qj=-Λβj-1-1fTj；
4 对于i=1,...,j-1，计算第5步(j=1时跳过第5步)；
5 计算：Uij(t)=Uij(t-1)+ViQj；
Vi=Vi+Uij(t-1)Vj；
6 计算：P(t)ψ=βd-1。
稳定性判断递推算法失败的一个原因是模型在修正过程中变为不稳定模型。解决的办法是在每次修正后判断模型的稳定性，若模型稳定，则接受修正，否则，将模型的修正量减半。
Jurry-AstrOm稳定判据
设线性定常离散系统的特征方程为
A(z)=a0zn+a1zn-1+......+an-1z+an=0
按下表建立Jurry表：
a0 a1 a2,....,an-1 an
-) an an-1 an-2,....,a1 a0 x an/a0
b0 b1 b2,.....,bn-1
-) bn-1 bn-2 bn-3,.......,b0 x bn-1/b0
,...........................................................
,..........................................................,
r0 r1
-) r1 r0 x r1/r0
q0
特征方程稳定的充分必要条件是：a0 >0,b0 >0,......,r0 >0,q0 >0
以上判据很容易编制成程序。
递推算法流程
1 指定模型结构，d=dim(θ);
2 初始化P(0)、θ(0)、ψ(0)、φ(0),L(0)= P(0)ψ(0),k=1，指定Λ及λ(0)
3 计算 e=y-φTθk-1,置参数修改步长Amy=1
4 预计算新参数θk-1=θk-1+LΛe×Amy
5 测试分母多项式是否稳定
6 若不稳定，令Amy=Amy/2，返回4
7 计算残差ε= y-φTθk
8 修正φ、ψ（包含需要计算的输出估计、辅助变量、滤波等）
9 计算增益向量L，修正P，λ=aλ+b
10 若数据用完，则结束；否则k=k+1,返回3。
其它
P(0)、θ(0)的初始值：P(0) =107I，θ(0)=0，P(0)充分大是为了保证有足够大的参数搜索区域，但θ(0)的选择还是要尽可能地接近真实参数。
数据复用：为了解决递推算法中信息利用不充分的问题，可采取数据复用的办法，即在数据用完后，以所得的参数估计作为下一轮递推估计的初始值，再次执行递推程序。新一轮的递推程序中，P(0)可选为P(0)=10 P(N)。
结果求期望：θ=
遗忘因子λ(t)：遗忘因子的作用有两个，一是对定常系统的辨识，在递推估计的初始阶段，λ(t)<1可加快参数向真值的运动，但在递推估计的结束阶段应保持λ(t)=1，否则会影响参数的收敛精度。二是对时变系统，遗忘因子可使参数估计跟随实际参数的变化，但要注意：1 参数变化频率应远小于系统响应频率；2 λ(t)的大小应适应参数变化频率的要求，λ(t)小，参数变化剧烈。
λ(t)的递推式通常为：λ(t)=aλ(t-1)+b，其中a、b应满足b<=1-a。λ(∞)=b/(1-a)。
参数a、b及λ(0)的选取决定了λ(t)的运动轨迹。在递推辨识中，通常选取λ(0)>0.8；a的大小控制着λ(t)的收敛速率，对于初始值和终值均相同的两条轨迹，a大则收敛慢，a小则收敛快。
第六讲 结构辨识
Hanker矩阵与积矩矩阵
Hanker矩阵：设g(1),g(2),.....g(L)为脉冲响应序列，定义Hanker矩阵为：
H(j,k)=，j>n0。Rank(H(j,k))= n0，任意k。对于有噪声的情况，可令Rg(k)为脉冲响应序列的自相关函数，用Rg(k)替代g(k)构造H(j,k)即可。
积矩矩阵：考虑差分方程：A()y(k)=B()u(k)+ w(k)，定义：
H*n=，L>n,Rank(H*n)=min(n+n0,2n)。因为U部分总是满秩。定义：Hn=H*Tn H*n，DR(n)=det(Hn)/ det(Hn+1)，若DR(n)较DR(n-1)显著增加，则n比n-1更接近n0，应选n。
残差方差
考虑伪线性回归模型：y(k)=φT(k)θ+w(k)，n0=dim(θ)。模型的输出残差为：
=y(k)-φT(k)=+w(k)，其中，=φT(k)θ-φT(k)，n=dim()。残差方差为：
V=T=(T++)，因为w(k)为白噪声，所以当L→∞时：
V=T+ **
其中T，因此，V=，n≥n0。
在实际计算中，V(n)的变化不可能在n0处有明显的拐点。
F检验法：
t=~F(n2-n1,L-n2)，若设定风险水平为α=0.05，则Fα(2,100)=3.09。据此，我们可通过假设检验来判断V(n)的变化：当n2-n1=2，L-n2>100时，若t>3.09，则V(n2)是显著大于V(n1)的，n2比n1更接近n0；若t<3.09，则V(n2)不显著大于V(n1)，因此，可取n0=n1。
残差自相关性模型的残差带有两方面的信息，一是参数误差信息，二是附加噪声信息。**式也告诉我们，当模型阶次大于等于实际系统阶次时，残差中主要包含附加噪声信息。因此，通过判断残差与附加噪声的统计性质是否一致，也可得出模型是否具有合适的阶次。
当附加噪声为白噪声时，可用残差的自相关函数来判断它的白色性。常用方法为：
设R(0),R(1),...为残差的自相关函数，ρ(k)=R(k)/R(0)，若：E[ε]=≈0，且|ρ(k)|≤，(k可取1~20)，则可认为残差为白噪声序列。
最终预报误差在辨识过程中，我们只能得到有限次(L)观测数据，由此得到的模型L对未来输出的估计误差如何是我们所关心的。这是用最终预报误差最小来确定模型阶次的出发点。记为根据模型L所得的输出预报。残差方差的有限数据估计为：=，最终预报误差可表示为：==E{[z(k)- ]2}，而：
E{[z(k)- ]2}=E{[hT(k)+w(k)]2}= E{T h(k)hT(k) }+=+
=Tr[(HHT)T] +
因为=（(HHT)-1Hw，所以，=dim(L)+。又因为=，
最后有：=。
AIC准则在极大似然原理中，我们通过极大化p(ZL|θ)来求取参数的极大似然估计，其中ZL是已经得到的有限次观测数据。设系统的真实参数为θ0，由极大似然原理，我们有以下关系式：
，+，
→，
2=2+2{+}
2-2=～χ2(n)，n=dim()，两边求期望，则有：
E{2}- E{2}=n。
如果考虑随机变量z，并定义极大似然估计的信息测度：I(θ0,)=E{}≥0，那么，使I(θ0,)最小的模型结构就是我们所要求的。我们有：
I(θ0,)=E{}=S(θ0)-S()=S(θ0)-E{logp(z|)}，将S()在θ0附近展开：
S()=S(θ0)+E{ +T+............}
因为S(θ)的最大值出现在θ0，所于=0，又因为=-M(Fisher信息矩阵)，因而：I(θ0,)= E{TM}，其中TM～χ2(n)，n=dim()，或者：E{2}- E{2}=n，注意到E{2}= E{2}，可得：
-E{2}=- E{2}+2n=AIC。现在，使I(θ0,)最小转变为使AIC最小。
在正态分布的假定下，-=+Llnσ+（L/2）ln2π，其中
→L/2，所以，AIC= Llnσ2+2n
逐步回归法对于一般的回归方程：y=，可以认为y是由xi构成的，但不同的xi对y的贡献是不同的。方程右边不包含变量z，是因为z对y的贡献为零。这一简单思想可用来确定模型结构。现在的问题是如何定义“贡献”及如何计算“贡献”。
定义：S总=，总体偏差平方和
S回=，回归平方和，
S剩=，剩余平方和三者之间的关系为：S总= S回+ S剩+
假定1 及均为零；2 =0，=0；
则S总= S回+ S剩，并且S剩== S总-，即S回=。
类似地，当我们在a1，a2，....，an中选取p个元来构造回归方程时，定义：=。
下面我们观察一下选元过程。定义：
XA==0，令C=XTX，则CA=0。
C==C0=
第一次消元，p=1：假定对第一列消元，则得到的第一个回归系数为：==，消元以后，C0转变为C1，其各元值为：c1(1,j)= c0(1,j)/c0(1,1)；
c1(i,j)=c0(i,j)- c0(i,1) c0(1,j)/c0(1,1)= c0(i,j)- c0(i,1) c1(1,j)。
c1(i,n+1)==Ei- c0(i,1) / c0(1,1)E1= Ei- c1(1,j) E1，i=2,...,n。
c1(n+1,j) ==Ej- c0(1,j)/c0(1,1) E1= Ej- c1(1,j) E1，j=2,...,n。
c1(n+1,n+1)= S总-E1= S1剩。
第二次消元，p=2：假定对第二列消元，则得到的第二个回归系数为：
==-E1=。此时，第一个回归系数改变为：=- c1(1,2)。消元以后，C1转变为C2，其各元值为：
c2(2,j)= c1(2,j)/c1(2,2)；j=3,...n；
c2(i,j)=c1(i,j)- c1(i,2) c1(2,j)/c1(2,2)。定义h22=1/ c1(2,2)，
c2(i,n+1)==- c1(i,2) / c1(2,2) =- c2(2,i)，i=3,...,n。
C2(n+1,n+1) = S1剩- c1(n+1,2) = S总-E1- [E2- c1(1,2) E1]= S总-E1- E2= S2剩。
一般地，当增加第p个参数时，原来已得到的参数的修正公式为：
=- cp-1(j,i)，也即：=- cp-1(j,i)。新补参数的计算公式为：
==-
定义偏回归平方和Qi：当增补参数和相应的xi到回归方程时，定义Qi=-为xi的偏回归平方和，以此来衡量xi对回归方程的贡献。
Qi=-=-
=--=Ei-
=cp-1(i,i)[ -]=()2 cp-1(i,i)=
= Cp(n+1,n+1)= Cp-1(n+1,n+1)- cp-1(n+1,i)= Cp-1(n+1,n+1)- Qi
注意Qi的计算可在第p次消元之前进行，这就允许我们在候选元中选择对回归方程的贡献最大的元。
为避免过估计，应当进行显著性检验。已经有人证明了：
统计量F= Qi/i(N-p-1)～Fα(1,N-p-1)，当置信度设为95%时，若N-p-1>150，且F>3.9时，可认为xi对回归方程的贡献是显著的。
将逐步回归法应用于结构辨识时，应注意：
1 设定的初始结构要大于系统的时延与阶次的和；
2 为了得到符合实际的模型，应保证各多项式的连续性（即若多项式的阶次为n，则指数小于n的参量的系数不应为零）。为此，对候补元的确定规定如下：
开始时，p=[-y(k-1),全部的uj]；在一次消元后，若已选中的元有[-y(k-1),...,-y(k-i),u(k-j),u(k-j-1),...,u(k-r)]，则后选元为p=[-y(k-i-1),u(k-j+1),u(k-r-1)]。
逐步回归法算法流程：
1 读入数据；
2 建立C0，令p=0；
3 p=p+1，确定候补元集合；
4 对候补元集合中的每个元计算Q值；
5 挑出Q值最大的元；
6 计算及F= Qi/i(N-p-1)，如果F>3.9，转7；否则转8；
7 进行第p次消元，转3
8 显示模型结构及参数，辨识结束。
系统辨识大作业
已知一系统为两输入单输出系统，观测数据受有色噪声污染，噪信比为N/S=0.1。系统经2000次采样，存放于文件T3T.TXT中。系统输入u1为7级M序列，u2为u1的63步移位序列。模型类可选为：A(q-1)y(k)=B1(q)u1(k)+ B2(q)u2(k)+w(k)/C(q-1)。要求编制程序，辨识出该模型的结构及参数。
作业文档要求：
描述问题；
选择辨识方法并简单说明所选方法中的结构辨识原理和参数估计原理；
程序流程图及程序清单；
说明程序中用到的一些技术，如数据标准化、UD分解、稳定性判断等；
结构搜索路线及各结构下的参数、残差；
给出最终结果：A(q-1)=
B1(q)=
B2(q)=
C(q-1)=
并给出选择此最终结果的理由；
用你的辨识结果来预报系统输出误差e(k)=y(k)-(k)，并画出e(101)-e(400)的曲线图。
课时：40
参考书：
1 徐南荣、宋文忠、夏安邦 系统辨识，东南大学出版社，1991
2 方崇智、萧德云 过程辨识，清华大学出版社，1988
3 Ljung L,And Soderstrom T,Theory and Practics of Recursive Identification,MIT.Press,Cambridge
4 [美]P.艾克霍夫 系统辨识-参数和状态估计，潘科炎等译，科学出版社，1977
5 [美]夏天长 系统辨识—最小二乘，熊光楞、李芳芸译，清华大学出版社，1983
期刊：
Automatica
Proc,IFAC Identification and System Parameter Estimation
IEEE Trans,On Automatic Control
自动化学报控制理论与应用
学习要求
1 培养独立学习一门新课程的能力，为今后学习和研究打下基础（要求大家尽量少依赖听课，多自学）。
2 掌握基本的辨识理论和辨识技术
3 能独立设计辨识实验，并编程计算
4 学习一些现代建模技术
考核办法
给出一个数据文件，通过编程对其进行辨识，并写出报告
第一讲 概论
实体与模型
实体：客观存在的事物及其运动状态，有时也称之为“系统”
模型：实体的一种简化描述。模型保持实体的一部分特征，而将其它特征忽略或者变化。不同的简化方法得到不同的模型。
模型分类：
直觉模型：地图、建筑模型、照片、软件演示文档等
物理模型：风洞、水力学模型、传热学模型、电力系统动态模拟模型等。（缩小的复制品）
数学模型：描述实体中一些关系和特征的数据模型。例如：投入/产出模型、热源与室温的关系模型等。
数学模型
数学模型还可分为：
图表模型：如阶跃响应、脉冲响应、频率响应、温度与热电偶输出关系表解析模型：代数方程、微分方程、差分方程、状态方程程序模型：神经网络仿真程序语言模型：模糊关系模型
获得数学模型的方法有：
经验总结法：模糊关系模型、静态线性关系模型
机理分析法：解析模型
实验法,图表模型
数据拟合法：解析模型、程序模型
用数据拟合法获得解析模型的过程即为系统辨识。
系统辨识
L,A,Zadeh[1962],辨识就是在输入和输出数据的基础上，从一组给定的模型类中，确定一个与所观测系统等价的模型。
这一定义给出了系统辨识的三要素：数据、模型类和准则。
数据：由观测实体而得。不唯一，受观测时间、观测目的、观测手段等影响。
模型类：规定了模型的形式。不唯一，受辨识目的、辨识方法等影响。
准则：规定了模型与实体等价的评判标准。不唯一，受辨识目的、辨识方法等影响。
系统辨识的三要素是评判数据拟合方法优劣的必要条件，只有在相同的三要素下，才可区分数据拟合方法的优劣；而在不同的三要素下，这种结论也会改变。（图1）
辨识目的明确模型应用的最终目的是很重要的，因为它将决定如何观测数据、如何选择三要素以及采用什么数据拟合方法等。而最根本的是它将影响辨识结果。
辨识目的主要取决于模型的应用。辨识模型应用有以下几个方面：
辨识目的 选择 准则
模型类
辨识实体 数据观测 数据 数据拟合 模型
图 1 系统辨识三要素
1 验证理论模型；要求：零极点、结构（阶次及时延）、参数都准确；模型类同理论模型。
2 设计常规控制器；要求：动态响应特性、零极点、时延准确；便于分析的模型类。
3 设计数字控制器；要求：动态响应特性、时延准确；便于计算机运算的模型类。
4 设计仿真/训练系统；要求：动态响应特性准确；便于模拟实现的模型类。
5 预报预测；要求：动态响应特性、时延准确；便于计算机运算的模型类。
6 监视过程参数，实现故障诊断；要求：参数准确；能直观体现被监视过程参数的模型类。
系统的定量与定性分析；要求：静态关系准确；模型简单，便于人脑判断。
辨识的一般步骤
我们将结合一个实际例子来说明辨识的一般步骤。
例：
上图为长网造纸的流程简图。对于造纸企业来说，质量控制就是要控制好成品纸的定量与水份。而纸的定量与水份与纸浆浓度D、纸浆流量F、车速V及蒸汽压力P都有关系：
G=f(D,F,V,P)
W=g(D,F,V,P)
为了采用计算机对上述过程进行控制，需要建立数学模型。
这就是系统辨识的第一步：明确辨识目的。
为了实现这个目标，我们要进一步了解系统的一些特点。
这是第二步：收集先验知识。
经过现场调查，我们发现：
1 车速调整存在同步困难，而不同步会引起断纸，因此，通常将车速设为恒定；
2 流量的改变到定量的改变存在约60秒的延迟，而响应过程只有约2秒；（保持浓度不变）
3 浓度的改变到定量的改变存在约120秒的延迟，而响应过程约80秒；（保持流量不变）
4 蒸汽压力的改变到水份的改变存在约45秒的延迟，而响应过程约60秒；（保持浓度与流量不变）
第三步：设计辨识试验。
辨识试验的目的是使采集到的数据能反映系统的动态特性，因此要对系统进行分块，设计对分块后系统施加的激励信号，设计数据采集时的采样频率。对于本例，其中的一个分块为流量和蒸汽压力对定量、水份的影响；
试验时，保持车速和纸浆浓度不变；
对流量和蒸汽压力，分别施加伪随机序列扰动，幅度以不引起断纸为限；
设定采样频率为2，试验时间为1000秒。
采集信号为：定量、水份、纸浆流量和蒸汽压力第四步：现场准备。（接线图）
现场准备要做以下几件事：
向企业领导申请试验时段；
准备扰动信号发生器，并通过预发信号，检验扰动信号是否准确；
测验现场信号的干扰情况，必要时设计模拟信号滤波器；
准备模数转换设备，调好信号的零迁和放大参数；
现场接线，将生产设备、试验设备与计算机连接。
第五步：数据采集。
将采集到的数据存盘，并编写数据说明文件；
第六步：数据预处理。
对采集到的原始数据进行变送器非线性校正、数字滤波、标准化、重抽样等加工，使数据适合辨识工具的处理，同时也应满足模型要求。
以上步骤为数据观测过程。
第七步：选择模型类。
选择模型类的工作有两部分：其一是选择应用模型，通常应依据辨识目的来选择；其二是选择参考模型，参考模型是便于进行结构辨识和参数估计的模型第八步：结构辨识与参数估计。
应用辨识理论和方法编制程序，对第六步所得的数据进行拟合，得到参考模型的阶次和参数。
第九步：模型检验。
对所得到的参考模型按评判准则进行检验，如不达要求，则分析问题所在，并返回到前期各相应步骤。
第十步：模型转换。
将参考模型转换为应用模型。
第十一步：应用评价。
从应用角度评价模型，如不符合应用要求，应分析问题所在，并返回到相应步骤。
下图描述了辨识各步骤之间的关系。
思考题：你认为系统辨识还有用吗？
明确辨识目的 收集先验知识
设计辨识试验
现场准备
数据采集
数据预处理
选择模型类
结构辨识
参数估计
不通过
模型检验 分析问题
通过
模型转换
不成功
应用评价 分析问题
成功
最终模型
系统辨识一般流程第二讲 辨识三要素一、数据本节介绍辨识数据的特点及获得适宜辨识的数据的方法。
随机过程X(t):在每一个时间点（t0）上，都是一个随机变量，其概率密度函数p(x,t)随时间变化。
平稳随机过程：在所有时间点上，概率分布都相同的随机过程，其概率密度函数p(x)不随时间变化。
各态遍历平稳随机过程：从整个时间轴上看，每个随机事件都会发生的平稳随机过程。其谱密度函数与概率密度函数类似。时间平均等于集合平均。
数字特征特征
随机过程
平稳随机过程
各态遍历平稳随机过程
均值
(期望值)
均方值
方差
=
自相关函数
协方差函数
=
=
=
一些关系
=
互相关函数
=
=
互协方差函数
=
=
其它
说明：离散计算时假设采样时间间隔为To，则时延τ=*To。
相关函数的性质
1 ≥0
2
3
4 若x(t)是周期为T的信号，则其自相关函数也是周期为T的信号。即：
5 若x(t)=y(t)+z(t)，且y(t)与z(t)互不相关（），则
6若x(t)=y(t)+z，其中，z是一个常数，则
7若x(t)均值为零，且不含有周期性成分，则当τ很大时，x(t)与x(t+τ)必然是互相独立的（不相关），因此，,τ充分大。
8 若x(t)均值为零，则。这是因为在通常情况下，等于向下平移，因此，当=0时，两者相等。
9 对于线性系统y(k)=G(z)u(k)，有
10 虽然x(t)是个随机过程，但却不是随机过程，而是一个确定性的时间函数。
Parseval定理与功率谱
Parseval定理：确定性信号x(t)的总能量为：
确定性信号x(t)的平均功率：
确定性信号x(t)的平均谱密度：
随机性信号x(t)的平均谱密度：
维纳—肯塔金关系式：
随机过程x(t)的谱密度与自相关函数构成一组傅立叶变换对：
定义互谱密度为互相关函数的傅立叶变换：
应用维纳—肯塔金关系式，可以证明，对于频率响应为的线性系统，在随机输入下的输出谱密度和互谱密度分别为：
输出谱密度关系告诉我们：要充分激励系统，就要使输入信号的频谱“宽”于系统频谱。
白噪声
如果一个零均值、平稳随机过程的谱密度为常数，我们称之为白噪声（由白色光联想而得）。白噪声有以下特点：
1
2 ，频谱宽度无限。
∞，τ=0
0，τ≠0
3 ，其中，为Dirac函数，即=
且
4 无记忆性，即t时刻的数值与t时刻以前的过去值无关，也不影响t时刻以后的将来值。从另一意义上说，即不同时刻的随机信号互不相关。
白噪声的用途：
1 作为系统输入时，有，τ=0,1,2,…，即为系统的单位脉冲响应。
2 作为被辨识系统输入时，可以激发系统的所有模态，可对系统充分激励；
3 作为被辨识系统输入时，可防止数据病态，保证辨识精度。
4 在辨识过程中，以输出估计误差是否具有白色性来判断辨识方法的优劣，也可用来判断模型的结构和参数是否合适。
5 产生有色噪声。
白噪声的产生方法：
1 （0，1）均匀分布白噪声：{，i=1,2,3…}，初值可取为：
2 正态分布白噪声：，
其中ξ为服从（0，1）均匀分布的白噪声。
有色噪声有色噪声是指每一时刻的噪声和另一时刻的噪声相关，因而其谱密度也不再是常数。在工业生产实际中，白噪声在物理上是不存在的，常见的往往是有色噪声。
有色噪声的表示定理：设平稳噪声序列{e(k)}的谱密度是ω的实函数，则必定存在一个渐近稳定的线性环节，使得在输入为白噪声序列的情况下，环节的输出是谱密度为的平稳噪声序列{e(k)}。
白噪声 线性环节 有色噪声
(成形滤波器)
{w(k)} H(z-1 ) {e(k)}
M序列（二位式最大长度伪随机序列）
例：4阶M序列 1 1 1 1 0 1 0 1 1 0 0 1 0 0 0 1 1 1 1 。。。。。
1 1 1 1 –1 1 –1 1 1 –1 –1 1 –1 –1 –1 1 1 1 1 。。。。。
特点：
周期性，周期长度为N=2n-1 (n为阶次)，是n个移位寄存器所能表示的最多状态数
M序列中某种状态连续出现的段称为游程。一个周期中有2n-1个游程，游程长度为1 ——n，但出现的概率是随机的。长度为1的游程有2n-2个，长度为2的游程有2n-3个，长度为3的游程有2n-4个，以此类推，最后，长度为n的游程有1个。
一个周期中1的个数比0的个数多1（例中有8个1，7个0）
若以-1代替0，则序列的自相关函数RM(τ)为
RM(τ)=1 τ=i*N,i=0,1,2,…
RM(τ)=1/N τ=i*N+k,i=0,1,2,…,k=1,2,…N-1
当N充分大时，M序列M1：{M(k)}与它的J步移位序列M2：{M(k+J)}在τ=0，1，…J-1和τ=J+1，…N时是不相关的，即
M序列的频谱：设为一步M序列信号的持续时间，N是一个周期的持续时间。M序列的频谱为：
其中：
上式告诉我们：
1 M序列的频谱不是光滑的曲线，而是线条谱。
2 M序列的直流分量（）与N2成正比，因此，加大N，可减少M序列中的直流分量。
3 M序列的频带为，因此，减少，可增加带宽。
4 谱线密度与N成正比，
应用,
将M序列作为扰动信号有以下好处：
1 幅值可取a,-a，容易选择，且当N充分大时，均值约等于0
2 在一个周期内，自相关函数RM(τ)近似为δ函数，因此，以M序列为输入的线性系统，其互相关函数序列等于脉冲响应序列（N大于过渡过程）
3 对于多输入单输出系统，可将同一M序列的不同移位序列（例：{M(k)}、{M(k+J)}、{M(k+2J)}等）作为各输入信号的扰动信号，当输出对各输入的过渡过程小于J时，可认为输入信号之间是互相正交的。
逆重复M序列：
将M序列与方波序列相乘，得逆重复M序列，其特点为：
周期=2N，一个周期中1的个数与-1的个数相等，自相关函数RM(τ)为
RM(τ)=1 τ=i*N*2,i=0,1,2,…
RM(τ)=-1 τ=i*N*2+1,i=0,1,2,…
RM(τ)=0 τ=i*N+k,i=0,1,2,…,k=1,2,…N-1
用作辨识系统的扰动信号的优点在于：具有白噪声的优点，而其幅值可控。
辨识试验设计设计原则 在安全的前提下，尽可能地激励系统；保持输入输出关系；适当解耦
明确目的与要求 模型用途了解辨识对象
划清要辨识系统的边界，选好输入/输出，从边界外连入的其它信号尽量保持稳定，并作为被辨识系统的噪声。
整体/局部 ——确定哪些输入需要叠加扰动信号，哪些输入要保持稳定输入/输出/噪声——确定是多输入还是单输入（耦合关系）、确定过渡过程是否有明显差异（时间常数）、了解噪信比的大小（滤波）及噪声类型（白色、有色）
值域范围——确定信号采集时是否需要零迁、放大安全工况——可叠加扰动信号的类型与幅值选择工况 生产负荷、试验时间、系统区域隔离、地理区域隔离、安全措施扰动信号设计要点,扰动信号频带应宽于过程的工作频带；持续时间为3-5个扰动信号周期；幅值由安全工况确定，对可中断生产的系统，且试验不引起原材料浪费的，可实施单独试验。
扰动信号类型——M序列/白噪声序列扰动信号幅值——由安全工况确定
M序列设计——级数n、步宽Ts、正交化
设系统过渡过程时间为T，最高工作频率为fmax，（通常，fmax=(3~15)/T）则以1.5*T< (2n-1)*Ts和1/Ts>3* fmax为原则，选择n 和Ts。
实际应用过程：
1 确认系统过渡过程T;
2 选择最高工作频率fmax=(3~15)/T；（设为fmax=k/T）
3 根据1/Ts>3* fmax，得出Ts=T/(3k)；
4 根据1.5*T< (2n-1)*Ts=(2n-1)*T/(3k)或4.5k=2n-1，求出n；
另一种设计方法：
1 确认系统过渡过程T，选择n；
2 根据n Ts=(0.1~0.8)T，算出Ts=(0.1~0.8)T/n。
对于多输入单输出系统，要求(2n-1-1)*Ts>1.5*∑(输出对各输入的过渡过程时间）。
白噪声序列——1/Ts>3* fmax
数据采集设计采样间隔——满足香农定理，设采样间隔为T0，1/T0>2* fmax。对于有扰动信号的系统，T0< Ts，且Ts应是T0的整数倍。
实验时间——数据长度：3-5个扰动信号周期，预实验1-2个扰动信号周期物理信号处理——信号零点选择、滤波器选择、放大器选择
（现场准备扰动信号发生器——测试、加限幅器、选择叠加方式、接线信号采集装置——测试、时间间隔标定、接线）
试验设计报告 上述内容及人员组织，报领导批准
数据预处理传感器非线性校正——根据传感器生产厂提供的静态非线性曲线，编制程序对相应的采样数据进行校正。
滤波——对于输入输出加相同的滤波器，输入输出关系不变重抽样——对于过渡过程较长的子系统，可用重抽样的方法，重新选择采样间隔，以降低模型阶次标准化——零均值、单位方差目的：降低舍入误差，防止病态方程，提高计算精度方法,y(k)+ay(k-1)=bu(k-1)+cu(k-2)+d-------------------(1)
∑y(k)/N+ ∑ay(k-1)/N= ∑bu(k-1)/N+ ∑cu(k-2)/N+d
=∑y(k)/N= ∑y(k-1)/N
=∑u(k-1)/N= ∑u(k-2)/N,d=(1+a) -(b+c)
Y(k)= y(k)- ,U(k)=u(k)-
Y(k)+aY(k-1)=bU(k-1)+cU(k-2)----------------------------------(2)
σ2y=∑Y(k)2/N,σ2u=∑U(k)2/N
Yσ(k)=Y(k)/ σy,Uσ(k)=U(k)/ σu
Yσ(k)+a Yσ(k-1)=b Uσ(k-1)* σu /σy +c Uσ(k-2) * σu /σy
B= b * σu /σy,C= c* σu /σy
Yσ(k)+a Yσ(k-1)=B Uσ(k-1)+C Uσ(k-2)-----------------------(3)
b= B*σy/ σu,c= C*σy/ σu ----反演公式二、辨识准则误差定义参数误差,J=∑/n,，适合于算法研究
方程误差,J=∑[e(k)]2/N e(k)= y(k)+y(k-1)-u(k-1)- u(k-2)，J是参数的线性函数；只能作一步预测；适合模型为：AY=B1U1+B2U2
输出误差,J=∑[e(k)]2/N,e(k)= y(k)-(k),(k)= -(k-1)+u(k-1)+u(k-2)，J是参数的非线性函数；当输入序列已知时，可作多步预测；适合模型为：
u(k) 过程 y(k) u(k) 过程 y(k)
B/A B/A
-1 /
e(k)
白化误差加权准则,J=∑[e(k) ]TW[e(k) ] /N
模型结构对误差的影响模型阶次越高，J值越小，导致模型过估计采样间隔对模型结构及误差的影响采样间隔增大一倍，模型阶次大致减少一半；
要进行不同采样步长下模型的比较，误差公式应改为：
J=∑[e(k) ]TW[e(k) ]*T02 /N,T0为采样步长三、数学模型对象特性静态模型 输入输出之间的关系不随时间变化动态模型 未来的输入输出关系受过去及现在的状态影响，模型的有效性不随时间变化动态建模 模型的有效性随时间变化，需要经常改变
集中参数/分布参数 集中参数模型只反映变量与时间的关系，用常微分方程描述；分布参数模型反映变量与时空之间的关系，用偏微分方程描述；
序列（权函数）/方程 序列（权函数）反映系统输出在标准输入下的响应；方程反映系统输入输出的关系，给定输入可计算出系统输出；
h(τ),y(t)=∫h(τ)u(t-τ)dτ
y(n)(t)+a1y(n-1)(t)+…+any(t)=b1u(n-1)(t-τ)+…+bmu(t-τ)
=AX+BU Y=CX+DU
线性/非线性 线性关系是：若y1=f(u1),y2=f(u2),则 ay1+by2=f(au1+bu2)
{hn(τ1,…,τn)},y(t)= ∑∫…∫hn(τ1,…,τn)u1(t-τ1)…un(t-τn)dτ1…dτn
=G(X,U)=GLX+GR(X,U) Y=F(X,U)
时变/时不变 模型的结构不边，但参数随时间变化。例如锅炉中的热交换器，由于灰尘的附着而改变热交换的效率
y(n)(t)+a1(t)y(n-1)(t)+…+an(t)y(t)=b1(t)u(n-1)(t-τ)+…+bm(t)u(t-τ)
确定型/随机型 若系统含有随机扰动，则称为随机型的。实际中很少没有随机扰动的，只是强度不同而已。扰动强度的度量用噪信比N/Sy=σW/ σX,(y=x+w)
单变量/多变量 Y(s)=G1(s)U1(s)+ G2(s)U2(s)
离散型/连续型 y(k)= ∑h(i)u(k-i)
A(q-1)y(k)=B(q-1)u(k)+w(k)
Y(z-1)=G(z-1)U(z-1)
X(k+1)=Ax(k)+Bu(k) y(k)=Cx(k)+Du(k)
有时延/无时延Y(s)=e-sτG1(s)U1(s)+ e-sτG2(s)U2(s)
y(k)+a1y(k-1)+…+any(k-n)=b1u(k-τ)+…+bmu(k-m-τ)
常用模型阶跃响应（函数、序列、曲线）
差分方程/微分方程连续传递函数/离散传递函数连续状态方程/离散状态方程一般描述/典型环节惯性环节 k/(1+Ts)n
微分环节 kTs/(1+Ts)n
积分环节 1/Tas(1+T0s)n
惯性+微分环节 k(1+aTs)/(1+Ts)n
模型间的相互转换
差分方程 离散传递函数 离散状态方程
阶跃响应
微分方程 连续传递函数 连续状态方程
差分方程与离散传递函数,A(q-1)y(k)= B(q-1)u(k) Y(z-1)=U(z-1)
微分方程与连续传递函数的关系与此类同。
差分方程、离散传递函数到阶跃响应：以阶跃信号为输入，通过计算机仿真得到阶跃响应。
阶跃响应到连续传递函数：
阶跃响应曲线转换为典型环节传递函数
1设系统传递函数为W(s)，在阶跃输入下，系统响应的拉氏变换为F(s)=W(s)/s，若y(kT) k=0,1,2,…为系统的阶跃响应序列，则F(s)的离散计算公式为F(s)=T∑e-skT y(kT)；
2 根据{y(kT)}的曲线形状，可选择一个典型环节传递函数w(s)，并令G(s)=w(s)/s为其阶跃响应的拉氏变换；
3 Jt=∫|g(t)-y(t)|dt,等价于 Js=|G(s)-F(s)|，或者 J=∑[G(sj)- F(sj)]2，sj充分小
4 对于充分小的s,选择sj>0，j=1,2,…,L（L大于G(sj)中参数个数）,可计算出F(sj)，同时G(sj)（以及J）中只含w(s)的未知参数；
5 令 =0，求出参数，即得G(s)，然后w(s)=s G(s)
离散传递函数与连续传递函数：
1 离散信号的保持
u(t) y(t)
G(s)
u(z)=Z[u(t)]=Z[u(kT)] y(z)=Z[y(t)]=Z[y(kT)]
u(kT) x(t) y*(t)
L(s) G(s)
u(z)=Z[u(kT)] X(s) y*(z)=Z[y*(t)]
G(z)= y*(z)/ u(z)=Z[G(s)L(s)] G(s)=Z-1[G(z)]/L(s)
若要求y*(z)= y(z)，则应满足y*(kT)= y(kT)，进而要求x(t)≈u(t)，这通过适当选择L(s)来实现。
L1：x(t)=u*(t)=，x(s)= =u*(s)，L(s)=1(或T)
当u(t)=δ(t)时，y*(t)= y(t)，故称为脉冲响应不变等效。G(z-1)=Z[G(s)]，G(s)=Z-1[G(z-1)]
L2：（零阶保持器）x(t)=u(kT)，kT≤t＜(k+1) T，L(s)=(1-e-Ts)/s
当u(t)=1(t)时，y*(t)= y(t)，故称为阶跃响应不变等效。G(z-1)=(1-z-1)Z[G(s)/s]，G(s)= Z-1[G(z-1)]/[(1-e-Ts)/s]
L3：x(t)=u(kT)*((k+1)T-t)/T+u((k+1)T)*(t-kT)/T，kT≤t＜(k+1)T，L(s)= (1/T+s)[(1-e-Ts)/s]2
斜坡响应不变等效。
Z变换与逆 z=eTs,s=ln(z)/T：（可用查表法进行）
G(z)=Z[G(s)]=G(s)| s=ln(z)/T，G(s)=Z-1[G(z)]=G(z)| z=eTs
近似变换：向后差分 s=
向前差分 s=
双线性变换 s=，
离散传递函数与离散状态方程
Y(z-1)=U(z-1)= H(z-1)U(z-1) X(k+1)=Ax(k)+Gu(k) y(k)=Cx(k)+Du(k)
从离散状态方程到离散传递函数,H(z-1)=C(zI-A)-1G+D
从离散传递函数到离散状态方程,设A(z-1)=1+a1z-1+......+anz-n,
B(z-1)=b0+b1z-1+......+bnz-n
则,0 1 0........,0 0 g1=b1-a1b0
0 0 1........,0 0 g2=b2- a1 g1-a2b0
A=,.....................................,G=,................................................
0 0 0,.........,0 1 gn-1=bn-1- a1 gn-2-..........- an-2 g1-an-1b0
-an –an-1 –an-2.........-a2 –a1 gn=bn- a1 gn-1-......................- an-1 g1-anb0
C=[1 0,...............0],D=b0
连续状态方程与连续传递函数的关系与此类同。
连续状态方程与离散状态方程:
=AX+BU Y=CX+DU X(k+1)=Fx(k)+Gu(k) y(k)=Hx(k)+Eu(k)
F=eAT≈I+AT，G=A-1(F-I)B，H=C,E=D (阶跃响应不变)
A=lnF/T≈(F-I)/T,B=(F-I)-1AG,C=H,D=E (阶跃响应不变)
辨识中的参考模型
ARMAX模型 A(q-1)y(k)= ∑Bi(q-1)/Fi(q-1)ui(k)+ D (q-1)/C(q-1)e(k)
AR自回归 X外生变量 MA滑动平均线性回归模型 y(k)=-a1y(k-1)-…-any(k-n)+b1u(k-1)+…+bmu(k-m)+ ν(k)
θ=[a1,…,an,b1,…,bm]T
φ=[-y(k-1),…,-y(k-n),u(k-1),…,u(k-m)]T
y(k)= θTφ+ν
伪线性回归模型 当φ中含有不可观测的变量时，通常用其估计值来代替。因此，含有待估计变量的线性回归模型称为伪线性回归模型。
预报误差模型 y=f(Y,U,θ) +ν，v为零均值的新息，f(Y,U,θ)为y的最优估计。
例：
考虑模型,y(k)=-ay(k-1)+bu(k-1)+e(k)+ce(k-1)
若e(k)、e(k-1)...,e(0)都是已知或可估计的，则e(k)的估计为：
(k)= y(k)+ay(k-1)-bu(k-1)- c(k-1)
则k时刻的输出估计为：
(k)= -ay(k-1)+bu(k-1)+ c(k-1) *
=-ay(k-1)+bu(k-1)+ c[y(k-1)+ay(k-2)-bu(k-2)- c(k-2)]
且有：c(k-1)= -acy(k-2)+bcu(k-2)+ cc(k-2)
两式相加得：
(k)+ c(k-1)= (c-a)y(k-1)+bu(k-1) **
此表明：若(0)是最优预报，则(1)、...、(k)都是最优预报。由此也得到预报误差模型为,f()：(c-a)y(k-1)+bu(k-1)- c(k-1)
v：y(k)- (k)= e(k)+ce(k-1)-c[y(k-1)- (k-1)]
==e(k)
或者,y(k)= (c-a)y(k-1)+bu(k-1)- c(k-1)+ e(k)
对于一般的差分模型：A()y(k)=B()u(k)+C()w(k)
其预报误差模型为：y(k)= (k)+ w(k)
其中,w(k)=y(k)-u(k)
(k)=(1-) y(k) +u(k)
请同学们将离散传递函数模型改写成预报误差模型，并用伪线性回归模型的形式表现之。
第三讲 参数估计的批量法最小二乘算法考虑差分方程：A()y(k)=B()u(k)+ w(k)，其中w(k)为白噪声。假定模型的结构已知（n,m），将其写成线性回归模型：
y(k)=-a1y(k-1)-…-any(k-n)+b1u(k-1)+…+bmu(k-m)+ w(k)
θ=[a1,…,an,b1,…,bm]T
φ(k)=[-y(k-1),…,-y(k-n),u(k-1),…,u(k-m)]T
y(k)=φT(k)θ+w(k)
若我们的观测数据可写出N个这样的等式，
YN=Φθ+WN,ΦT=[φ(1），。。。，φ(N)]
ΦT YN=ΦTΦθ+ΦT WN
θ=(ΦTΦ)-1ΦT YN-(ΦTΦ)-1ΦT WN
θLS=(ΦTΦ)-1ΦTYN,条件1：E(ΦT WN)=0 （w(k)白色）
条件2,ΦTΦ可逆以上结果等同于求使：J=Σ(y(k)-φT(k)θ)2最小的θ，因此称为最小二乘算法。
最小二乘算法的MATLAB程序：
clear all; %最小二乘法for MISO
load t2; %读入数据矩阵
tt=t2;
clear t2;
ll=size(tt); %得到矩阵维数200x2
r=ll(2)-1;
m(1)=input('order of A(z)');
tao(1)=0;
for i=1:r
i
m(i+1)=input('order of B(z)');
tao(i+1)=input('timedelay of B(z)');
end
if m(1)<max(m),m(1)=max(m);end
n=m(1)+max(tao);
lll=ll(1)-n;
in1=r+1;
kn=0;
for k=1:in1
for i=1:lll
for j=1:m(k)
jtao=j+tao(k);
ff(i,j+kn)=tt(i+n-jtao,k);
if k==1,ff(i,j+(k-1)*n)=-ff(i,j+(k-1)*n);end
end;end;
kn=kn+m(k);
end;
for i=1:lll
yy(i)=tt(i+n,1);
end;
fa=ff'*ff;
fay=ff'*yy';
zz=inv(fa);
zta=zz*fay
m，tao
加权最小二乘：若对各次观测数据加不同的权，即求使Jα=Σαk(y(k)-φT(k)θ)2最小的θ,则得到参数的加权最小二乘估计：
θLSα=(ΦTΑΦ)-1ΦTΑYN
参数估计的一般性质无偏性
如果E()=0，（=-） 或E()=，则称估计为无偏的。
有效性如果无偏估计满足Cov()=M-1，则称估计为有效的。其中：
称为Fisher信息矩阵，其逆M-1称为Crammer-Rao 下界。在一般情况下，有Crammer-Rao不等式：
Cov()=E{(-)(-)T}≥M-1
对于z=Hθ+w，若噪声w为零均值、协方差阵Σw=E(WWT)的正态分布噪声，即：
W～N(0，Σw)，则输出信号Z～N(E(Hθ)，Σw)，于是，
，因此：
，于是，M=E(HTΣw-1H)，其Crammer-Rao不等式为：Cov()≥E(HTΣw-1H)-1。
有效估计也称为最小方差估计、马尔可夫估计。
一致性若估计为渐近无偏的（E()=0），且Var()=0，则称为一致估计。Var()=E{(-)T(-)}
最小二乘、加权最小二乘估计性质估计性质
最小二乘
加权最小二乘
无偏性
= (ΦTΦ)-1ΦT WN，当:
1 w(k)为零均值独立随机过程（白噪声）时；或者
2 E(ΦT WN)=0时，有：E()=0。无偏条件较严格。
= (ΦTΑΦ)-1ΦTΑWN
无偏条件同最小二乘估计。
有效性
Cov()=E{(ΦTΦ)-1ΦTWNWTNΦ(ΦTΦ)-1}
=(ΦTΦ)-1ΦTΣwΦ(ΦTΦ)-1≥M-1,若W是服从正态分布的独立随机过程，则Σw=I，因而Cov()=I(ΦTΦ)-1= M-1。
Cov()=(ΦTΑΦ)-1ΦTΑΣwΑΦ(ΦTΑΦ)-1，当Α=Σw-1时，Cov()=(ΦTΣw-1Φ)-1= M-1，此时不要求W是服从正态分布的独立随机过程
一致性
若W是服从正态分布的独立随机过程，则Cov()==0,其中收敛于一个有界非奇异矩阵。而Var()只是Cov()的对角线上元的和，因而也有：Var()=0。
同最小二乘估计
噪声方差估计
e= YN-=Φθ+WN-Φ=Φθ+WN-Φ{(ΦTΦ)-1ΦTYN }=WN-Φ(ΦTΦ)-1ΦT WN
=(I-Φ(ΦTΦ)-1ΦT) WN=B WN,因为BT=B,B2=B
所以，E(eTe)=E(WN TB WN)= E(WN T WN- WN TΦ(ΦTΦ)-1ΦT WN)
=N- (TrΦ(ΦTΦ)-1ΦT)= (N-dimθ)
=Σ(y(k)-)2/(N-dimθ)
在有色噪声环境下，最小二乘估计是有偏的。下面的一些算法对最小二乘进行改进。
广义最小二乘考虑差分方程：A()y(k)= B()u(k)+ w(k)，其中w(k)为白噪声。假定模型的结构已知（n,m）。
如果噪声模型C()已知，显然用C()对输入/输出数据进行滤波，则可得到满足最小二乘估计无偏条件的模型：A()(k)= B()(k)+ w(k)，其中：
(k) = C()y(k)，(k) = C()u(k)。在C()未知时，我们可考虑采用迭代估计的方法去求得。
广义最小二乘的计算步骤为：
1 令C0()=1，i=0下标表示迭代次数；ξ0=1000000
2 计算(k) = Ci()y(k)，(k) = Ci()u(k)；i=i+1;
用最小二乘估计A()(k)= B()(k)+ w(k)中的参数；
用估计模型()、()以及各时刻的观测数据，计算出残差(k):
(k)= ()y(k)- ()u(k)
计算ξi=Σ2(k)及ξ=ξi-1-ξi，如果ξ小于一定数，则结束辨识。否则转下一步。
对于噪声模型C()(k)= w(k)，用最小二乘估计出参数，得到更新的Ci()后，返回2。
以上算法的每一次循环中都要进行滤波和两次求逆。下面的算法将在计算工作量上有所改进。
偏倚校正算法
仍考虑差分方程：A()y(k)= B()u(k)+ w(k)，其中w(k)为白噪声。
令(k)= w(k)，则A()y(k)= B()u(k)+ (k)。分别写成回归模型：
Y=Φθ+，=C+W，组合起来有Y=[Φ，]+W，其最小二乘解为：
=，利用分块矩阵求逆公式，有：
=(ΦTΦ)-1ΦTY-(ΦTΦ)-1ΦT
=D-1TMY
M=I-Φ(ΦTΦ)-1ΦT D=TM
须要注意，的求取仍然是一个迭代过程。偏倚校正算法的计算步骤为：
1 令C()=1，用最小二乘法求=θLS=(ΦTΦ)-1ΦTY，并保留Φ、Γ=(ΦTΦ)-1ΦT以及M=I-ΦΓ。
2 计算=Y-Φ，并依据=C+W构造，计算D=TM。
3 计算=D-1TMY，并计算=Γ及=θLS-。
4 若参数已收敛，则结束辨识，否则转2。
以上算法的一次循环中没有滤波，且只有一次求逆。如果将第3步中的计算改为：
=(T)-1T，则还可省去D的计算。（这一改进由夏天长首先给出。）
辅助变量法分析最小二乘法中，在Y=Φθ+W的各项乘上ΦT,然后利用ΦT W的期望值为零，得到参数的无偏估计。受此启发，若在Y=Φθ+的各项乘上ΨT，使其满足以下两个条件：1.ΨT的期望值为零；2,ΨTΦ可逆，则也可得到参数的无偏估计。
下面讨论辅助变量的选取：
设模型为A()y(k)= B()u(k)+ (k)，若u(k)与(k)不相关：
a 选取辅助模型D()z(k)= F()u(k)，用z(k)、u(k)构造Ψ；
b 若系统的纯时延τ已知，则可用u(k-τ)、u(k)构造Ψ；
c 用u(k)、u(k)构造Ψ
d (k)= D()w(k),若D()的阶次n已知，则可用y(k-n)、u(k)构造Ψ；
e 先求出最小二乘解θLS=(ΦTΦ)-1ΦTY,然后依据()z(k)= ()u(k)计算出输出估计z(k)，再用z(k)、u(k)构造Ψ；
多步最小二乘考虑差分方程：A()y(k)= B()u(k)+ w(k)，其中w(k)为白噪声。模型可改写为C()A()y(k)= C()B()u(k)+ w(k)或：
D()y(k)=F()u(k)+w(k)，其中:
D()=C()A()，F()= C()B()---------**。
此模型可用最小二乘解出D()、F()。这是第一步。
第二步可有两种不同方法：
a 解同次幂方程组由**式，可得关系：A()F()=B()D()，两边分别展开，并按的同次幂相等规则，可列出na+nb+nc个方程：
F=Hθ，其中θ=[a1,...ana,b1,...bnb]T，F=[f1,...fnb+nc,0,...0]
0 0 0...........,0 1 0 0...........,0 0
-f1 0 0..........,0 -d1 1 0.............0 0
H= -f2 –f1 0...........,0 -d2 -d1 1.............0 0
,...................................................
,...............................................-d1 1
,.............,-f2 -f1,.......................-d2 –d1
,.............,-f3 -f2...............................
,....................................................
-fnb+nc................................................
0 -fnb+nc...............,-dna+nc.......................
0 -fnb+nc..........,0 -dna+nc.......................
....................................................
0 0..........,-fnb+nc 0........................-dna+nc
na列 nb列用最小二乘可求得θ的无偏估计，即A()、B()。此法也可用来求C()。
B 传递函数等价降阶
由**式，可有,F()/D()=B()/A()，此说明两个传递函数是等价的。对F()/D()施加激励信号{u(k)}，可得输出z(k)= F()/D()u(k)。用最小二乘法处理{ z(k)，u(k)}，选择合适的阶次，可得A()、B()的无偏估计。
相关最小二乘设模型为A()y(k)= B()u(k)+ (k)，若u(k)与(k)不相关，则用u(k-j)乘模型中的各项并求期望，得：A()Ruy(j)= B()Ru(j)，用最小二乘法可得A()、B()的无偏估计。
第四讲 辨识原理
随机逼近法设广义误差e(k)是参数估计值θ的函数，参数辨识问题可通过极小化e(k)的方差来实现。即求参数θ使下列准则函数最小：J(θ)=1/2E{e2(k)}。J(θ)的负梯度为：
=E{e(k)}。如果可求解=0，则可求得参数的估计。但当e(k)的分布未知时，实际上是不可求解的。
在计算数学中，求二次函数的极小值常采用迭代法。首先给出参数的一个估计值，以二次函数在该参数估计值处的负梯度为修正方向，选取适当的步长后，修正参数估计值，直到收敛。仿此，我们有：θ(k+1)= θ(k)+ρ(k)。如果在求时不求期望，则得到一个随机的迭代算法，称之为随机逼近法：
例：考虑线性回归模型：y(k)=φT(k)θ+e(k)，其中e(k)是零均值随机噪声。J(θ)=1/2E{[ y(k)-φT(k)θ]2},=E{φ(k)[ y(k)-φT(k)θ]}。如果假定e(k)是各态遍历的，则梯度为零可用E{φ(k)[ y(k)-φT(k)θ]}=φ(k)[ y(k)-φT(k)θ]=0来代替，由此得到了最小二乘法。
应用随机逼近法，可得：
θ(k+1)= θ(k)+ρ(k)φ(k)[y(k)-φT(k)θ(k)]。
ρ(k)的选取应保证迭代收敛，可选取满足如下条件的ρ(k),
ρ(k)>0且；；，例如ρ(k)=1/k；
若以准则函数的二阶导数（即海赛矩阵）之逆来参与选择修正方向，则称为牛顿法：
θ(k+1)= θ(k)+ρ(k)R(k)-1φ(k)[y(k)-φT(k)θ(k)]，其中：
R(k)==E[φ(k)φT(k)]=R(k-1)+ρ(k)[φ(k)φT(k)-R(k-1)]
牛顿法的优点是收敛速度快。
模型参考自适应辨识方法（简略）
考虑线性回归模型：y(k)=φT(k)θ+e(k)，其中e(k)是零均值随机噪声。以输出估计误差为反馈信号，以PI调节器的方式来修正参数：
θ=θI+θP,θI(k)= θI(k-1)+Pφ(k)ε(k)，θP(k)=Qφ(k)ε(k)
其中P为对称正定矩阵，Q满足P/2+Q>0
ε(k)为广义误差，最简单的取法为ε(k)= y(k)-φT(k)θ(k-1)
极大似然法输出z是一个随机变量，它的概率密度p(z|θ)取决于参数θ。当获得观测序列ZL={z(1),z(2),...,z(L)}T时，由该观测序列组成的联合概率密度p(ZL|θ)应当取得最大值。（当一个随机事件发生了，我们有理由相信，外部条件一定处于使这随机事件发生的概率最大时的状态。）那么，θ的极大似然估计就是使p(ZL|θ)|θML=max的参数估计值。
由于p(ZL|θ)中ZL已知，因而它只是参数θ的函数，故称它为θ的似然函数。有时也记作L(ZL|θ)。极大似然原理可用下列等价的表示方式：
，，，
求解极大似然估计的下一步是要给出p(ZL|θ)的具体描述。
独立观测设z(1),z(2),...,z(L)是一组在独立观测条件下获得的随机序列，即各观测值是互相独立的，则p(ZL|θ)可简化为：p(ZL|θ)= p(z(1)|θ)p(z(2)|θ)...p(z(L)|θ)
= p(z(k)|θ)，其对数似然函数为：L(ZL|θ)=ln p(z(k)|θ)。
设z(k)～N(m,σ)，即：p(z(k)|θ)=，则其负对数似然函数为：
- L(ZL|θ)=+Llnσ+（L/2）ln2π
当σ已知时，准则函数就是：J(θ)=
当σ未知时，可先由min(J(θ))求及Jmin，再由min(- L(ZL|θ))求。
，
二、非独立观测若z(1),z(2),...,z(L)是非独立观测条件下获得的随机序列，即观测值z(k)是在已有观测z(1),z(2),...,z(k)的基础上得到的，则p(ZL|θ)应按条件概率的乘法规则写成：
p(ZL|θ)= p(z(L)|ZL-1,θ)p(ZL-1 |θ)，以此类推，有：
p(ZL|θ)= p(z(k)|Zk-1,θ)。
设z(k)～N(mk,σk)，并取mk为z(k)的条件均值：(k)=E[z(k)| Zk-1]，即：p(z(k)| Zk-1,θ)=，则其负对数似然函数为：- L(ZL|θ)= +（L/2）ln2π。
当σk≡σ已知时，准则函数就是：J(θ)=
当σk≡σ未知时，可先由min(J(θ))求及Jmin，再由min(- L(ZL|θ))求。
例：考虑模型：z=+，其中w～N(0,σ2I),则～N(,σ2I)，
p(|θ)=，
J=lnL(θ)=
，，，，
，i=1,2,...,n
，i=1,2,...,n
，i=1,2,...,n
，i=1,2,...,n
，→
直接求解以上非线性方程组是困难的，如果已有参数的某个估计，并用其构成如下预滤波器：，，，则上面前两个方程组可写成：
，，这等价于辅助变量法。对于噪声模型，我们已有：
，定义：，则w=C/Dv=Cv-(D-1)w，再次构造预滤波器：v*=v/D，w*=w/D,
则，，这也等价于辅助变量法。
预报误差估计法对于模型结构已知的系统，如果我们获得参数的某个估计，则可对输出进行预报，那么用预报误差的大小来衡量参数估计的优劣也是合理的。定义输出估计误差：
e(k)=y(k)-(k)，将收敛于e(k)的协方差阵。通常采用的预报误差有：
J1(θ)=Tr(ΛD(θ)),J2(θ)=log det(D(θ))。可以证明，当e(k)～N(0,Σe)时，若J1(θ)中的Λ取为Λ=LΣe-1，则预报误差估计等价于极大似然估计；
例：考虑模型：z=+。利用已经获得的参数估计和噪声w的到k-1时刻为止的估计，我们有输出估计：，并且可得噪声估计：=，由此可见，以为输出估计误差时，预报误差法与极大似然法的准则函数是等价的。
若记e(k)= (k)，或写成Ce=Dw，那么对Ce的估计可写成Ce=(D-1)w，或者，
于是，我们有输出估计：，那么输出估计误差（也作为噪声估计）：
Bayse估计
参数是与观测值关联的随机变量。在给定的观测数据下，使得验后的条件概率密度p(θ|ZL)最大的参数最有可能是真实参数。根据全概率公式：
p(θ|ZL)= p(z(L)|ZL-1,θ)p(θ|ZL-1)/p(z(L)|ZL-1)
= p(θ|Z0)
其中p(θ|ZL-1)为参数的验前概率密度函数，p(z(L)|ZL-1,θ)为数据的条件概率密度函数，p(z(L)|ZL-1)是与参数无关的项。根据验后条件概率密度p(θ|ZL)，有两种方法可求得参数估计，一是极大验后参数估计法，另一是条件期望参数估计法，它们统称为Bayes方法。
1 极大验后参数估计法参数的极大验后估计必定满足：=0=+
若令上式右边第一项为零，所得的参数估计即为极大似然估计。可见，极大验后参数估计比极大似然估计多考虑了参数的验前概率知识，因而一般情况下极大验后参数估计优于极大似然估计。但是参数的验前概率密度不易获得，因此，极大似然估计比极大验后参数估计用得更普遍。
2 条件期望参数估计法直接以参数的条件期望作为参数估计值是合乎逻辑的。(L)=E{θ|ZL}。可以证明，条件期望参数估计等价于极小化参数估计误差的方差,minE{[θ-(L)]T[θ-(L)]|ZL}
因此也称为最小方差估计。
例,考虑线性回归模型：y(k)=φT(k)θ+e(k)，其中e(k)是零均值随机噪声，
e(k)～N(0,)。设参数θ在数据Z0条件下的验前概率密度为：
p(θ|Z0)=，其中N=dimθ。那么，参数θ在数据ZL-1条件下的验前概率密度为：
p(θ|ZL-1)=。由噪声e(k)～N(0,)可知，
p(z(L)|ZL-1,θ)=，于是有：
log p(θ|ZL)=const--
=const-
若令，则log p(θ|ZL)= const-，显然极大验后参数估计应该是，而从正态分布的形式来看，它也是参数的条件期望估计。
第五讲 参数估计的广义递推算法广义递推算法的一般结构假定对象模型可用伪线性回归方程描述：y(k)=φT(k)θ+w(k)。设是参数θ的一个估计，e(k,)是广义误差，它是参数估计误差的一种描述。定义二次准则函数：
V(t,)=，其中Λ是对误差矢量各分量的加权矩阵，α(k,t)是对k=1,...,t 这个时间串上所有e(k,)的加权序列，满足：
α(k,t)=λ(t)α(k,t-1),k<t
α(t,t)=1，λ(t)≤1
当λ(t)≡1时，α(k,t)≡1。当λ(t)<1时，称为遗忘因子。记在时刻t使V(t,)最小化的参数估计为：t，则应有：。又若已知在时刻t-1处的最优估计t-1，则将V(t,)在t-1处展开成台劳级数：
V(t,)=
V(t,t-1)+[-t-1]+[-t-1]T[-t-1]+o(|-t-1|)
由可得：
+[t-t-1]= o(|t-t-1|)，若d*d矩阵可逆，则：
t=t-1-{}-1+ o(|t-t-1|)
定义记号：=-，则：
=-=--
=λ(t)-
=λ(t)+ -
假定：
1）=0，t-1为时刻t-1处的最优估计
2）=0
3）≈=R(t-1,t-2)
t=t-1+ R-1(t,t-1)，R(t,t-1)= λ(t)R(t-1,t-2)+ *
若定义：P(t)= R-1(t,t-1)，则P-1(t)=λ(t)P-1(t-1) +，
于是P(t)=，t=t-1+P(t)
此为广义递推算法的基本结构。如果再考虑辅助变量算法，则广义递推算法的一般结构为：P(t)=，t=t-1+P(t)
广义误差的选择假定对象模型可用伪线性回归方程描述：y(t)=φT(t)θ+w(t)。为了构造广义误差，我们先来构造输出估计。先定义记号：
t=[，，，]T，，，，)的构造中为广义输出,为输入观测值，为白色噪声估计，为有色噪声估计。对于广义输出，可选的例子有：
=[-y(t-1),-y(t-2),...,-y(t-n)]T,即系统输出的观测值。
=[-(t-1|t-1),-(t-2|t-1),...,-(t-n|t-1)]T,即系统输出的验前估计值。
=[-(t-1|t),-(t-2|t),...,-(t-n|t)]T,即系统输出的验后估计值。
有色噪声估计可选为0或：
=[-(t-1),...,-(t-r)]T,其中(t)=y(t)-，，0，0) t-1,或
=[-(t-1),...,- (t-r)]T,其中(t)=y(t)-，，0，0) t；
白色噪声估计可选为0或验前输出估计误差或验后输出估计误差：
=[(t-1),...,(t-r)]T,其中(k)=y(k)-，，，)t,
=[(t-1),...,(t-r)]T,其中(k)=y(k)-，，，)t-1
广义误差e(k,)通常为验前输出估计误差，即e(k,)=(k)，尽管如此，由于、及的不同选取，可构造的算法多达27种。
递推算法收敛性分析主要有三种分析方法：
ODE 法：将递推算法转化为伴随微分方程，通过研究微分方程的稳定性来判断递推算法的收敛性。
鞅理论：引进随机Lyapunov函数，并用鞅理论分析这个函数的收敛性，进而得到参数估计的收敛性。
间接法：通过研究R(t)t的收敛性，间接得到参数估计的收敛性。
ODE 法的启发式讨论令t=R(t,t-1)，并假定λ(t)=1，则：
t=t-1+ ，=+(-) **
t-t-1=dθ，=dt，-=dR，f(θ)= ,g(θ)= ，则有：
f(θ)，g(θ)-
ODE法的应用例请参考Ljung的著作。
Martingale收敛定理先介绍一些概念：
σ-代数：设S由集合X的某些子集组成，（X本身为S的一个元素），若S对“差”及可列和封闭，则称S为σ-代数。
递增的σ-代数序列：设S={At,t∈T}为σ-代数序列，若At(，(t∈T；
At(As，t≤s，则称S为递增的σ-代数序列。
随机变量序列向σ-代数的转化：对所有的t，若随机变量序列Xt可以从σ-代数序列At观测得到，即E{ Xt | At }= Xt，则称随机变量序列Xt可转化为σ-代数序列At。换言之，包含在At中的信息足以完全确定Xt。
Martingale：若①{ At }为递增的σ-代数序列；②{ Xt }可转化为σ-代数序列At；
③E(|Xt |)<∞；④Xs=E(Xt| As,s<t)，则称序列{ Xt }为对于{ At }的Martingale。
定理：设{ Xt }为对于{ At }的向量Martingale，且对于所有的t，均有||E{ Xt Xt T}||<∞，则存在一具有有界协方差的向量随机变量X，使得Xt→X,(a.s.)，Xt=E(X| At),(t>0。
Martingale收敛定理：设{Tn}、{αn}、{βn}均为非负随机变量序列，它们可可转化为σ-代数的递增序列Fn，使得E[Tn|Fn-1]≤Tn-1-αn-1+βn-1，若 a.s.(几乎处处)，则
Tn→T (a.s.)，且 a.s,
对于算法：
t=t-1+ R-1(t,t-1)，R(t,t-1)=R(t-1,t-2)+ *
定义随机变量：η(t)=e(t)[1-ψTR-1(t)ψ]，若η(t)可表示为：
η(t)=H-1(z-1)[ψT ]+w(t)，=t-θ0
其中w(t)满足,E[w(t)|Ft-1]=0，E[w2(t)|Ft-1]=σ2，e(t)-w(t)∈Ft-1，
Ft-1是由w(0),...,w(t-1),ψ(0),...,ψ(t-1)构成的σ-代数。
则当H-1(z-1)-1/2严格正实时→0 a.s.
以上结论的证明需引用Martingale收敛定理，其中：
Tt=[T R(t)+2αtβt]/t，αt=η(t)- w(t)+ψT/2，βt=-ψT
典型算法
P(t)=，t=t-1+P(t)
针对Ay=Bu+Cw的方程误差法
——RELS递推增广最小二乘计算式：φ(-y,u,)=[-y(k-1),...,-y(k-n),u(k-1),...,u(k-m),(k-1),...,(k-r)]T;
(k)=y(k)-(k),(k)=φT(-y,u,),
=[a1,...,an,b1,...,bm,c1,...cr]T;
误差结构图,
e(k)=Cw(k)/A
u(k) 过程 y(k)
B/A
(z-1) (z-1)
( z-1)-1 (k)
分析：
(k)= y-u-(-1)(k)=(-)u+[e-(-1)(k)]
若e=w，且u与w不相关，则→w意味着：→，→
但不能保证→A，→B，→C
针对Ay=Bu+Cw的输出误差法计算式：φ(-,u,)=[-(k-1),...,-(k-n),u(k-1),...,u(k-m),(k-1),...,(k-r)]T;
(k)=y(k)-(k),(k)=φT(-,u,),
=[a1,...,an,b1,...,bm,c1,...cr]T;
误差结构图,
e(k)=Cw(k)/A
u(k) 过程 y(k)
B/A
(z-1) -1(z-1) (k)
( z-1)-1
分析：
(k)= u+e-[u+(-1)(k)]=( -)u+[e- (-1)(k)]
若e=w，且u与w不相关，则→w意味着：→,→(-1)+1
若e=w，且u与w不相关，则→w意味着：→，→
针对Ay=Bu+w/D的方程误差法计算式：φ(-y,u,-ξ)=[-y(k-1),...,-y(k-n),u(k-1),...,u(k-m),-ξ(k-1),...,-ξ(k-r)]T;
ξ(k)=y(k)-φT(-y,u,0),(k)=y(k)-(k),(k)=φT(-y,u,-ξ),
=[a1,...,an,b1,...,bm,d1,...dr]T;
误差结构图,
e(k)=w(k)/AD
u(k) 过程 y(k)
B/A
(z-1) (z-1)
( z-1) (k)
分析：
(k)= [y-u]= [(-)u+e]
若e=w，且u与w不相关，则→w意味着：→，→DA
但不能保证→A，→B，→D
针对y=Bu/A+Cw的输出误差法计算式：φ(-,u,)=[-(k-1),...,-(k-n),u(k-1),...,u(k-m),(k-1),...,(k-r)]T;
(k)=y(k)-(k),(k)=φT(-,u,),
=[a1,...,an,b1,...,bm,c1,...cr]T;
误差结构图,
e(k)=Cw(k)
u(k) 过程 y(k)
B/A
(z-1) (k)
( z-1)-1
分析：
(k)= u+e- [u+(-1)(k)]=( -)u+[e- (-1)(k)]
若e=Cw，且u与w不相关，则→w意味着：→,→C
针对y=Bu/A+w/D的输出误差法计算式：φ(-y,u,-ξ)=[-y(k-1),...,-y(k-n),u(k-1),...,u(k-m),-ξ(k-1),...,-ξ(k-r)]T;
ξ(k)=y(k)-φT(-y,u,0),(k)=y(k)-(k),(k)=φT(-y,u,-ξ),
=[a1,...,an,b1,...,bm,d1,...dr]T;
误差结构图,
e(k)=w(k)/D
u(k) 过程 y(k)
B/A
(z-1) ( z-1) (k)
分析：
(k)= [y-u]= [(-)u+e]
若e=w，且u与w不相关，则→w意味着：→，→D
递推算法的实现
UD分解算式P(t)=，t=t-1+P(t)对舍入误差敏感，易导致病态。令P(t)=U(t)D(t)UT(t),P(t-1)=U(t-1)D(t-1)UT(t-1)。定义：d=dim(θ)
f=UT(t-1)ψ，g= D(t-1)f，β=λ(t)I+ψTP(t-1)ψΛ=λ(t)I+fTgΛ=λI+[]Λ
P(t)=[U(t-1)D(t-1)UT(t-1)- U(t-1)D(t-1)UT(t-1)ψΛβ-1ψTU(t-1)D(t-1)UT(t-1)]
=U(t-1)[D(t-1)-gΛβ-1gT]UT(t-1)
D(t-1)-gΛβ-1gT===+
=[ D(t-1)-gΛβ-1gT]-=-+DdededT-gΛβ-1gT
定义：Md=-DdededT+gΛβ-1gT，注意到Md的最后一行和最后一列全为零，因而有：
=Dd-VTddΛβ-1Vdd，=1，=- VTdiΛβ-1Vdd，其中：
g=Vd==。定义βk=λI+[]Λ,βd=λI+[]Λ=β，注意到gi=Difi=VTid,所以：βk+1=βk+gk+1Tgk+1Λ=βk+ Vdk+1VTdk+1Λ，或者写成：
βk-1=βk+1-1+βk-1Vdk+1VTdk+1Λβk+1-1，各项右乘Vdk+1得：
βk-1Vdk+1=βk+1-1Vdk+1+βk-1VdiVTdiΛβk+1-1Vdk+1=βk+1-1Vdk+1+βk-1Vdk+1[Dk+1-]，于是有：βk+1-1=βk-1。定义：
Vd-1=，我们有：=-Vd-1Λβ-1Vdd+，因而,
Md=Vd-1Λβd-1VTd-1+Vd-1Λβd-1VddVTddβd-1ΛVTd-1=Vd-1Λ[βd-1+βd-1-1VddVTddβd-1Λ]VTd-1
= Vd-1Λβd-1-1VTd-1。至此，我们把化成与相同的形式，即=-Vd-1Λβd-1-1VTd-1，
=D(t-1)-gΛβ-1gT=-VdΛβd-1VTd，仿此，可求出所有的和：
=Diβi-1βi-1,=-giΛβj-1-1fTj。以下是UD分解的计算流程：
初始化U(0)、D(0)，t=0，
1 t=t+1；计算f=UT(t-1)ψ，g= D(t-1)f，β0=λ(t)I；
2 对于j=1,...,d=dim(θ)，计算下列3-5步；
3 计算：βj=βj-1+fTjgjΛ；
Djj(t)= Djj(t-1)βj-1βj-1；
Vj=gj；
Qj=-Λβj-1-1fTj；
4 对于i=1,...,j-1，计算第5步(j=1时跳过第5步)；
5 计算：Uij(t)=Uij(t-1)+ViQj；
Vi=Vi+Uij(t-1)Vj；
6 计算：P(t)ψ=βd-1。
稳定性判断递推算法失败的一个原因是模型在修正过程中变为不稳定模型。解决的办法是在每次修正后判断模型的稳定性，若模型稳定，则接受修正，否则，将模型的修正量减半。
Jurry-AstrOm稳定判据
设线性定常离散系统的特征方程为
A(z)=a0zn+a1zn-1+......+an-1z+an=0
按下表建立Jurry表：
a0 a1 a2,....,an-1 an
-) an an-1 an-2,....,a1 a0 x an/a0
b0 b1 b2,.....,bn-1
-) bn-1 bn-2 bn-3,.......,b0 x bn-1/b0
,...........................................................
,..........................................................,
r0 r1
-) r1 r0 x r1/r0
q0
特征方程稳定的充分必要条件是：a0 >0,b0 >0,......,r0 >0,q0 >0
以上判据很容易编制成程序。
递推算法流程
1 指定模型结构，d=dim(θ);
2 初始化P(0)、θ(0)、ψ(0)、φ(0),L(0)= P(0)ψ(0),k=1，指定Λ及λ(0)
3 计算 e=y-φTθk-1,置参数修改步长Amy=1
4 预计算新参数θk-1=θk-1+LΛe×Amy
5 测试分母多项式是否稳定
6 若不稳定，令Amy=Amy/2，返回4
7 计算残差ε= y-φTθk
8 修正φ、ψ（包含需要计算的输出估计、辅助变量、滤波等）
9 计算增益向量L，修正P，λ=aλ+b
10 若数据用完，则结束；否则k=k+1,返回3。
其它
P(0)、θ(0)的初始值：P(0) =107I，θ(0)=0，P(0)充分大是为了保证有足够大的参数搜索区域，但θ(0)的选择还是要尽可能地接近真实参数。
数据复用：为了解决递推算法中信息利用不充分的问题，可采取数据复用的办法，即在数据用完后，以所得的参数估计作为下一轮递推估计的初始值，再次执行递推程序。新一轮的递推程序中，P(0)可选为P(0)=10 P(N)。
结果求期望：θ=
遗忘因子λ(t)：遗忘因子的作用有两个，一是对定常系统的辨识，在递推估计的初始阶段，λ(t)<1可加快参数向真值的运动，但在递推估计的结束阶段应保持λ(t)=1，否则会影响参数的收敛精度。二是对时变系统，遗忘因子可使参数估计跟随实际参数的变化，但要注意：1 参数变化频率应远小于系统响应频率；2 λ(t)的大小应适应参数变化频率的要求，λ(t)小，参数变化剧烈。
λ(t)的递推式通常为：λ(t)=aλ(t-1)+b，其中a、b应满足b<=1-a。λ(∞)=b/(1-a)。
参数a、b及λ(0)的选取决定了λ(t)的运动轨迹。在递推辨识中，通常选取λ(0)>0.8；a的大小控制着λ(t)的收敛速率，对于初始值和终值均相同的两条轨迹，a大则收敛慢，a小则收敛快。
第六讲 结构辨识
Hanker矩阵与积矩矩阵
Hanker矩阵：设g(1),g(2),.....g(L)为脉冲响应序列，定义Hanker矩阵为：
H(j,k)=，j>n0。Rank(H(j,k))= n0，任意k。对于有噪声的情况，可令Rg(k)为脉冲响应序列的自相关函数，用Rg(k)替代g(k)构造H(j,k)即可。
积矩矩阵：考虑差分方程：A()y(k)=B()u(k)+ w(k)，定义：
H*n=，L>n,Rank(H*n)=min(n+n0,2n)。因为U部分总是满秩。定义：Hn=H*Tn H*n，DR(n)=det(Hn)/ det(Hn+1)，若DR(n)较DR(n-1)显著增加，则n比n-1更接近n0，应选n。
残差方差
考虑伪线性回归模型：y(k)=φT(k)θ+w(k)，n0=dim(θ)。模型的输出残差为：
=y(k)-φT(k)=+w(k)，其中，=φT(k)θ-φT(k)，n=dim()。残差方差为：
V=T=(T++)，因为w(k)为白噪声，所以当L→∞时：
V=T+ **
其中T，因此，V=，n≥n0。
在实际计算中，V(n)的变化不可能在n0处有明显的拐点。
F检验法：
t=~F(n2-n1,L-n2)，若设定风险水平为α=0.05，则Fα(2,100)=3.09。据此，我们可通过假设检验来判断V(n)的变化：当n2-n1=2，L-n2>100时，若t>3.09，则V(n2)是显著大于V(n1)的，n2比n1更接近n0；若t<3.09，则V(n2)不显著大于V(n1)，因此，可取n0=n1。
残差自相关性模型的残差带有两方面的信息，一是参数误差信息，二是附加噪声信息。**式也告诉我们，当模型阶次大于等于实际系统阶次时，残差中主要包含附加噪声信息。因此，通过判断残差与附加噪声的统计性质是否一致，也可得出模型是否具有合适的阶次。
当附加噪声为白噪声时，可用残差的自相关函数来判断它的白色性。常用方法为：
设R(0),R(1),...为残差的自相关函数，ρ(k)=R(k)/R(0)，若：E[ε]=≈0，且|ρ(k)|≤，(k可取1~20)，则可认为残差为白噪声序列。
最终预报误差在辨识过程中，我们只能得到有限次(L)观测数据，由此得到的模型L对未来输出的估计误差如何是我们所关心的。这是用最终预报误差最小来确定模型阶次的出发点。记为根据模型L所得的输出预报。残差方差的有限数据估计为：=，最终预报误差可表示为：==E{[z(k)- ]2}，而：
E{[z(k)- ]2}=E{[hT(k)+w(k)]2}= E{T h(k)hT(k) }+=+
=Tr[(HHT)T] +
因为=（(HHT)-1Hw，所以，=dim(L)+。又因为=，
最后有：=。
AIC准则在极大似然原理中，我们通过极大化p(ZL|θ)来求取参数的极大似然估计，其中ZL是已经得到的有限次观测数据。设系统的真实参数为θ0，由极大似然原理，我们有以下关系式：
，+，
→，
2=2+2{+}
2-2=～χ2(n)，n=dim()，两边求期望，则有：
E{2}- E{2}=n。
如果考虑随机变量z，并定义极大似然估计的信息测度：I(θ0,)=E{}≥0，那么，使I(θ0,)最小的模型结构就是我们所要求的。我们有：
I(θ0,)=E{}=S(θ0)-S()=S(θ0)-E{logp(z|)}，将S()在θ0附近展开：
S()=S(θ0)+E{ +T+............}
因为S(θ)的最大值出现在θ0，所于=0，又因为=-M(Fisher信息矩阵)，因而：I(θ0,)= E{TM}，其中TM～χ2(n)，n=dim()，或者：E{2}- E{2}=n，注意到E{2}= E{2}，可得：
-E{2}=- E{2}+2n=AIC。现在，使I(θ0,)最小转变为使AIC最小。
在正态分布的假定下，-=+Llnσ+（L/2）ln2π，其中
→L/2，所以，AIC= Llnσ2+2n
逐步回归法对于一般的回归方程：y=，可以认为y是由xi构成的，但不同的xi对y的贡献是不同的。方程右边不包含变量z，是因为z对y的贡献为零。这一简单思想可用来确定模型结构。现在的问题是如何定义“贡献”及如何计算“贡献”。
定义：S总=，总体偏差平方和
S回=，回归平方和，
S剩=，剩余平方和三者之间的关系为：S总= S回+ S剩+
假定1 及均为零；2 =0，=0；
则S总= S回+ S剩，并且S剩== S总-，即S回=。
类似地，当我们在a1，a2，....，an中选取p个元来构造回归方程时，定义：=。
下面我们观察一下选元过程。定义：
XA==0，令C=XTX，则CA=0。
C==C0=
第一次消元，p=1：假定对第一列消元，则得到的第一个回归系数为：==，消元以后，C0转变为C1，其各元值为：c1(1,j)= c0(1,j)/c0(1,1)；
c1(i,j)=c0(i,j)- c0(i,1) c0(1,j)/c0(1,1)= c0(i,j)- c0(i,1) c1(1,j)。
c1(i,n+1)==Ei- c0(i,1) / c0(1,1)E1= Ei- c1(1,j) E1，i=2,...,n。
c1(n+1,j) ==Ej- c0(1,j)/c0(1,1) E1= Ej- c1(1,j) E1，j=2,...,n。
c1(n+1,n+1)= S总-E1= S1剩。
第二次消元，p=2：假定对第二列消元，则得到的第二个回归系数为：
==-E1=。此时，第一个回归系数改变为：=- c1(1,2)。消元以后，C1转变为C2，其各元值为：
c2(2,j)= c1(2,j)/c1(2,2)；j=3,...n；
c2(i,j)=c1(i,j)- c1(i,2) c1(2,j)/c1(2,2)。定义h22=1/ c1(2,2)，
c2(i,n+1)==- c1(i,2) / c1(2,2) =- c2(2,i)，i=3,...,n。
C2(n+1,n+1) = S1剩- c1(n+1,2) = S总-E1- [E2- c1(1,2) E1]= S总-E1- E2= S2剩。
一般地，当增加第p个参数时，原来已得到的参数的修正公式为：
=- cp-1(j,i)，也即：=- cp-1(j,i)。新补参数的计算公式为：
==-
定义偏回归平方和Qi：当增补参数和相应的xi到回归方程时，定义Qi=-为xi的偏回归平方和，以此来衡量xi对回归方程的贡献。
Qi=-=-
=--=Ei-
=cp-1(i,i)[ -]=()2 cp-1(i,i)=
= Cp(n+1,n+1)= Cp-1(n+1,n+1)- cp-1(n+1,i)= Cp-1(n+1,n+1)- Qi
注意Qi的计算可在第p次消元之前进行，这就允许我们在候选元中选择对回归方程的贡献最大的元。
为避免过估计，应当进行显著性检验。已经有人证明了：
统计量F= Qi/i(N-p-1)～Fα(1,N-p-1)，当置信度设为95%时，若N-p-1>150，且F>3.9时，可认为xi对回归方程的贡献是显著的。
将逐步回归法应用于结构辨识时，应注意：
1 设定的初始结构要大于系统的时延与阶次的和；
2 为了得到符合实际的模型，应保证各多项式的连续性（即若多项式的阶次为n，则指数小于n的参量的系数不应为零）。为此，对候补元的确定规定如下：
开始时，p=[-y(k-1),全部的uj]；在一次消元后，若已选中的元有[-y(k-1),...,-y(k-i),u(k-j),u(k-j-1),...,u(k-r)]，则后选元为p=[-y(k-i-1),u(k-j+1),u(k-r-1)]。
逐步回归法算法流程：
1 读入数据；
2 建立C0，令p=0；
3 p=p+1，确定候补元集合；
4 对候补元集合中的每个元计算Q值；
5 挑出Q值最大的元；
6 计算及F= Qi/i(N-p-1)，如果F>3.9，转7；否则转8；
7 进行第p次消元，转3
8 显示模型结构及参数，辨识结束。
系统辨识大作业
已知一系统为两输入单输出系统，观测数据受有色噪声污染，噪信比为N/S=0.1。系统经2000次采样，存放于文件T3T.TXT中。系统输入u1为7级M序列，u2为u1的63步移位序列。模型类可选为：A(q-1)y(k)=B1(q)u1(k)+ B2(q)u2(k)+w(k)/C(q-1)。要求编制程序，辨识出该模型的结构及参数。
作业文档要求：
描述问题；
选择辨识方法并简单说明所选方法中的结构辨识原理和参数估计原理；
程序流程图及程序清单；
说明程序中用到的一些技术，如数据标准化、UD分解、稳定性判断等；
结构搜索路线及各结构下的参数、残差；
给出最终结果：A(q-1)=
B1(q)=
B2(q)=
C(q-1)=
并给出选择此最终结果的理由；
用你的辨识结果来预报系统输出误差e(k)=y(k)-(k)，并画出e(101)-e(400)的曲线图。
