《控制系统辨识算法与仿真》全套教学课件

系统辨识与建模全套可编辑PPT课件

线性系统理论是基础，解决系统的模型描述和基础知识，即线性系统一般可描述为：

最优控制解决在某一性能指标约束下，如何解算最优输入u(t)；

最优估计主要解决状态变量X的估计和预测。

系统辨识目的系统辨识目的上述问题解决的先决条件：

模型中的A、B、C、D已知。

亦即系统的结构和参数已知，也就是要知道系统的传递函数、或是脉冲传递函数、或是差分方程、或是系统的频率特性。

那么，如何获取系统的结构和参数？

系统辨识目的：

如何获取系统的模型及其参数？学习要求

毕业要求指标点课程教学目标2.问题分析2.2能够应用数学、自然科学和工程科学知识，对自动化领域复杂工程问题中的控制对象及控制过程建立数学模型。目标1：从系统输入输出数据入手，运用适当的系统辨识方法，如阶跃响应法、脉冲响应法、相关分析法、最小二乘法以及其改进算法、极大似然法等，建立起控制对象或控制过程的数学模型。4.研究/现代工具的使用4.1能够对自动化工程问题存在的各类物理现象进行观测和分析，明确其中的关联因素和本质特征，并对自动化系统或控制对象进行合理建模或描述。目标2：能根据实际问题选择合适的辨识方法进行系统模型的建立和仿真分析。掌握建模过程中数据准备和处理、辨识方法的使用以及利用一定的实验手段对建模过程和结果进行验证；初步掌握用Matlab进行系统建模与仿真分析设计的方法。学习要求

培养独立学习一门新课程的能力，为今后学习和研究打下基础（要求大家尽量少依赖听课，多自学）。掌握基本的辨识理论和辨识技术能独立设计辨识实验，并编程计算学习一些现代建模技术

考核方式及成绩评定（1）考核方式本课程为考试课程，期末考试为开卷笔试。（2）成绩评定总评成绩=笔试60%+过程考核40%。课程内容第一章绪论第二章辨识相关基础知识第三章线性系统经典辨识方法第四章最小二乘参数辨识与极大似然辨识第五章模型结构辨识第六章基于遗传算法的参数辨识及其应用第七章基于差分进化算法的参数辨识第一讲

绪论

系统模型建立系统数学模型的基本方法系统辨识的基本概念系统辨识步骤辨识的发展和应用1.1系统数学模型及建模方法系统：按某种相互依赖关系联系在一起的客体的集合。

将研究的对象看成是一个系统，从整体行为上对系统进行研究是现代科学技术中的一种十分重要的方法论。这里研究的“对象”是抽象的，重要的是其输入、输出关系，只要处理“输入输出的因果关系”就可以将许多不同的问题进行统一的处理。模型：把关于实际系统的本质的部分信息简缩成有用的描述形式。

模型可以用来揭示系统的运动规律，是系统的一种客观写照或缩影，是分析、预测和控制系统行为特征的有力工具。对实际系统而言，模型一般不可能考虑所有因素。如果要求模型越精确，模型就会变得越复杂；相反，如果适当降低模型的精度要求，只考虑主要因素而忽略次要因素，模型就可以简单一些。在建立实际系统的模型时，存在着精确性和复杂性的矛盾，找出两者的折中解决办法往往是建立实际系统模型的关键。1.1系统数学模型及建模方法实体与模型

实体：客观存在的事物及其运动状态，有时也称之为“系统”模型：实体的一种简化描述。模型保持实体的一部分特征，而将其它特征忽略或者变化。不同的简化方法得到不同的模型。建模：凡是用一个数学模型去表示系统的某种因果关系，都属于建模的范畴，建模涵盖了将近一个世纪系统模型化的过程。

模型的表现形式：直觉模型。它指系统的特性以非解析形式直接存储在人脑中，靠人的直觉控制系统的变化。例如：司机对汽车的驾驶，指挥员对战斗的指挥，依靠的就是这类直觉模型。物理模型。它是根据相似原理，把实际系统加以缩小的复制品，或是实际系统的一种物理模拟。例如，建筑模型，风洞、水洞模型，传热学模型，电力系统动态模拟等，均是物理模型。图表模型。它以图形或表格的形式来表现系统的特性。如阶跃响应、脉冲响应和频率特性等。图表模型也称为非参数模型。数学模型。它用数学结构的形式来反映实际系统的行为特性。当模型的结构阶次和参数确定后，数学模型也就确定了。数学模型还可分为：参数模型：用数学表达式描述的模型。如：代数方程、微分方程、差分方程、状态方程非参数模型：如阶跃响应、脉冲响应、频率响应、温度与热电偶输出关系表关系：非参数模型可通过实验获得；可转化为参数模型。数学模型的分类数学模型的分类（2）非线性模型。非线性模型用来描述输入输出变量之间不满足线性关系的非线性系统，一般不满足叠加原理。在讨论线性和非线性问题时，要注意2点区别：系统线性是指模型的输出关于输入变量是线性的，而不是模型的输出关于参数是线性的，如果模型的输出关于参数是线性的，则称之为参数空间线性。如果模型经过适当的数学变换，将本来是非线性的模型转换成线性模型，则原模型称作本质线性，否则原模型称作非本质线性。163数学模型的分类（3）动态模型。动态模型用来描述系统处于过渡过程时各状态变量之间的关系，其一般为时间的函数。（4）静态模型。静态模型用来描述系统处于稳态时（各状态变量的各阶导数均为零）各状态变量之间的关系，一般不是时间的函数。（5）确定性模型。由确定性模型所描述的系统，当状态确定后，其输出响应是唯一确定的。（6）随机性模型。由随机性模型所描述的系统，当状态确定后，其输出响应仍然是不确定的。（7）宏观模型。宏观模型用来研究事物的宏观线性，一般用联立方程或积分方程描述。（8）微观模型。微观模型用来研究事物内部微小单元的运动规律，一般用微分方程或差分方程描述。建立模型的目的（1）深入研究系统的一种手段感性认识—提炼—理性，探讨发展规律（2）分析、设计、仿真的基础对系统进行设计、分析，需要数值仿真（3）预测利用模型预测系统工况和性能建立系统数学模型的基本方法1机理建模(理论建模，白箱建模)

