【2-1】 设某物理量Y 与X1、X2、X3的关系如下：Y=θ1X 1+θ2X 2+θ3X 3 由试验获得的数据如下表。

试用最小二乘法确定模型参数θ1、θ2和θ3 X1: 0.62 0.4 0.42 0.82 0.66 0.72 0.38 0.52 0.45 0.69 0.55 0.36 X2: 12.0 14.2 14.6 12.1 10.8 8.20 13.0 10.5 8.80 17.0 14.2 12.8 X3: 5.20 6.10 0.32 8.30 5.10 7.90 4.20 8.00 3.90 5.50 3.80 6.20 Y: 51.6 49.9 48.5 50.6 49.7 48.8 42.6 45.9 37.8 64.8 53.4 45.3解：MATLAB 程序为：Clear all;A= [0.6200 12.000 5.20000.4000 14.2000 6.10000.4200 14.6000 0.32000.8200 12.1000 8.30000.6600 10.8000 5.10000.7200 8.2000 7.90000.3800 13.0000 4.20000.5200 10.5000 8.00000.4500 8.8000 3.90000.6900 17.0000 5.50000.5500 14.2000 3.80000.3600 12.8000 6.2000];B=[51.6 49.9 48.5 50.6 49.7 48.8 42.6 45.9 37.8 64.8 53.4 45.3]';C=inv(A'*A)*A'*B=[0.62 12 5.2;0.4 14.2 6.1;0.42 14.6 0.32;0.82 12.1 8.3;0.66 10.8 5.1;0.72 8.2 7.9;0.38 13 4.2;0.52 10.5 8;0.45 8.8 3.9;0.69 17 5.5;0.55 14.2 3.8;0.36 12.8 6.2]公式中的A 是ΦN, B 是YN ，运行M 文件可得结果：在matlab 中的运行结果：C=29.59032.44660.4597【2-3】 考虑如下模型)()(3.03.115.0)(2121t w t u z z z z t y ++-+=---- 其中w(t)为零均值、方差为1的白噪声。

根据模型生成的输入/输出数据u(k)和y(k),分别采用批处理最小二乘法、具有遗忘因子的最小二乘法（λ=0.95）和递推最小二乘法估计模型参数（限定数据长度N 为某一数值，如N=150或其它数值），并将结果加以比较。

