《系统辨识与建模》.doc

系统辨识与建模课时：40参考书：1 徐南荣、宋文忠、夏安邦 系统辨识，东南大学出版社，19912 方崇智、萧德云 过程辨识，清华大学出版社，19883 Ljung L. And Soderstrom T. Theory and Practics of Recursive Identification, MIT.Press, Cambridge4 美P.艾克霍夫 系统辨识-参数和状态估计，潘科炎等译，科学出版社，19775 美夏天长 系统辨识最小二乘，熊光楞、李芳芸译，清华大学出版社，1983期刊：1 Automatica2 Proc. IFAC Identification and System Parameter Estimation3 IEEE Trans. On Automatic Control4 自动化学报5 控制理论与应用学习要求1 培养独立学习一门新课程的能力，为今后学习和研究打下基础（要求大家尽量少依赖听课，多自学）。2 掌握基本的辨识理论和辨识技术3 能独立设计辨识实验，并编程计算4 学习一些现代建模技术考核办法给出一个数据文件，通过编程对其进行辨识，并写出报告第一讲 概论实体与模型实体：客观存在的事物及其运动状态，有时也称之为“系统”模型：实体的一种简化描述。模型保持实体的一部分特征，而将其它特征忽略或者变化。不同的简化方法得到不同的模型。模型分类：直觉模型：地图、建筑模型、照片、软件演示文档等物理模型：风洞、水力学模型、传热学模型、电力系统动态模拟模型等。（缩小的复制品）数学模型：描述实体中一些关系和特征的数据模型。例如：投入/产出模型、热源与室温的关系模型等。数学模型数学模型还可分为：图表模型：如阶跃响应、脉冲响应、频率响应、温度与热电偶输出关系表解析模型：代数方程、微分方程、差分方程、状态方程程序模型：神经网络仿真程序语言模型：模糊关系模型获得数学模型的方法有：经验总结法：模糊关系模型、静态线性关系模型机理分析法：解析模型实验法： 图表模型数据拟合法：解析模型、程序模型用数据拟合法获得解析模型的过程即为系统辨识。系统辨识L. A. Zadeh1962 ：辨识就是在输入和输出数据的基础上，从一组给定的模型类中，确定一个与所观测系统等价的模型。这一定义给出了系统辨识的三要素：数据、模型类和准则。 数据：由观测实体而得。不唯一，受观测时间、观测目的、观测手段等影响。 模型类：规定了模型的形式。不唯一，受辨识目的、辨识方法等影响。 准则：规定了模型与实体等价的评判标准。不唯一，受辨识目的、辨识方法等影响。系统辨识的三要素是评判数据拟合方法优劣的必要条件，只有在相同的三要素下，才可区分数据拟合方法的优劣；而在不同的三要素下，这种结论也会改变。（图1）辨识目的明确模型应用的最终目的是很重要的，因为它将决定如何观测数据、如何选择三要素以及采用什么数据拟合方法等。而最根本的是它将影响辨识结果。辨识目的主要取决于模型的应用。辨识模型应用有以下几个方面： 辨识目的 选择 准则 模型类 辨识实体 数据观测 数据 数据拟合 模型 图 1 系统辨识三要素1 验证理论模型；要求：零极点、结构（阶次及时延）、参数都准确；模型类同理论模型。2 设计常规控制器；要求：动态响应特性、零极点、时延准确；便于分析的模型类。3 设计数字控制器；要求：动态响应特性、时延准确；便于计算机运算的模型类。4 设计仿真/训练系统；要求：动态响应特性准确；便于模拟实现的模型类。5 预报预测；要求：动态响应特性、时延准确；便于计算机运算的模型类。 6 监视过程参数，实现故障诊断；要求：参数准确；能直观体现被监视过程参数的模型类。10 系统的定量与定性分析；要求：静态关系准确；模型简单，便于人脑判断。辨识的一般步骤我们将结合一个实际例子来说明辨识的一般步骤。例：上图为长网造纸的流程简图。对于造纸企业来说，质量控制就是要控制好成品纸的定量与水份。而纸的定量与水份与纸浆浓度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 02 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 ，频谱宽度无限。，=00，03 ，其中，为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个游程，游程长度为1n，但出现的概率是随机的。长度为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充分大时，均值约等于02 在一个周期内，自相关函数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=(315)/T）则以1.5*T3* fmax为原则，选择n 和Ts。 实际应用过程：1 确认系统过渡过程T;2 选择最高工作频率fmax=(315)/T；（设为fmax=k/T）3 根据1/Ts3* fmax，得出Ts=T/(3k)；4 根据1.5*T1.5*(输出对各输入的过渡过程时间）。白噪声序列1/Ts3* fmax数据采集设计采样间隔满足香农定理，设采样间隔为T0，1/T02* fmax。对于有扰动信号的系统，T00，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)=Zu(t)=Zu(kT) y(z)=Zy(t)=Zy(kT) u(kT) x(t) y*(t) L(s) G(s) u(z)=Zu(kT) X(s) y*(z)=Zy*(t) G(z)= y*(z)/ u(z)=ZG(s)L(s) G(s)=Z-1G(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)=ZG(s)，G(s)=Z-1G(z-1)L2：（零阶保持器）x(t)=u(kT)，kTt(k+1) T，L(s)=(1-e-Ts)/s当u(t)=1(t)时，y*(t)= y(t)，故称为阶跃响应不变等效。G(z-1)=(1-z-1)ZG(s)/s，G(s)= Z-1G(z-1)/(1-e-Ts)/sL3：x(t)=u(kT)*(k+1)T-t)/T+u(k+1)T)*(t-kT)/T，kTt(k+1)T，L(s)= (1/T+s)(1-e-Ts)/s2 斜坡响应不变等效。2 Z变换与逆 z=eTs， s=ln(z)/T：（可用查表法进行） G(z)=ZG(s)=G(s)| s=ln(z)/T，G(s)=Z-1G(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-anb0C=1 0 .0， D=b0连续状态方程与连续传递函数的关系与此类同。连续状态方程与离散状态方程:=AX+BU Y=CX+DU X(k+1)=Fx(k)+Gu(k) y(k)=Hx(k)+Eu(k)F=eATI+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, , bmT =-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)+ cy(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)-cy(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, , bmT (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)-1T YN-(T)-1T WNLS=(T)-1TYN, 条件1：E(T WN)=0 （w(k)白色） 条件2: T可逆以上结果等同于求使：J=(y(k)-T(k)2最小的，因此称为最小二乘算法。最小二乘算法的MATLAB程序：clear all; %最小二乘法for MISOload t2; %读入数据矩阵 tt=t2;clear t2;ll=size(tt); %得到矩阵维数200x2r=ll(2)-1;m(1)=input(order of A(z);tao(1)=0;for i=1:rim(i+1)=input(order of B(z);tao(i+1)=input(timedelay of B(z);endif m(1)max(m),m(1)=max(m);endn=m(1)+max(tao);lll=ll(1)-n;in1=r+1;kn=0;for k=1:in1for i=1:lllfor 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);endend;end;kn=kn+m(k);end;for i=1:lllyy(i)=tt(i+n,1);end;fa=ff*ff;fay=ff*yy;zz=inv(fa);zta=zz*faym，tao 加权最小二乘：若对各次观测数据加不同的权，即求使J=k(y(k)-T(k)2最小的,则得到参数的加权最小二乘估计：LS=(T)-1TYN参数估计的一般性质1 无偏性 如果E()=0，（=-） 或E()=，则称估计为无偏的。2 有效性如果无偏估计满足Cov()=M-1，则称估计为有效的。其中： 称为Fisher信息矩阵，其逆M-1称为Crammer-Rao 下界。在一般情况下，有Crammer-Rao不等式：Cov()=E(-)(-)TM-1对于z=H+w，若噪声w为零均值、协方差阵w=E(WWT)的正态分布噪声，即：WN(0，w)，则输出信号ZN(E(H)，w)，于是， ，因此： ，于是，M=E(HTw-1H)，其Crammer-Rao不等式为：Cov()E(HTw-1H)-1。有效估计也称为最小方差估计、马尔可夫估计。3 一致性若估计为渐近无偏的（E()=0），且Var()=0，则称为一致估计。Var()=E(-)T(-)最小二乘、加权最小二乘估计性质估计性质最小二乘加权最小二乘无偏性= (T)-1T WN，当:1 w(k)为零均值独立随机过程（白噪声）时；或者2 E(T WN)=0时，有：E()=0。无偏条件较严格。= (T)-1TWN无偏条件同最小二乘估计。有效性Cov()=E(T)-1TWNWTN(T)-1=(T)-1Tw(T)-1M-1,若W是服从正态分布的独立随机过程，则w=I，因而Cov()=I(T)-1= M-1。Cov()=(T)-1Tw(T)-1，当=w-1时，Cov()=(Tw-1)-1= M-1，此时不要求W是服从正态分布的独立随机过程一致性若W是服从正态分布的独立随机过程，则Cov()=0,其中收敛于一个有界非奇异矩阵。而Var()只是Cov()的对角线上元的和，因而也有：Var()=0。同最小二乘估计噪声方差估计e= YN-=+WN-=+WN-(T)-1TYN =WN-(T)-1T WN =(I-(T)-1T) WN=B WN， 因为BT=B ,B2=B所以，E(eTe)=E(WN TB WN)= E(WN T WN- WN T(T)-1T WN) =N- (Tr(T)-1T)= (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;3 用最小二乘估计A()(k)= B()(k)+ w(k)中的参数；4 用估计模型()、()以及各时刻的观测数据，计算出残差(k): (k)= ()y(k)- ()u(k) 5 计算i=2(k)及=i-1-i，如果小于一定数，则结束辨识。否则转下一步。6 对于噪声模型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)-1TY-(T)-1T=D-1TMYM=I-(T)-1T D=TM须要注意，的求取仍然是一个迭代过程。偏倚校正算法的计算步骤为：1 令C()=1，用最小二乘法求=LS=(T)-1TY，并保留、=(T)-1T以及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)-1TY ，然后依据()z(k)= ()u(k)计算出输出估计z(k)，再用z(k)、u(k)