机理清楚不适合复杂系统“白箱”建模机理分析(化学,物理,物料、能量平衡,传热传质)机理建模实例—电动机的数学模型他励直流电动机在电流连续时，电压U和转速n之间的传递函数：其中：机电时间常数电磁时间常数建立系统数学模型的基本方法e拟合,统计分析外特性等价适合复杂系统建模机理不清“黑箱”建模2系统辨识（试验建模、黑箱建模、数据建模）建立系统数学模型的基本方法2系统辨识（试验建模、黑箱建模、数据建模）实际系统---输入/输出数据---数学模型优点：无需先验信息、专业知识要求少，建模时间短，提高动态特性模型适用范围：较复杂系统建立系统数学模型的基本方法3机理分析和系统辨识相结合（灰箱建模）-机理已知的部分采用机理建模，机理未知的部分采用辨识建模-利用机理建模确定模型的结构，利用辨识建模确定模型的参数。1.2辨识的定义、内容和步骤L.A.Zadeh[1962]：辨识就是在输入和输出数据的基础上，从一组给定的模型类中，确定一个与所观测系统等价的模型。这一定义给出了系统辨识的三要素：数据、模型类和准则。系统辨识的三要素数据：由观测实体而得。不唯一，受观测时间、观测目的、观测手段等影响。模型类：规定了模型的形式。不唯一，受辨识目的、辨识方法等影响。准则：规定了模型与实体等价的评判标准。不唯一，受辨识目的、辨识方法等影响。

系统辨识的三要素是评判数据拟合方法优劣的必要条件，只有在相同的三要素下，才可区分数据拟合方法的优劣；而在不同的三要素下，这种结论也会改变。（图1）

数据拟合选择辨识目的准则模型类辨识实体数据观测数据模型系统辨识的三要素辨识算法的基本原理要素1要素1要素2要素3批处理递推系统辨识的内容和步骤第一步：验前知识和辨识目的

明确模型应用的最终目的是很重要的，因为它将决定如何观测数据、如何选择三要素以及采用什么数据拟合方法等。而最根本的是它将影响辨识结果。辨识目的主要取决于模型的应用。将决定：模型类型、精度要求、辨识要求、辨识方法

模型的要求用于设计或预报（理论模型验证、过程参数检测、故障检测）：模型准确度要求高。用于控制：模型准确度要求低，确定性系统的自适应控制要求又高一些。用经典理论控制：需要非参数模型（脉冲响应、阶跃响应、频率特性）用现代控制理论（二次型性能指标的最优控制、Kalman滤波）：需要状态空间模型。第二步：实验设计（1）选择变量输入变量应该能够设置输出变量应该能够测量与我们感兴趣的现象有关输入信号的基本要求：输入信号必须包含有足够丰富的频率分类。通常选阶跃信号，脉冲信号，伪随机信号，自然干扰等（2）确定实验期限长短和采样间隔实验时间:受干扰、漂移、经济性等限制，一般取系统主时间常数的10倍。(3)确定开环或闭环辨识一般做开环辨识本身固有反馈系统，必须做闭环辨识注意：可辨识性可辨识性：所用模型是否等价于真实系统，是否能得到参数的唯一解。(4)输入/输出数据的产生、检测、存储选择产生信号的设备（信号发生器、计算机）输入、输出数据的检测、存储。(5)离线或在线辨识（一）离线辨识(5)离线或在线辨识（一）离线辨识

要求：系统的模型结构和阶次预先确定(1)过程：系统模型及阶次n选定后，记录下系统全部的I/O数据，然后再用参数估计方法，辨识系统的模型参数。(2)特点：需存储数据量大，计算量大，辨识精度较高。事后数据处理方法，不能用于实时控制系统。(5)离线或在线辨识（二）在线辨识(5)离线或在线辨识（二）在线辨识要求：模型结构和阶次事先确定(1)过程：系统模型及阶次n选定后，先获取一小部分数据，估计系统模型参数，再获取新的I/O数据，采用递推修正算法获得新的参数估计值，重复上述过程，直至系统运行停止。(2)特点：数据量小，计算量小，辨识精度稍低。是一种在线数据处理方法，用于实时控制系统。第三步：确定模型结构模型形式：静态-----动态线性------非线性时变-----非时变确定性----随机集中参数---分布参数时域------频域参数------非参数模型阶次的选择纯滞后时间的估计第四步：参数估计模型的未知部分是以未知参数的形式出现参数估计的方法是本课程的重点第五步：模型验证实际测量输出和模型的计算输出比较模型参数应当保证两个输出之间在选定意义上的接近若不一致，经修改稿模型结构的假设，修改实验设计，重复实验1.3系统辨识技术发展和应用技术发展：20世纪30-50年代，经典控制理论

建模占主要地位，针对单输入单输出系统20世纪50年代后，现代控制理论

对建模提出新要求，引入了系统辨识技术，使用多种对象：单输入/输出、多输入/多输出、连续/离散、定常/慢时变1.3系统辨识技术发展和应用系统辨识应用：航天、航空、航海、交通运输、化工、冶金、核能、电力、机械、楼宇自动化等工程技术部门，社会境界系统等方面也有使用。1.3系统辨识技术发展和应用系统辨识研究动向

系统辨识和反馈控制相结合的自校正调节器闭环系统的辨识问题多变量系统的结构辨识时变系统的参数辨识辨识算法的收敛性分析小样本下的建模问题如何把系统辨识应用于分布参数系统，模糊系统，机器人工程，智能系统等1.3辨识中常用的误差准则辨识时所选用的误差准则是辨识的3个要素之一，是用来衡量模型接近实际系统的标准。误差准则也称为等价准则、损失函数、准则函数、误差准则函数等。它通常被表示为误差的泛函数，记为

1.3辨识中常用的误差准则输出误差准则输入误差准则1.3辨识中常用的误差准则广义误差准则课堂作业1建立系统数学模型的基本方法有哪些？

2系统辨识的定义和三要素分别是什么？

3写出系统辨识的步骤系统辨识与建模

第一章绪论回顾建立系统数学模型的基本方法1机理建模(理论建模，白箱建模)2系统辨识（试验建模、黑箱建模、数据建模）3机理分析和系统辨识相结合（灰箱建模）第一章绪论回顾系统辨识的三要素数据模型类准则系统辨识步骤第二章数学预备知识2.3控制系统的时间响应分析（选讲）

2.2随机过程的定义及其数字特征2.1随机变量的定义及其数字特征54/175（一）随机变量的定义

设(S,F,P)是一概率空间。对于sS，X(s)是一个取实数值的单值函数。若对于任意实数x，{s:X(s)<x}是一随机事件，亦即{s|X(s)<x}F，则称X(s)为随机变量，并将它简记为X。2.1随机变量及其分布称为随机变量X的概率分布函数或分布函数55/175随机变量的分布函数与性质随机变量的全部可能取值是有限个或可列无限多个，称X为离散型随机变量。随机变量的全部可能取值是不可列的，则称X为连续型随机变量56/175概率分布函数有如下性质：性质1