解：1、批处理最小二乘法M文件如下：clear all;close all;a=[1 -1.3 0.3]';b=[1 0.5]';d=1;%对象参数na=length(a)-1;nb=length(b)-1;%na,nb为A,B阶次N=200;%观测数据组uk=zeros(d+nb,1);%创建d+nb行1列的零列向量，给输入赋初值0yk=zeros(na,1);%输出初值x1=1;x2=1;x3=1;x4=0;%移位寄存器初值S=1;%方波初值W=randn(N,1);%产生均值为零方差为1的白噪声序列w(k)theta=[a(2:na+1);b];%theta为列向量，[-1.3 0.3 1 0.5]',为对象参数真值for k=1:Nphi(k,:)=[-yk;uk(d:d+nb)]';%phi为行向量，组成phi矩阵y(k)=phi(k,:)*theta+W(k);%采集输出数据IM=xor(S,x4);%进行异或运算，产生逆M序列if IM==0u(k)=-1;elseu(k)=1;endS=not(S);%产生方波M=xor(x3,x4);%进行异或运算，产生M序列x4=x3;x3=x2;x2=x1;x1=M;%寄存器移位for i=d+nb:-1:2uk(i)=uk(i-1);enduk(1)=u(k);for i=na:-1:2yk(i)=yk(i-1);endyk(1)=y(k);endthetae=inv(phi'*phi)*phi'*y'%计算参数估计值thetaeJ=y*y'-2*(phi'*y')'*thetae+thetae'*phi'*phi*thetae%求最小估计误差平方和clear all;close all;a=[1 -1.3 0.3]';b=[1 0.5]';d=1;%对象参数na=length(a)-1;nb=length(b)-1;%na,nb为A,B阶次N=200;%观测数据组uk=zeros(d+nb,1);%创建d+nb行1列的零列向量，给输入赋初值0yk=zeros(na,1);%输出初值x1=1;x2=1;x3=1;x4=0;%移位寄存器初值S=1;%方波初值W=randn(N,1);%产生均值为零方差为1的白噪声序列w(k)theta=[a(2:na+1);b];%theta为列向量，[-1.3 0.3 1 0.5]',为对象参数真值for k=1:Nphi(k,:)=[-yk;uk(d:d+nb)]';%phi为行向量，组成phi矩阵y(k)=phi(k,:)*theta+W(k);%采集输出数据IM=xor(S,x4);%进行异或运算，产生逆M序列if IM==0u(k)=-1;elseu(k)=1;endS=not(S);%产生方波M=xor(x3,x4);%进行异或运算，产生M序列x4=x3;x3=x2;x2=x1;x1=M;%寄存器移位for i=d+nb:-1:2uk(i)=uk(i-1);enduk(1)=u(k);for i=na:-1:2yk(i)=yk(i-1);endyk(1)=y(k);endthetae=inv(phi'*phi)*phi'*y'%计算参数估计值thetaeJ=y*y'-2*(phi'*y')'*thetae+thetae'*phi'*phi*thetae%求最小估计误差平方和运行结果如下：thetae =-1.30220.29820.97090.5091J =186.0817与真值比较参数a1 a2 b0 b1真值-1.3 0.3 1 0.5估计-1.3022 0.2982 0.9709 0.50912、递推最小二乘法clear all;close all;a=[1 -1.3 0.3]';b=[1 0.5]';d=1;%对象参数na=length(a)-1;nb=length(b)-1;%na,nb为A,B阶次N=200;%观测数据组uk=zeros(d+nb,1);%创建d+nb行1列的零列向量，给输入赋初值0yk=zeros(na,1);%输出初值u=randn(N,1);%输入采用白噪声序列w=randn(N,1);%产生均值为零方差为1的白噪声序列w(k)theta=[a(2:na+1);b];%theta为列向量，[-1.3 0.3 1 0.5]',为对象参数真值thetae_1=zeros(na+nb+1,1);%thetae的chuzhiP=10^6*eye(na+nb+1);%P初值for k=1:Nphi=[-yk;uk(d:d+nb)];%phi为行向量，组成phi矩阵y(k)=phi'*theta+w(k);%采集输出数据K=P*phi/(1+phi'*P*phi);thetae(:,k)=thetae_1+K*(y(k)-phi'*thetae_1);P=(eye(na+nb+1)-K*phi')*P;thetae_1=thetae(:,k);for i=d+nb:-1:2uk(i)=uk(i-1);enduk(1)=u(k);for i=na:-1:2yk(i)=yk(i-1);endyk(1)=y(k);end%画参数随采样时刻的估计结果plot([1:N],thetae);xlabel('k');ylabel('参数估计a,b');legend('a1','a2','b0','b1');axis([0 N -2 2]);运行结果：图1 递推最小二乘法参数估计图中可以看出，辨识过程很不稳定，参数波动较大，因为在试验中将输入和噪声幅值一样大，噪声干扰相对输入较大，所以应该调整输入和噪声的幅值分别另：u=25*randn(N,1),50*randn(N,1),100*randn(N,1)所得结果如下图1）u=25*randn(N,1)2)u=50*randn(N,1)3)u=100*randn(N,1)综上可以看出只有当输入幅值较大时系统辨识的效果会更好，参数波动也不是很大。

3、遗忘因子最小二乘法clear all;close all;a=[1 -1.3 0.3]';b=[1 0.5]';d=1;%对象参数na=length(a)-1;nb=length(b)-1;%na,nb为A,B阶次N=200;%观测数据组uk=zeros(d+nb,1);%创建d+nb行1列的零列向量，给输入赋初值0yk=zeros(na,1);%输出初值u=1*randn(N,1);%输入采用白噪声序列xi=randn(N,1);%产生均值为零方差为1的白噪声序列w(k)thetae_1=zeros(na+nb+1,1);%thetae的chuzhiP=10^6*eye(na+nb+1);%P初值lambda=0.95%遗忘因子for k=1:Ntheta(:,k)=[a(2:na+1);b];%对象参数真值phi=[-yk;uk(d:d+nb)];%phi为行向量，组成phi矩阵y(k)=phi'*theta(:,k)+xi(k);%采集输出数据K=P*phi/(lambda+phi'*P*phi);thetae(:,k)=thetae_1+K*(y(k)-phi'*thetae_1);P=(eye(na+nb+1)-K*phi')*P/lambda;thetae_1=thetae(:,k);for i=d+nb:-1:2uk(i)=uk(i-1);enduk(1)=u(k);for i=na:-1:2yk(i)=yk(i-1);endyk(1)=y(k);end%画参数随采样时刻的估计结果subplot(1,2,1);plot([1:N],thetae(1:na,:));hold on;plot([1:N],theta(1:na,:),'k:');xlabel('k');ylabel('参数估计a');legend('a1','a2');axis([0 N -2 2]);subplot(1,2,2)plot([1:N],thetae(na+1:na+nb+1,:));hold on;plot([1:N],theta(na+1:na+nb+1,:),'k:'); xlabel('k');ylabel('参数估计b');legend('b0','b1');axis([0 N -0.5 2]);运行结果：改变输入信号的幅值，使u=12.5*randn(N,1),25*randn(N,1),50*randn(N,1),100*randn(N,1) 1)u=12.5*randn(N,1)2)u=25*randn(N,1)3)u=50*randn(N,1)4)4)u=100*randn(N,1)比较可知：只有当输入幅度较大时系统辨识结果更好【3-2】 设有被控过程：)()2.11()()6.07.11(1221k u z z k y z z ----+=+-给定期望传递函数的分母多项式为)08.06.01()(211---+-=z z z A m ，试按照极点配置方法设计控制系统，使期望输出无稳态误差，并写出控制表达式u(k)。