系统建模与仿真第4次课

12026/9/5(1)乘同余法首先，用递推同余式产生正整数序列{xi}，即M为2的方幂，即M=2k，k为大于2的整数；A≡3或A≡5（mod8），且A不能太小；初值x0取正奇数。再令则{ξi}是伪随机数序列，循环周期可达2k-2。回想22026/9/5(2)混合同余法首先，用递推同余式产生正整数序列{xi}，即M为2的方幂，即M=2k，k为大于2的整数；A=2n+1，其中2≤n≤34；初值x0非负整数，C为正整数。再令则{ξi}是伪随机数序列，循环周期可达2k。32026/9/5线性系统g(τ)正常输入X(t)y(t)+yw(t)延迟τ乘法器积分器Kg(τ)白噪声Xw(t)Xw(t-τ)含有正常输入时的系统辨识模拟方块图42026/9/5M序列的产生及性质M序列是伪随机二位式序列的一种形式，它含有白噪声的性质，不仅能够确保有较好的辨识效果，并且工程上又易于实现。X1X2X3X4移位脉冲XOR输出52026/9/5M序列的性质（1）由n级移位寄存器产生的M序列的最大周期为N=2n-1。（2）M序列中，状态“0”或“1”持续出现的段称为游程。游程中“0”或“1”的个数称为游程长度。由n级移位寄存器产生的M序列的游程总数2n-1，“0”“1”各二分之一;并且长度为1的游程占总数的1/2，有2n-2个;并且长度为2的游程占总数的1/4，有2n-3个;以这类推，长度为i（1≤i≤n-2）的游程占总数的1/2i，有2n-i-1个;长度为n-1的游程只有1个，为“0”的游程;长度为n的游程只有1个，为“1”的游程;11110001001101062026/9/5（3）全部M序列均含有移位可加性，即2个彼此移位等价的相异M序列，按位模2相加所得到的序列仍为M序列，并与原M序列等价。10011010111100010011010110101111000100110101111100010011010111100072026/9/5则式（2.87）可写为

（2.88）

N维输出向量2n+1维参数向量N维噪声向量N×(2n+1)维测量矩阵最小二乘法：82026/9/5（2.94）

最小二乘估计要求残差的平方和为最小，即按照指标函数

（2.95）

为最小来确定估值。求J对的偏导数并令其等于0可得的最小二乘估计

（2.98）

J为极小值的充分条件是

（2.99）

即矩阵为正定矩阵，或者说矩阵是非奇异的。

92026/9/5[例2.1]已知某一单输入单输出线性系统的差分方程形式为但其参数，，为未知数，且为不相关的随机序列。经过辨识试验，测得5组输入输出数据为102026/9/5试求出其最优参数预计。112026/9/5[解]

令最优参数估计为，令输出的最优估计为。测量矩阵为122026/9/5该矩阵的转置为两者之积为132026/9/5的特征值为，，。由于它的特征值均为正数，所以为正定矩阵，满足残差二次型取最小的充分条件，其中。142026/9/5根据残差二次型取最小的必要条件可得最优参数估计为152026/9/5矩阵的逆为于是162026/9/5最后求得即最优参数估计为172026/9/52.5.1.3最小二乘辨识中的输入信号问题

当矩阵的逆阵存在时，式（2.98）才有解。一般地，如果是随机序列或伪随机二位式序列，则矩阵是非奇异的，即存在，式（2.98）有解。现在从矩阵必须是正定的这一要求出发，来讨论对的要求。在这里为了方便起见，假定是均值为0的随机过程。

可以推出矩阵为正定的必要条件是：为持续激励信号。（推导过程略）

随机序列或伪随机二位式序列都可以作为测试信号。

182026/9/5由于M序列对时间是离散的，而输入需要对时间持续，因此在实际应用中，总是把状态为“0”和“1”的M序列变换成幅度为+a和-a的二电平序列，其中“0”对应高电平+a，“1”对应低电平-a。这种对时间持续的序列称为二电平M序列二电平M序列的产生11110001001101002△+a-a6△16△t192026/9/52.5.1.4最小二乘预计的概率性质如果ξ(k)是不有关随机数序列，且均值为0。1）无偏性2）一致性3）渐进正态性性辅助变量法、广义最小二乘法如果ξ是均值为0且服从正太分的白噪声向量，则最小二乘参数预计值服从正态分布。202026/9/52.5.2一种不需矩阵求逆的最小二乘法

设系统的微分方程模型为

（2.141）

令

（2.142）

（2.143）

212026/9/5则式（2.141）可以写为

（2.144）

取

222026/9/5则有

（2.145）

上式中矩阵的阶数越大，所包含的信息量就越多，系统参数估计的精度就越高。为了获得满意的辨识结果，矩阵的阶数常常取得相当大。这样，在用式（2.146）计算系统参数的估计值时，矩阵求逆的计算量很大。本节介绍一种算法来代替矩阵求逆，在不降低辨识精度的前提下，可以使辨识速度有较大提高。具体算法如下。

系统的最小二乘辨识结果为

（2.146）

232026/9/5首先设系统的阶次为0，则有242026/9/5则和均为常数，即

（2.147）

252026/9/5由式（2.146）可得

（2.148）

设

（2.149）

262026/9/5若系统阶次为n时已经求出，则系统阶次数为n+1时有

（2.150）

式中

（2.151）

式中：为列向量；为一标量。

272026/9/5由分块矩阵求逆公式可得（2.152）

式中

282026/9/5则

这时，仿照上述方法容易求出

式中