是x的单调非减函数，即对于，有性质2

为非负的，且满足性质3

随机变量在区间内的概率为性质4F(x)是右连续的，即随机变量的分布函数与性质57/175随机变量的概率密度函数58/175随机变量的概率密度函数59/175例2

设随机变量X的分布函数为求：（1）系数A;(2)X落在区间（0.3,0.7)内的概率；（3）X的概率密度。解：（1）由于分布函数是右连续的，所以A=1（2）（3）随机变量的概率密度函数60/175练习：随机变量X服从柯西分布求：（1）系数A和B;(2)落在区间(－1,1)内的概率；（3）X的概率密度。61/175练习：随机变量X服从柯西分布求：（1）系数A和B;(2)落在区间(－1,1)内的概率；（3）X的概率密度。解：（1）根据（2）（3）62/175随机变量数学期望与方差63/175随机变量数学期望与方差64/175随机变量数学期望与方差【例】投掷一枚骰子，出现的点数是个离散型随机变量，其概率分布为如下。计算数学期望和方差X=xi123456P(X=xi)=pi1/61/61/61/61/61/6解：数学期望为：方差为：例．若X

N(µ,σ2)，求E(X)。解：X的概率密度为：特别地，若X

N(0,1)，则E(X)=0。67/175例：均匀分布：设X在[a,b]上均匀分布，求其数学期望E(X)。

解：由定义，所以

即[a,b]上的均匀分布的期望恰为区间中点。例设随机变量在上服从均匀分布,及解根据随机变量函数数学期望的计算公式,有求例设随机变量在上服从均匀分布,及解根据随机变量函数数学期望的计算公式,有求完2.2随机过程的定义及其数字特征随机变量：取值不能预先知道，随实验结果的不同而变化。（只能事先预计，以多大的概率取什么值）过程：

确定过程：可用一个或几个时间的函数描述

随机过程：没有确定规律，不能用一个或几个确定的时间函数描述。随机过程定义

设有一簇用时间t描述的函数x(t),{x1(t),x2(t),……,xk(t),…},t∈T，如果对应于每一个固定的t1∈T，x(t1)都是随机变量，那么就称x(t)为随机过程。其中xk(t)是随机信号，为x(t)的一个样本函数，或一个实现x(t1)是随机过程x(t)的一个状态，是随机变量x1(t),x2(t),……,xk(t),…的集合称为随机信号总体随机过程的数字特征1数学期望（均值函数）设随机过程x(t)的所有样本函数有n个在t1时刻，随机信号的平均值：在任意时刻，随机信号的平均值：随机过程的摇动中心随机过程的数字特征2二阶原点矩随机过程的数字特征3方差函数（又称二阶中心矩）反映偏离摆动中心的程度随机过程的数字特征4自相关函数当t1=t2时刻画随机过程两个不同时刻的依赖关系随机过程的数字特征5自协方差函数当t1=t2时随机过程的数字特征6互相关函数度量两个不同随机过程x(t1)和y(t2)之间的相互关系随机过程的数字特征7互协方差函数描述两个随机过程x(t1)和y(t2)的线性相关程度说明两个随机过程x(t)和y(t)不相关，则实例1求，数学期望，二阶原点距，自相关函数，自协方差函数实例1求，数学期望，二阶原点距，自相关函数，自协方差函数实例1求，数学期望，二阶原点距，自相关函数，自协方差函数实例1求，数学期望，二阶原点距，自相关函数，自协方差函数实例1求，数学期望，二阶原点距，自相关函数，自协方差函数实例2其中：其中：时间序列1时间序列2计算：互相关函数、互协方差函数。判别两个时间序列的相关性实例2其中：其中：时间序列1时间序列2计算：互相关函数、互协方差函数。判别两个时间序列的相关性练习时间序列计算2.2平稳随机过程随机过程分类：

马尔科夫过程：未来状态的确定只与某一时刻有关，与所取时刻之前的特性无关。

平稳随机过程：统计特性（数字特征）不随时间t的变化而变化。（期望值恒定）

非平稳随机过程：统计特性（数字特征）随时间t的变化而变化。（期望值不恒定）2.2.1平稳随机过程定义给定随机过程x(t),如果则称x(t)为平稳随机过程。90/175解例题91/17592/17593/175课堂练习设随机过程服从均匀分布试判别X(t)的平稳性。其中94/175解设随机过程服从均匀分布试判别X(t)的平稳性2.2.2各态历经性

求随机过程x(t)的数学期望，必须给出的x(t)n个样本函数，固定时间t，求算术平均。特点：固定t=t1，进行大量的观测实验，获得近似值实际工程应用存在困难2.2.2各态历经性平稳随机过程数学期望是常值数字特征与时间无关自相关函数数字特征与计时起点无关2.2.2各态历经性平稳随机过程数学期望是常值数字特征与时间无关自相关函数数字特征与计时起点无关2.2.2各态历经性多个样本在空间上的平均一次抽样对时间的平均代替各态历经性定义平稳随机过程的每条实验曲线（样本函数实现）x1(t),x2(t),……,xk(t),…的统计特性是彼此相同的，随机信号总体的统计性质可以用一条记录曲线的统计性质来表示，这叫各态历经性。2.2.2各态历经性一个样本计算平稳随机过程的数学期望和自相关函数数学期望自相关函数两个平稳随机过程之间互相关函数2.2.2各态历经性判断

是否具有遍历性，其中均匀分布于(0，2)。解：即X(t)具有遍历性。2.2.2各态历经性103/175103104/175（3）所以X（t）为平稳随机过程105/175(4)106/175(5)(6)X(t)具有各态历经性，或称为该随机过程具有可遍历性各态历经定理的重要价值在于它从理论上给出了如下保证：一个平稳过程X(t)，若0<t<+∞，只要它满足各态历经性条件，便可以根据“以概率1成立”的含义，从一次试验所得到的样本函数x(t)来确定该过程的均值和自相关函数。估计m与的方法1得样本函数的N个值将上面的积分表示为和式23根据这两个估计式，可以算出各不同数值时相关函数的一系列近似值，从而可以作出相关函数的近似曲线。109/175109作业2.3控制系统的时间响应分析

时域分析是指在时间域内对系统的性能进行分析，是通过系统在典型信号作用下的时域响应，来建立系统的结构、参数与系统的性能的定量关系。

系统的时域响应由两部分组成：瞬态响应和稳态响应（这是从稳定性角度分析）常用的典型输入信号常用的典型输入信号常用的典型输入信号常用的典型输入信号常用的典型输入信号一阶系统的瞬态性能分析典型二系统的瞬态性能分析典型二系统的瞬态性能分析单位阶跃响应函数为典型二系统的瞬态性能分析单位阶跃响应函数为典型二系统的瞬态性能分析单位阶跃响应函数为典型二系统的瞬态性能分析单位阶跃响应函数为此时二阶系统就近似于一个惯性系统。

单位阶跃的瞬态响应无超调，无振荡，过渡过程比临界阻尼时长。

对于不允许产生振荡的控制系统，应工作在过阻尼状态，它的瞬态响应指标类似一阶系统。系统辨识与建模

第三章线性系统辨识经典方法3-1引言3-2

阶跃响应法3-3脉冲响应法（重点）3-4频率特性法3-1引言

系统辩识是随着控制理论的发展而发展的，早期的辨识工作是与经典控制理论的发展相适应而展开的。而近代系统辨识的发展也是在近代控制理论发展的基础上提出来的。

在经典控制理论中，所分析研究的是单输入、单输出系统，经常用到的系统模型是频率响应、传递函数等。所以，早期系统辩识工作的主要内容也就是寻求描述单变量系统的频率特性、传递函数，常常建立系统的非参数模型，也就是用曲线或一组采样值来表示系统的特性，或是用一些实验数据，求出系统的参数模型，用微分（差分）方程，频率响应函数或传递函数来描述系统。发展到目前，己有许多不同的辨识方法，也有不同的分类方法，从所研究对象模型形式分，可分为：经典辨识法（常指非参数模型辨识）要求对象：线性、未知结构、适于复杂系统。阶跃响应法脉冲响应法频率响应法相关分析法现代辨识法（常指参数模型辨识）要求：先有结构再辨识参数，若无结构则需先辨识结构。最小二乘法梯度校正法极大似然法

如果系统不是线性系统，则考虑在工作点附近近似看成为线性的加以考虑。这一章我们主要先从经典的辨识方法或者说从非参数模型辨识方法入手，来了解系统辨识的方法。

对非参数模型的辨识方法（或称经典辨识法），所获得的系统模型是非参数模型，它是在假设系统为线性的前提下，不必事先确定模型的具体结构。前者获得时域响应，后者获得频域响应。

在经典辨识时，为了进行系统辨识，必须对被识别对象进行激励，激励系统时所采用的信号称为测试信号，而测试信号可分为两大类。非周期测试信号：脉冲信号、阶跃信号、斜坡信号；周期测试信号：正弦信号。所谓经典的辨识方法即：通过给系统加一信号，观察其响应。对同一系统根据不同信号观察其响应；观察同一信号作用在不同系统中的响应，从而找到系统的特性规律。现代辨识方法现代辨识方法，是预先假设一种模型的结构，通过极小化模型与系统之间的误差准则函数，确定模型的参数。3.2阶跃响应法由于阶跃响应曲线与经典控制理论中对控制系统提出的时域性能指标有直接联系，如：线性连续系统的时域特性主要有两种：阶跃响应曲线冲激响应曲线上升时间峰值时间超调量调整时间因而更为常用。在被识对象上人为地施加一个已经确定的阶跃性扰动，测定出对象的响应，一定是个随时间而变化的曲线，然后根据该响应曲线，推求出被识对象的传递函数，这就是阶跃响应法。被测对象阶跃响应法，即为：阶跃响应法的辨识流程为阶跃响应法的辨识流程为阶跃响应曲线的获取方法有两种：阶跃响应直接测定阶跃响应间接测定3.2.1阶跃响应直接测定(1)当系统有非线性时，在工作点附近建立阶跃响应模型：对于一般非电量对象，可在线测试阶跃响应曲线，即将扰动u(t)加到稳态输入U0上，记录输入、输出值，直到被测对象进入一个新的稳定状态为止。(U0,Y0)为工作点，认为在工作点附近为线性的。幅值为A(2)在线实验，用户不希望轻易加信号。若向上的幅值和向下的幅值不同得多，可认为非线性很强，可将幅值再减小点。需反复实验：若噪声为白噪声，多次后均值为0。∴加入信号的幅值A要尽可能小，应在系统或质量允许范围内，但又不能太小否则可能被噪声淹没。3.2.2阶跃响应的间接测定

对于不具有自平衡能力的对象，（如：包含有积分环节等），若仍用阶跃输入作测试信号，则会引起输出量的持续变化，以致超出正常运行允许的界限。这时经常采用加方波信号，（实际上加的方波信号可以看成一个正阶跃和一个带延时的负阶跃的叠加）。若加上单位幅值方波，则有：△△的响应曲线分段作图方波信号的响应曲线（实验测量的曲线）进一步可得例题：以方波信号为输入信号激励某系统，采集到该系统在采样时间点的输出为同样也可从阶跃响应曲线上求得参数n,

k,T,τ,b。若系统为：或：在进行阶跃响应辨识法建模的过程中，要注意以下的测试要点：（1）使对象在某一稳态下运行，突然改变其输入量，使它达到另外一个稳态。阶跃扰动的取值范围为：额定值的5%-20%。（2）一般应反复测试几次，取平均值，以消除干扰带来的偶然性误差。并且应在不同工况（最大、平均负荷、最小负荷）下，重复实验。（3）在同一平衡工况下，施加反向阶跃输入信号，由正、负向响应特性的比较，可检验对象的非线性特性。（4）准确记录加入阶跃输入的计时起点。3.2.3由阶跃响应求系统的传递函数由阶跃响应确定系统的传递函数假设传递函数的结构已知利用少量特征参数确定传递函数的参数试探法阶跃曲线形状分析：一阶惯性环节、二阶惯性环节、具有纯滞后的一阶惯性环节确定传递函数的参数，并检查数据拟合情况如拟合不理想，则重新试探y(t)y0u0u(t)0.63y0/u0T2T3T0.63y0/u0处的t值为时间常数T值。(1)如:有滞后的一阶惯性环节的传递函数为：需要确定参数:k,T,τΔu为输入的幅值(2)取t1，t2

时刻的y(t1)和y(t2)，且有t2>t1>τ,则有：3.2.3由阶跃响应求系统的传递函数3.2.3由阶跃响应求系统的传递函数二阶环节阶跃响应：3.3脉冲响应法

这种方法是用施加狭脉冲信号，近似地获取冲激响应曲线。