292026/9/5(2.153)302026/9/5这样，就能够按照式（2.146）辨识出阶次为n+1时系统的参数。由于这一过程只涉及矩阵相乘和矩阵与向量相乘等运算，因此计算量较小，而矩阵求逆的精度不变。因此说，本节算法在不损失辨识精度的前提下提高了辨识速度，这一算法特别合用于阶次未知状况下的系统辨识。312026/9/5为了实现实时控制，必须采用递推算法，这种辨识办法重要用于在线辨识。另外，某些自适应控制办法，例如自校正控制办法也要用到递推最小二乘法。其基本原理是一边辨识一边控制，循环递推。2.6递推最小二乘法分析

递推最小二乘算法是动态系统实时辨识中使用得最多的。对其算法进行进一步研究含有较大的实用意义。322026/9/5图4.1动态系统递推最小二乘在线辨识过程原理图332026/9/5首先将式（2.88）改记为其中342026/9/5如果再获得一对新的观测值，则有下面的递推公式这里无视其推导过程，直接给出公式。352026/9/5（2.154）（2.155）（2.156）递推最小二乘算法362026/9/5其中为标量，它表示在已测得的基础上又多测得一个输出。另外表示在的最后一行再补充一行，所以其列数与的列数相等，均为2n+1列。372026/9/5（1）设为N的初始值，则可算出初值

（2）假定，c是充分大的常数，

I为单位矩阵

递推之后能得到较好的参数估计。

，则经过若干次为了进行递推计算，需要给出

和的初值

和

，有两种给出初值的办法。

382026/9/5（2.154）（2.155）（2.156）其中的公式（2.155）所表达的可以理解为递推最小二乘增益系数。公式（2.155）不具有递推性。公式（2.154）和（2.156）具有递推性。注释392026/9/52.6.1递推最小二乘法举例[例2.2]已知某一单输入单输出线性系统的差分方程形式但其参数，，为未知数，且为不相关的随机序列。实时测得4组输入输出数据为402026/9/5假设试用递推最小二乘法递推两步，分别求得和。412026/9/5[解]对于这个问题有，先令。422026/9/5于是432026/9/5442026/9/5由式(2.155)和（2.156）可得452026/9/5因此由式（2.154）得462026/9/5至此，已求得了第一步递推值。472026/9/5令，则。482026/9/5492026/9/5502026/9/5至此，已求得了第二步递推值。512026/9/52.6.2递推最小二乘法程序构造522026/9/5532026/9/5542026/9/5552026/9/52.7递推最小二乘辨识的在线算法考虑被控对象的在线辨识和控制问题。为不相关的随机其中序列。重要目的是设计状态反馈控制律，使得输出趋向于零。但主要问题是模型参数、未知。562026/9/5整个闭环系统原理以下图所示。572026/9/5

将被控对象模型变换成状态空间模型，令状态变量为：则有：这是一种一阶线性离散系统。582026/9/5暂不考虑干扰信号，设状态反馈控制律为：其中为反馈增益系数。可得闭环系统方程其特性方程为：592026/9/5解出闭环极点为：假设盼望闭环极点为：模不大于1则有：602026/9/5于是在理论上状态反馈控制律为：可见状态反馈控制律不但依赖于系统的状态，而且还依赖于模型参数、。但是这两个参数却是未知的，在这种情况下，我们通过最小二乘法辨识出这两个参数和，最后得到实际的控制612026/9/5为了仿真研究的方便，我们在解算模型的输出时，假设这两个参数为、。在假设这两个参数已知的理想情形下，我们先对闭环系统进行仿真研究。622026/9/5假设状态初值为：随机干扰的范畴为：632026/9/5编写下列M文献x0=-2.3;x(1)=x0;fork=1:1:20x(k+1)=0.3*x(k)+0.2*rand(1)-0.1;endt=0:1:length(x)-1;stem(t,x);axis([0,21,-3,1]);绘制枝干图随机干扰642026/9/5运行成果以下：652026/9/5可见在抱负状况下，闭环系统是渐近稳定的。而实际状况下，被控系统的参数是未知的，需要通过最小二乘法进行辨识。下面给出实际闭环系统仿真流程图662026/9/5开始给定真实参数、解算模型输出之用设定较大的正常数给定估计参数初值、不可为零否则控制律发散672026/9/5设定预计参数向量的初始化给定系统初始状态，即初始输出相当于682026/9/5设置设置矩阵P的初始化，相当于692026/9/5状态反馈控制律系统模型1表示1×1矩阵，即标量随机数，范围为0～1702026/9/5712026/9/5722026/9/5结束否？是结束否732026/9/5编制M文献以下：a1=2;%真实参数b0=1.4;%真实参数c=10000;%很大的正常数a1_m=0;%预计参数初值b0_m=0.5;%预计参数初值sita_m(:,1)=[a1_m,b0_m]';%预计参数向量初值y(1)=-2.3;%系统初始状态（即初始输出）z=c^2;%很大的正常数之平方P(:,:,1)=diag([z,z],0);%P阵的初始化742026/9/5forN=1:1:20u(N)=(sita_m(1,N)+0.3)*y(N)/sita_m(2,N);y(N+1)=-a1*y(N)+b0*u(N)+0.2*rand(1)-0.1;fai(:,N+1)=[-y(N),u(N)]';W=P(:,:,N)*fai(:,N+1);r=1/(1+W'*fai(:,N+1));K(:,N+1)=W*r;P(:,:,N+1)=P(:,:,N)-K(:,N+1)*W';sit