冲激响应的双边傅里叶变换是频域传递函数或系统频域响应，且任何信号经过线性时不变系统，其输出为输入信号与脉冲响应的卷积，因此，脉冲响应法是最常用的系统辨识方法。3.3脉冲响应法但是，在工程上是不可能现实的，也就是说不可能获得准确的冲激响应h(t)。但是可以证明，只要用一个波形面积为1，且持续时间极短的（<<线性系统的调整时间Ts）狭脉冲，就认为获得了近似的冲激响应，输出信号h(t)。

如果输入尖脉冲信号(如用锤子以适当的力量敲击)，其波形面积有可能不等于1(S≠1)，则只要把所得到的输出信号y(t)÷S，同样可以认为得到了h(t)。当然，理论上直接测定线性系统的冲激响应h(t)，要求系统加入单位冲激信号，也就是作用时间为无限短，而幅值为无限大。1t脉冲响应也可以由阶跃响应获得：为采样时间，应充分小。

脉冲响应法分直接方法和间接方法。

直接法的辨识过程与阶跃响应法相同，其区别是在系统的输入端输入单位脉冲信号，在输出端测量的是脉冲响应值，最终通过脉冲响应曲线实现系统的辨识。

间接法是在系统的输入端输入具有各态历经性的随机信号，利用随机过程的相关原理，实现线性系统的辨识。

本节将主要介绍间接辨识法--相关分析法。相关分析法的优势：可以充分激励系统的动态特性抗干扰能力强辨识精度高可在线辨识相关分析法识别对象动态特性的原理，是建立在随机过程理论基础上，这里的随机过程是平稳随机过程，且具有各态历经性。相关分析法讲解思路：理论基础：维纳-霍夫方程简化输入信号为白噪声伪随机二位式---M序列怎么产生M序列性质相关系数怎么计算M序列辨识的步骤对于单输入、单输出系统（SISO），设x(t)是输入变量，y(t)是输出变量可以用传递函数(频域)或卷积公式(时域)描述根据线性系统理论，输入和输出之间的因果关系传递函数卷积公式或输入x(t)是随机过程，则y(t)也是随机过程探讨：输入随机过程x(t)的自相关函数Rxx(t)输入x(t)和输出y(t)的互相关函数Rxy(t)关系？3.2.1维纳-霍夫方程SISO系统可由下图表示：

g(τ)为系统的脉冲响应函数，即为我们需要求解的。依据线性系统的卷积定理有：

用x(t)乘以上式，有将t用代替设x(t)为均值0的平稳随机过程，则y(t)亦为均值0的平稳随机过程。任取时刻t的卷积公式为取时间的平均，即上式即为维纳-霍夫方程。它是相关分析法识别线性对象动态特性的重要依据。维纳-霍夫方程与卷积公式对比维纳-霍夫方程给出了输入的自相关函数，输入、输出的互相关函数、脉冲响应函数，三者之间的关系

由上面我们可以看出，一个具有脉冲响应函数为g(t)的被识对象，如果其输入量是信号x(t)的自相关函数,则其响应就等于输入信号x(t)与相应的输出信号y(t)

之间的互相关函数。重写维纳霍夫方程：若方程中Rxy(·)及Rx(·)已知，则解上述方程可得g(τ)但一般情况下，上述方程极难求解。只有在某些特殊情况，维纳霍夫方程才可解。特殊情况：当x(t)为白噪声信号时，有代入维纳霍夫方程后，可得可见，

g(τ)的求解，只需计算Rxy(·)即可。若观测时间Tm充分大，则有由于X,Y是记录的数据序列，则有为了便于计算互相关函数，进行以下变换：

如果输入为白噪声，只要测出输入信号与输出信号之间的互相关函数Rxy(τ)，就可求出被识对象的脉冲响应函数g(τ)。

且白噪声的能量分布在很广的频率范围内，不影响生产过程的正常运行，不会使被测对象过分偏离正常状态。3.2.2伪随机二位式序列

白噪声作为测试信号存在缺陷

首先，白噪声在物理上不易实现。

其次，是要得到精确的互相关函数和自相关函数，需要无限长时间的观测数据及长时间的积分，这不仅难以做到，而且会产生信号零点漂移，记录仪器零点漂移等新问题。

在实际工作中常采用近似的白噪声，即：伪随机噪声，其自相关函数近似于白噪声的自相关函数。3.2.2伪随机二位式序列定义：伪随机信号是一种均值为0、自相关函数、且具有重复周期T的平稳随机过程。即有：

0±T±2T±3T伪随机噪声由白噪声截断而来，是一个周期性信号。伪随机噪声及自相关函数计算自相关函数和互相关函数，设可见计算Rxx(τ)、Rxy(τ)只需计算一个周期即可。是否一定能成立？不一定成立，但是可近似

选择合适的周期T，仍能近似保持这一关系，精度可满足工程需要。伪随机噪声信号作为输入信号，则有：选择适当的截断周期，使g(τ)在τ<T时已衰减至零。得到了与白噪声作为输入的相同辨识结果。是多个脉冲响应之和则

伪随机信号具有周期性，它的产生方法有多种，最简单的方法就是将一个随机信号取其中一段（长度为T），使之在其他时间段内按该段重复，直至无穷。

实践中采用较多的是二位式伪随机序列，它的赋值不是随机的，产生和计算都很方便。离散白噪声：连续白噪声等间隔采样而成的随机序列。显然，该噪声具有连续白噪声相同的统计特性，即二位式：离散随机变量取值只有两种数值，没有第三种取值情况。例如，-1和1.离散二位式白噪声：序列中元素一般取为1和-1，也是真正意义上的白噪声，用其作为输入信号，辨识g(τ)的结果与连续白噪声的结果完全一致。例：某离散二位式噪声

1111-1-1-11-1-111-11-1···主要性质:(1)-1和1出现的次数相等；

(2)总游程数为(N+1)/2,且-1和1出现的游程相等，最多相差1个。(N为序列长度)

(3)其自相关函数为注：状态“1”和“-1”连续出现的段叫游程。176/175例如m序列:

000111101011001,p=15共有８个游程：长度为４的游程有一个：即“1111”长度为３的游程有一个：即“000”长度为２的游程有两个：即“00”和“11”长度为１的游程有四个：即“0”和“１M序列举例

二位式伪随机信号可以人为地产生和复制，比连续式信号容易产生和重现，在计算互相关函数时，可以将乘法简化为取正反向的运算。最常用的二位式伪随机序列为M序列。

M序列信号可以由专用多级脉冲触发器组成的电路产生，也可以由在线控制计算机产生。

M序列生成电路由N级脉冲触发器组成。其中，在K0级脉冲触发器的输出处给出反馈信号，与第n级的输出进行异或（即模2和）后，送入第1级触发器，下图为其原理图：a11时钟脉冲置初始状态X0X1+a22X2+a33X3+an-1Xn-1+nXn图中：X0,X1,……,Xn-1,Xn为二元序列，是n个脉冲触发器的输出，只要n个触发器的初始值不均为0，则随着时钟脉冲的逐个输入，在n级输出处会给出一个一定规律的循环交变的信号序列。M序列中各元素之间满足下式关系：这个序列为二进制逻辑，0、1交替给出,所对应的电路信号为-V和+V,故又称双电平信号。2）为模二和，即逻辑异或（相同为0，不同为1）；3）只要a1、a2、a3

、……、an-1选择适当，即可得到一个最大长度的循环序列，故称为最大长度伪随机信号序列——

M序列。M序列的周期长度为N=2n-1。4）任意一个脉冲触发器输出均为M序列。其中：

1）a1、a2、a3

、……、an-1为反馈系数,取值为0或1,若取0

则表示此输出无反馈,取1则为有反馈。异或（模２和）相同为0，相异为1

【例】设一个M序列由4个脉冲触发器组成（即N=4），在第3个脉冲触发器输出端引出反馈（即K0=3）。下面为产生这个M序列的电路图和它的时序变化，设对4个触发器预置初值为全“1”。11111234时钟预置初值输出123456789101112131415161718192011110111001100011000010000101001110001101011010110101101111011110111001100011000

【练习】设一个M序列由4个脉冲触发器组成（即N=4），在第3个脉冲触发器输出端引出反馈（即K0=3）。下面为产生这个M序列的电路图和它的时序变化，设对4个触发器预置初值为全“1”。10111234时钟预置初值u1u2u3u4输出u(t)123456789101112131415161718192010110101101011011110111101110011000110000100001010011100011010110101101000011000

【练习】设一个M序列由4个脉冲触发器组成（即N=4），在第2个脉冲触发器输出端引出反馈（即K0=2）。下面为产生这个M序列的电路图和它的时序变化，设对4个触发器预置初值为全“1001”。10011234时钟预置初值u1u2u3u4输出u(t)123456789101112131415161718192010011100111011110111001110011100111011110111001110011100111011110111001110011100程序实现M(0)M(1)M(2)M(3)M(4)伪代码：

在实际应用中，可根据M序列周期长度的不同，选择移位寄存器的级数n,只要适当地选择反馈节点，便可以得到期望的序列长度N。在移位寄存器中，决定反馈通道是关键。（即ai的取值）a11时钟脉冲置初始状态X0X1+a22X2+a33X3+an-1Xn-1+nXn二电平M序列及其性质工程实际：将M序列转变成电平信号，“0”取为av,“1”取为-av。移位脉冲周期为Δ，则该二电平M序列的周期为NΔ。性质：(1)周期长度N=2n-1（反映在时间为T=N×△）(2)在一个周期中，1的状态个数为0的状态个数为二电平M序列及其性质(3)游程总数为

显然，当N→∞时，mx→0.(4)均值mx在一个周期NΔ内，其均值mx为(5)自相关函数Rx(τ)

不加证明给出Rx(τ)的计算结果。当N很大时，上式可写为Rx(τ)的波形如下：从图中可知，Rx(τ)由两部分构成：

三角脉冲分量，记为Rx1(τ);

直流分量，记为Rx2(τ);则有：其中:的波形如下：当Δ很小时，Rx1(τ)可认为是脉冲函数，则有可见，M序列具有白噪声序列的数字特性。3.6.4L序列的产生及性质M序列的不足：信号中含有直流分量，给互相关运算的数据处理带来不便；辨识结果难以补偿被测系统运行中缓慢的非随机性干扰.解决策略：逆M序列----L序列L序列产生的方法011110010000113.2.3.用M序列辨识线性系统的脉冲响应

g(τ)以M序列信号为输入，探讨g(τ)与Rxy(τ)的关系当Δ很小时，Rx1(τ)可认为是脉冲函数，有则有：记：，C为常数3.2.3.用M序列辨识线性系统的脉冲响应

g(τ)Rxy(τ)可根据输入输出数据序列计算：可见，只需将Rxy(τ)曲线向上平移C，即可得g(τ)。C由Rxy(τ)的稳态值确定，一般用目测法就能估算出来。其图形如图示。C(2)公式法求g(τ)两边积分有(此处T=N△)代入(1)式，可解得g(τ)为公式法求g(τ)公式组上述公式组，每次计算可得到g(τ)的一个离散值。测试原理图这种线路称为“串行法”，每改变一个延时，得到脉冲响应曲线的一点，而每测一点，输入的序列至少要有一个周期以上的长度。取得一条曲线，要经过许多次试验，颇费时间。并行辨识法系统乘法器乘法器白噪声正常输入r(t)设定积分时间

由脉冲响应求传递函数，方法很多，在这里只举几种常用的方法。

一阶系统，其传递函数为：其中k和T可从曲线中获得：3.2.4基于脉冲响应曲线的系统辨识

2）二阶系统，其传递函数可描述为：其中ξ和ω0可从脉冲响应曲线上直接获得：3)一种通用求法：设系统的传递函数为：当特征方程：当有n个单根对应的脉冲响应：时，则传递函数可写成：当特征方程式具有重根时，传递函数可写成：为单根，s0

为r阶重根，则对应的脉冲响应为：利用脉冲响应获得来确定3.2.5利用M序列辨识的步骤1）估计系统过渡过程时间（也称调整时间）

和系统的最高工作频率

作为选择M序列参数的依据。过渡过程时间通过简单的阶跃响应曲线来估计，最高截止频率则通过给系统施加不同周期的正弦信号或方波信号，观测系统输出来估计。2）精心选择M序列的参数

，当系统频率特性接近低通滤波特性时，M序列的参数t应满足：3）M序列的循环周期必须大于系统的过渡过程时间，以保证时间大于

后脉冲响应衰减接近于零。4）

确定M序列的幅值。M序列的幅度a不能选择过大，以免系统进入非线性区或影响系统正常生产，但也不能过小，以保证一定的信噪比。3.2.5利用M序列辨识的步骤5）采集数据，当M序列刚加上时，系统输出在一段时间内是非平稳的，一般从第二个循环周期开始采集数据。令

，并将

记为

，系统输出

记为6）计算互相关函数。如果测得r+1个周期的采样值，用后r个周期的采样值进行计算：考虑到M序列为二电平伪随机序列，则7)计算脉冲响应值（一）第一次预备试验

这次试验的目的是粗估一下被测系统的调整时间Ts，以及最高工作频率fM。1.粗估TS

在系统正常运行工况下，叠加一个窄的脉冲，或幅度不大的阶跃扰动，用示波器或记录仪观测其响应幅值的变化，注意记录Ym和0.1Ym的t值。在实际工程应用中，在正式进行系统辨识时，需先做两次预备实验，以获取用于参数辨识的相关参数。窄的脉冲幅度不大的阶跃2.测估最大工作频率fM

可以借用专用仪器，或在线计算机系统施加不同周期的低频f低振幅的矩形波信号，用示波器或记录仪观察对矩形波扰动信号的响应,，最后确定该系统的最大工作频率fM。

随着信号频率f增加，响应曲线震荡会逐渐减弱…，直到刚好看不出振荡时，也就是系统对此激励不再反应了，此时所对应的频率为fM。f<<fMf<fMf>fM确定Δ(移位脉冲周期)和N

(序列长度)（二）选择参数适合的PRBS序列作为测试信号

为获得良好的辨识效果，测试信号应能充分激励起被测系统在其通频带范围内所有频率下的动态行为。也就是要求测试信号的频谱能全部覆盖待辨识系统的全部重要工作频率。

为了满足这个要求，往往需要选用短的时钟周期Δ。但是，从

序列频谱Su(τ)上谱线高度接近一致的区段(即ω=0～2π/3Δ)来看,应完全覆盖系统的通频带（即ω=0～2πfM）。

另外,为使PRBS

序列相对于被辨识的系统来说,可以近似地当作白噪声,因此应取：当然，Δ越小M序列信号的自相关函数Ru(τ)

的图形越接近于冲激函数序列，辨识结果便会越接近冲激响应。但是，Δ也不宜太小，在M序列幅值a

受到限制时，过宽的频率范围，就会减少重要的频率区间的有效功率。为使所得冲激响应g(t)能够在T=NΔ之内结束，防止Ru(τ)曲线上出现重叠现象，∴要求T=NΔ>Ts

。即（Ts为调整时间）通常取

如果系统的过渡过程可以看出有振荡，但在频率再升高，振荡减少后又有回升现象，则说明系统中存在一个或多个谐波。这种情况往往要记录Δωb值。则应选择Δ值使之在每个峰内的Δωb有两条以上的功率谱线，即：上述多条要求，最好均能兼顾到。（三）第二次预备试验在Δ和N已选定的基础上，再作一次预备试验来确定a的大小。为了提高信噪比，应选取足够大的M序列信号电平a。但是，由于不能影响系统的正常运行（如：超出线性范围，或产品误差超出生产允许的限度），故a又不能取得过大。试验中将选定的Δ和N=2n-1参数的M序列施加到系统的输入端，令a由小逐渐变大，直到│y(t)│<允许的限度为止，以确定适合的a值。y(t)a被测系统（四）正式试验将参数为Δ、N及a的M序列施加到被测系统上，同时采集k、u(k)、y(k)=Y0+y’(k)+η(k)，k=0,1,2,3,…,γ。要知道，M序列u(k)

实质上是周期信号，施加到被测系统后，系统输出端产生的响应

y(k)

也将是周期波形。不过在t=0～Ts的过渡过程初始阶段内，会出现一些非周期衰减分量,使

y(k)

不能满足平稳的特性，也就是不能满足:均值=常数均方差=常数相关函数与t无关与τ有关∴正式试验中在施加M序列扰动后至少一个周期再开始记录数据k,u(k),y(k)。为了提高辨识精度，在预扰动之后要测取γ=3～4个整周期的数据。小结：粗估Ts（系统的调整时间）和fM（最高工作频率）；一、利用M序列作为输入信号辨识被测系统的脉冲响应的步骤：选择M序列的参数Δ和N，一般满足在Δ和N基础上经过实验确定M序列的幅值a；施加M序列，在第一个周期（N·Δ>TS）之后采样k、x(k)、y(k)计算出获得（通常取）举例：用相关分析法测定和计算常压加热炉的炉膛温度与燃料量之间的动态关系。3.2.6相关分析法的抗干扰性g(t)

x(t)ξ(t)Y’(t)y(t)+

+可测输入过程噪声可测输出则x(t)与y(t)的互相关函数为：因为和不相关，故有即有：总的来说，如果输入为白噪声，则有：

结论：相关分析法不受噪声的影响（白噪声）或者说容易去除噪声（具有平稳特性的噪声）的影响，即相关分析法具有抗干扰性。

采用相关法测定（或称辨识）对象的动态特性的主要优点是：(1)辩识试验可以在被识对象正常的运行状态下进行;(2)只要输入测试信号与干扰噪声是统计独立的，辨识工作可不受干扰影响；(3)相关法辨识不需要有关于被识对象的验前知识。3-3

频率特性法频率特性是描述动态系统的非参数模型频率特性的测取正弦波法矩形波法

给定为施加任何一种周期信号u(t)波形时，均可获得频率响应，但是必须在系统已经进入稳态之后，再测量输出y(t)的值，才有意义。（否则输出中将含有非周期过渡过程含量，致使结果不准）。从频率特性得到系统传递函数幅相特性对数频率特性3-3

频率特性法单一正弦波法在待测系统输入端加上某个频率的正弦信号，记录输出达到稳态后输出的振荡波形对于线性系统，得到的是一个与输入同频率的、但幅值与相位发生变化的正弦波，根据幅值比和相位移，可得到线性系统的频率特性使用正弦波法可测出系统的带宽，当增加输入正弦信号频率至

max时，系统输出幅值将趋近于零此法对缓慢响应过程非常费时，此时可利用线性系统符合叠加原理的特点，采用组合正弦信号3-3

频率特性法组合正弦波法在待测系统输入端加上频率、幅值均已知的组合正弦波在稳态下测取输出组合波，再利用傅氏变换对输出组合波作分解3-3

频率特性法矩形波法若难以产生正弦波，则可以使用矩形波输入信号输入信号的傅立叶级数分解：高次谐波可以忽略3-4

频率特性法幅相特性求传递函数一阶惯性环节3-4

频率特性法幅相特性求传递函数二阶环节幅相特性分布在两个象限φ'=π/2,r(w')≤r(0)/23-4

频率特性法幅相特性求传递函数带纯滞后的一阶惯性环节3-4

频率特性法对数频率特性求传递函数一阶环节或一阶滞后环节3-4

频率特性法对数频率特性求传递函数二阶环节或二阶滞后环节1经典辨识算法包含哪几种？2给出脉冲响应辨识法的辨识步骤3相关分析法的数学理论基础是什么？基于该理论基础，最终是为了获得什么？4维纳霍夫方程是什么？其输入是什么？输出是什么？在输入为白噪声时，维纳霍夫方程的形式是？5给出基于M序列的辨识过程。系统辨识与建模第4章基于最小二乘法的参数辨识

本章的学习目的1、掌握最小二乘参数辨识方法的基本原理2、掌握常用的最小二乘辨识方法和极大似然法3、熟练应用最小二乘参数辨识方法进行模型参数辨识4、能够编程实现最小二乘参数辨识4.1输入输出模型随机模型确定性模型1.确定性模型

)(kG

)(ku)(kyn阶差分方程描述：引进单位延迟算子设多项式

则差分方程可写为1.确定性模型

)(zG

)(ku)(ky系统脉冲传递函数描述：如果u(k)和y(k)的初始条件为零，即当。设y(k)和u(k)的变换为Y(z)和U(z)，则（4-1）式的z变换为

也可以写成时2.随机模型观测值可表示为：整理得：2.随机模型相当于采样m组数据，得到2.随机模型例题：通过试验确定热敏电阻阻值和温度间的关系

当测量没有任何误差时，仅需2个测量值。每次测量总是存在随机误差。根据最小二乘的准则有根据求极值的方法，对上式求导4.2最小二乘法批量处理方法单输入单输出(SISO)随机模型的n阶差分方程为

)(kG

)(ku)(kx)(kv)(ky单输入单输出(SISO)随机模型的n阶差分方程为引入向量表示形式：系统输入输出的最小二乘格式为

现在分别测出个输出值和输入值：及。则可写出N个方程：上述N个方程可写成下列向量－矩阵形式4.2.1一般最小二乘算法

式中为N个输出值组成的向量；所组成的n维向量所组成的维向量；所组成的N维噪声，即或

为输出值所组成的阵块；为输入值所组成的矩阵块。即一般最小二乘的矩阵形式为式中

为维测量矩阵，为维参数向量。因此，是一个含有个未知参数的N个方程组成的联立方程组。如果，则方程组是不定的，不能唯一地确定参数向量。如果，则当测量误差时，就能准确地解出参数向量，即如果测量误差不等于零，则

从上式可看出，随机测量噪声对参数的估计值有影响，为了尽量减小对的估值的影响，应该取，即方程数目大于未知数数目。在这种情况下，不能用解方程的方法求，而要采用数理统计的方法求的估值。这样可减小对的估值的影响。这种给定测量向量和测量矩阵求参数估值的问题，就是系统参数的辨识问题。式中写出差分方程估计式的某一行，得设表示的最优估值，表示的最优估值，则有设表示与之差，通常称为残差。由上式得

把分别代入

，可得残差把这些残差写成向量形式：最小二乘法估计要求残差的平方和为最小，即按照指标函数为最小确定估值。可按来求的最小二乘法估计值。即由此式用左乘等号的两边，得显然，当矩阵存在时，

才有解。如果是随机序列或伪随机二位式序列，则矩阵是非奇异的，即存在，

有解。J为极小值的充分条件是

因为有解与正定等价，所以可以保证正定来确定对输入序列的要求。则

因此，要求正定，根据正定矩阵的性质，必须保证正定。这个条件称为阶持续激励条件。通常，输入序列采用随机序列或M序列时，它们都满足这个持续激励条件。显然，若为常值序列时，为奇异阵，不满足持续激励条件。4.2.2加权最小二乘算法4.2.2加权最小二乘算法4.2.2加权最小二乘算法4.2.2加权最小二乘算法4.2.3正则化最小二乘算法4.2.3正则化最小二乘算法例4.2考虑仿真对象选择如下的辨识模型进行一般的最小二乘参数辨识。

Matlab程序为：clearclcu=[-1,1,-1,1,1,1,1,-1,-1,-1,1,-1,-1,1,1];%系统辨识的输入信号为一个周期的M序列

y=zeros(1,16);%定义输出观测值的长度

fork=3:16

y(k)=-1.5*y(k-1)-0.7*y(k-2)+u(k-1)+0.5*u(k-2);%用理想输出值作为观测值

end

subplot(3,1,1)%画三行一列图形窗口中的第一个图形

stem(u)%画出输入信号u的经线图形

subplot(3,1,2)%画三行一列图形窗口中的第二个图形

i=1:1:16;%横坐标范围是1到16，步长为1

plot(i,y)%图形的横坐标是采样时刻i,纵坐标是输出观测值z,图形格式为连续曲线

subplot(3,1,3)%画三行一列图形窗口中的第三个图形

stem(y),gridon%画出输出观测值z的经线图形，并显示坐标网格

u,y%显示输入信号和输出观测信号%L=14%数据长度

Matlab程序为：

%L=14%数据长度

HL=[-y(2)-y(1)u(2)u(1);-y(3)-y(2)u(3)u(2);-y(4)-y(3)u(4)u(3);-y(5)-y(4)u(5)u(4);-y(6)-y(5)u(6)u(5);-y(7)-y(6)u(7)u(6);-y(8)-y(7)u(8)u(7);-y(9)-y(8)u(9)u(8);-y(10)-y(9)u(10)u(9);-y(11)-y(10)u(11)u(10);-y(12)-y(11)u(12)u(11);-y(13)-y(12)u(13)u(12);-y(14)-y(13)u(14)u(13);-y(15)-y(14)u(15)u(14)]%给样本矩阵HL赋值

ZL=[y(3);y(4);y(5);y(6);y(7);y(8);y(9);y(10);y(11);y(12);y(13);y(14);y(15);y(16)]%给样本矩阵zL赋值

%calculatingparameters%计算参数

c1=HL'*HL;c2=inv(c1);c3=HL'*ZL;c=c2*c3%计算并显示

%DISPLAYPARAMETERS

a1=c(1),a2=c(2),b1=c(3),b2=c(4)%从中分离出并显示a1、a2、b1、b2

%End

运行结果为：u=-11-11111-1-1-11-1-111y=000.50000.25000.52502.11254.30126.47316.19883.2670-0.9386-3.1949-4.6352-6.2165-5.5800-2.5185HL=001.0000-1.0000-0.50000-1.00001.0000-0.2500-0.50001.0000-1.0000-0.5250-0.25001.00001.0000-2.1125-0.52501.00001.0000-4.3012-2.11251.00001.0000-6.4731-4.3012-1.00001.0000-6.1988-6.4731-1.0000-1.0000-3.2670-6.1988-1.0000-1.00000.9386-3.26701.0000-1.00003.19490.9386-1.00001.00004.63523.1949-1.0000-1.00006.21654.6352