文档介绍：该【过程辨识程序集合 】是由【zkusha】上传分享，文档一共【19】页，该文档可以免费在线阅读，需要了解更多关于【过程辨识程序集合 】的内容，可以使用淘豆网的站内搜索功能，选择自己适合的文档，以下文字是截取该文章内的部分文字，如需要获得完整电子版，请下载此文档到您的设备，方便您编辑和打印。1
%% 1统计近似抽样法 %%
clear
clc
%***生成（0,1）均匀分布的随机序列***%
%-------------初始化-------------%
A=179;
M=2^35;
x(1)=11;
n=111;
%-------生成均匀分布随机序列------%
for i=1:n
x(i+1)=mod(A*x(i),M);
end
x1=x/M;
for i=2:(n+1)
u(i-1)=x1(i);
end
plot(u);
%--计算生成的均匀分布的均值和方差--%
sum=u(1);
for i=2:111
sum=sum+u(i);
end
e=sum/111
sum=(u(1)-e)^2;
for i=2:111
sum=sum+(u(i)-e)^2;
end
d2=sum/111
d1=sqrt(d2)
%****生成标准正态分布***%
%-------------初始化-------------%
miu=0;
sigma=1;
%-------生成正态分布随机序列------%
for i=1:100
sum=u(i);
for m=(i+1):(i+11)
sum=sum+u(m);
end
elta(i)=miu+sigma*(sum-6);
end
figure(2)
plot(elta);
sum=elta(1);
%--计算生成的正态分布的均值和方差--%
for i=2:100
sum=sum+elta(i);
end
E=sum/100
sum=(elta(1)-E)^2;
for i=2:100
sum=sum+(elta(i)-E)^2;
end
D2=sum/100
D1=sqrt(D2)
%2变换抽样法
clear
clc
%用混合同余法生成（0,1）均匀分布的随机序列--%
M=2^35;
c1=1;
c2=*(2^35);
A1=2^12+1;
A2=2^7+1;
x1(1)=20867350019;
x2(1)=**********;
for i=1:100
x1(i+1)=mod(x1(i)*A1+c1,M);
x2(i+1)=mod((x2(i)*A2+c2)*A2+c2,M);
end
xi1=x1/M;
xi2=x2/M;
subplot(2,1,1);
plot(xi1);
subplot(2,1,2);
plot(xi2);
%---------用变换抽样法生成两个相互独立的标准正态分布---------%
for i=2:101
elta1(i-1)=(-2*log(xi1(i)))^(1/2)*cos(2*pi*xi2(i));
elta2(i-1)=(-2*log(xi1(i)))^(1/2)*sin(2*pi*xi2(i));
end
figure(2);
subplot(2,1,1);
plot(elta1);
19
subplot(2,1,2);
plot(elta2);
%-------------求生成的两个标准正态分布的均值----------------%
sum1=elta1(1);
sum2=elta2(1);
for i=2:100
sum1=sum1+elta1(i);
sum2=sum2+elta2(i);
end
E1=sum1/100
E2=sum2/100
%-------------求生成的两个标准正态分布的方差----------------%
sum1=(elta1(1)-E1)^2;
sum2=(elta2(1)-E2)^2;
for i=2:100
sum1=sum1+(elta1(i)-E1)^2;
sum2=sum2+(elta2(i)-E2)^2;
end
D1(1)=sum1/100;
D1(2)=sqrt(D1(1));
D2(1)=sum2/100;
D2(2)=sqrt(D2(1));
D1
D2
%--M序列，逆M序列--%
clear
clc
%-------------初始化------------%
M(1)=1;M(2)=0;M(3)=0;
M(4)=1;M(5)=1;M(6)=0;
M(7)=0;M(8)=0;M(9)=1;
P=9;
n=2^9-1;
%-------产生一个周期的M序列------%
for i=10:n
M(i)=xor(M(i-4),M(i-9));
end
figure(1);
subplot(2,1,1);
plot(M)
%---统计产生的M序列中0、1的个数---%
num1=0;
num0=0;
for i=1:n
if M(i)==1
num1=num1+1;
end
if M(i)==0
num0=num0+1;
end
end
fprintf('逻辑0出现的次数为：%\n',num0);
fprintf('逻辑1出现的次数为：%\n',num1);
%-------统计游程-------%
you1=0;
you2=0;
you3=0;
you4=0;
you5=0;
you6=0;
you7=0;
you8=0;
you9=0;
you=1;
sum_you_0=0;
sum_you_1=0;
M(512)=2; %为了循环判断能够到达第511项故假设x(512)=2,之所以设其值为2是要其与x(511)项必须不同否则影响游程的判断；
for i=1:n
%i
a=M(i)-M(i+1);
if a==0
you=you+1;
continue;
else
if you==1
you1=you1+1;
if M(i)==0
sum_you_0=sum_you_0+1;
else
3
sum_you_1=sum_you_1+1;
end
end
if you==2
you2=you2+1;
if M(i)==0
sum_you_0=sum_you_0+1;
else
sum_you_1=sum_you_1+1;
end
end
if you==3
you3=you3+1;
if M(i)==0
sum_you_0=sum_you_0+1;
else
sum_you_1=sum_you_1+1;
end
end
if you==4
you4=you4+1;
if M(i)==0
sum_you_0=sum_you_0+1;
else
sum_you_1=sum_you_1+1;
end
end
if you==5
you5=you5+1;
if M(i)==0
sum_you_0=sum_you_0+1;
else
sum_you_1=sum_you_1+1;
end
end
if you==6
you6=you6+1;
if M(i)==0
sum_you_0=sum_you_0+1;
else
sum_you_1=sum_you_1+1;
end
end
if you==7
you7=you7+1;
if M(i)==0
sum_you_0=sum_you_0+1;
else
sum_you_1=sum_you_1+1;
end
end
if you==8
you8=you8+1;
if M(i)==0
sum_you_0=sum_you_0+1;
disp('长度为P-1即8的游程为“0”游程');
else
sum_you_1=sum_you_1+1;
disp('长度为P-1即8的游程为“1”游程');
end
end
if you==9
you9=you9+1;
if M(i)==0
sum_you_0=sum_you_0+1;
disp('长度为P即8的游程为“0”游程');
else
sum_you_1=sum_you_1+1;
4
disp('长度为P即9的游程为“1”游程');
end
end
you=1;
end
end
fprintf('“0”游程的总数为：%\n',sum_you_0);
fprintf('“1”游程的总数为：%\n',sum_you_1);
fprintf('游程数为1的游程个数为：%\n',you1);
fprintf('游程数为2的游程个数为：%\n',you2);
fprintf('游程数为3的游程个数为：%\n',you3);
fprintf('游程数为4的游程个数为：%\n',you4);
fprintf('游程数为5的游程个数为：%\n',you5);
fprintf('游程数为6的游程个数为：%\n',you6);
fprintf('游程数为7的游程个数为：%\n',you7);
fprintf('游程数为8的游程个数为：%\n',you8);
fprintf('游程数为9的游程个数为：%\n',you9);
%-----------逆M序列-----------%
s(1)=1;
for i=1:(n-1)
if s(i)==1
s(i+1)=0;
else
s(i+1)=1;
end
end
subplot(2,1,2);
plot(s);
for i=1:n
nM(i)=xor(M(i),s(i));
end
figure(2)
plot(nM);
%--------将逆M序列的幅值变为-A，A-------%
A=2;
for i=1:n
if nM(i)==0
nM_A(i)=-A;
else
nM_A(i)=A;
end
end
figure(3)
plot(nM_A);
% 基本最小二乘法一次计算程序 %
clear
clc
L=400;
%-------产生4阶M序列------%
M=[1 0 0 1];
for i=5:(L+3)
M(i)=xor(M(i-4),M(i-3));
end
u=M;
%------产生服从N（0,1）的白噪声------%
v=randn(1,L+3);
%------产生z------%
a1=-;
a2=;
b1=1;
b2=;
z=[1 0];
for i=3:(L+3)
z(i)=b1*u(i-1)+b2*u(i-2)+v(i)-a1*z(i-1)-a2*z(i-2);
end
%------产生HL------%
na=2;
nb=2;
for j=1:(na+nb)
for i=1:L
if j<=na
HL(i,j)=-z(i+na-j);
else
5
HL(i,j)=u(i+na+nb-j);
end
end
end
%-----得到参数估计值------%
theta=inv(HL'*HL)*HL'*(z(3:402))';
disp(' 基本最小二乘法一次完成算法的计算结果');
disp('参 数 a1 a2 b1 b2');
fprintf('真 值 % % % %\n',a1,a2,b1,b2);
fprintf('估计值 % % % %\n',theta(1),theta(2),theta(3),theta(4));
% 递推最小二乘法算法 %
clear
clc
L=480;
%------产生服从N（0,1）的白噪声------%
v=randn(1,L+3);
%----------产生4阶逆M序列-----------%
M=[1 1 0 1];
for i=5:(L+3)
M(i)=xor(M(i-4),M(i-1));
end
s(1)=1;
for i=1:(L+2)
if s(i)==1
s(i+1)=0;
else
s(i+1)=1;
end
end
for i=1:(L+3)
nM(i)=xor(M(i),s(i));
end
u=nM;
%------产生z------%
a1=-;
a2=;
b1=1;
b2=;
z=[1 0];
for i=3:(L+3)
z(i)=b1*u(i-1)+b2*u(i-2)+v(i)-a1*z(i-1)-a2*z(i-2);
end
%---------递推部分-------%
P=(10^6)*eye(4);
theta=*ones(4,1);
na=2;
nb=2;
a_1(1)=theta(1);
a_2(1)=theta(2);
b_1(1)=theta(3);
b_2(1)=theta(4);
for k=3:(L+2)
for i=1:(na+nb) %生成h向量
if i<=na
h0(i)=-z(k-i);
else
h0(i)=u(k+na-i);
end
end
%Theta_mid=Theta;
%P_mid=P;
h=h0';
K=P*h*inv(h'*P*h+1);
theta=theta+K*(z(k)-h'*theta);
P=(eye(4)-K*h')*P;
a_1(k-1)=theta(1); %存储每一次计算a1,a2,b1,b2的值
a_2(k-1)=theta(2);
b_1(k-1)=theta(3);
b_2(k-1)=theta(4);
end
disp(' 加权最小二乘法递推算法的辨识结果');
disp('参 数 a1 a2 b1 b2');
fprintf('真 值 % % % %\n'
7
,a1,a2,b1,b2);
fprintf('估计值 % % % %\n',theta(1),theta(2),theta(3),theta(4));
x=length(a_1);
i=1:x;
figure(1)
plot(i,a_1,i,a_2,i,b_1,i,b_2);
legend('a_1','a_2','b_1','b_2');
title('参数估计值变化过程');
%加权最小二乘法一次计算程序
clear
clc
L=400;
%-------产生4阶M序列------%
M=[1 0 0 1];
for i=5:(L+3)
M(i)=xor(M(i-4),M(i-3));
end
u=M;
%------产生服从N（0,1）的白噪声------%
v=randn(1,L+3);
%------产生z------%
a1=-;
a2=;
b1=1;
b2=;
z=[1 0];
for i=3:(L+3)
z(i)=b1*u(i-1)+b2*u(i-2)+v(i)-a1*z(i-1)-a2*z(i-2);
end
%------产生HL------%
na=2;
nb=2;
for j=1:(na+nb)
for i=1:L
if j<=na
HL(i,j)=-z(i+na-j);
else
HL(i,j)=u(i+na+nb-j);
end
end
end
%------产生加权矩阵-------%
miu=;
for i=1:L
for j=1:L
if i==j
I(i,j)=miu^(L-i);
else
I(i,j)=0;
end
end
end
I(i)=miu^(L-i);
%-----得到参数估计值------%
theta=inv(HL'*I*HL)*HL'*I*(z(3:402))';
disp(' 加权最小二乘法一次完成算法的计算结果');
disp('参 数 a1 a2 b1 b2');
fprintf('真 值 % % % %\n',a1,a2,b1,b2);
fprintf('估计值 % % % %\n',theta(1),theta(2),theta(3),theta(4));
% 广义最小二乘法 %
clear
clc
L=483;
%------产生e------%
v=randn(1,L);
e(1)=v(1);
e(2)=v(2);
c1=0;
c2=0;
for i=3:L
e(i)=c1*e(i-1)+c2*e(i-2)+v(i);
end
%----产生4阶逆M序列
7
M=[1 1 0 1];
for i=5:L
M(i)=xor(M(i-4),M(i-1));
end
s(1)=1;
for i=1:L
if s(i)==1
s(i+1)=0;
else
s(i+1)=1;
end
end
for i=1:L
nM(i)=xor(M(i),s(i));
end
u=nM;
%------产生z------%
a1=-;
a2=;
b1=1;
b2=;
z=[1 0];
for i=3:L
z(i)=b1*u(i-1)+b2*u(i-2)+e(i)-a1*z(i-1)-a2*z(i-2);
end
%------产生zf------%
zf(1)=1;
zf(2)=0;
for i=3:L
zf(i)=z(i)-0*z(i-1)-0*z(i-2);
end
%------产生uf------%
uf(1)=u(1);
uf(2)=u(2);
for i=3:L
uf(i)=u(i)-0*u(i-1)-0*u(i-2);
end
%-----递推部分初始化----%
P=(10^6)*eye(4);
theta=*ones(4,1);
PE=10*eye(2);
thetaE=[;];
na=2;
nb=2;
a_1(1)=theta(1);
a_2(1)=theta(2);
b_1(1)=theta(3);
b_2(1)=theta(4);
for i=3:L-3
for k=1:(na+nb) %生成h向量
if k<=na
h0(k)=-zf(i-k);
else
h0(k)=uf(i+na-k);
end
end
h=h0';
K=P*h*inv(h'*P*h+1);
theta=theta+K*(z(i)-h'*theta);
P=(eye(4)-K*h')*P;
he=[-e(i-1);-e(i-2)];
KE=PE*he*inv(1+he'*PE*he);
thetaE=thetaE+KE*(e(i)-he'*thetaE);
PE=(eye(2)-KE*he')*PE;
a_1(i-1)=theta(1);
a_2(i-1)=theta(2);
b_1(i-1)=theta(3);
b_2(i-1)=theta(4);
c_1(i-1)=thetaE(1);
c_2(i-1)=thetaE(2);
end
%-------------%
disp(' 广义最小二乘法的辨识结果');
disp('参 数 a1 a2 b1 b2 c1 c2');
fprintf('真 值 % % % % % %\n',a1,a2,b1,b2,c1,c2);
fprintf('估计值 % % % % % %\n',theta(1),theta(2),theta(3),theta(4),thetaE(1),thetaE(2));
9
x=length(a_1);
i=1:x;
figure(1)
plot(i,a_1,i,a_2,i,b_1,i,b_2,i,c_1,i,c_2);
legend('a_1','a_2','b_1','b_2','c_1','c_2');
title('参数估计值变化过程');
% 偏差补偿最小二乘算法
clear
clc
L=483;
%------产生服从N（0,1）的白噪声------%
v=randn(1,L);
%----------产生4阶逆M序列-----------%
M=[1 1 0 1];
for i=5:L
M(i)=xor(M(i-4),M(i-1));
end
s(1)=1;
for i=1:L
if s(i)==1
s(i+1)=0;
else
s(i+1)=1;
end
end
for i=1:L
nM(i)=xor(M(i),s(i));
end
u=nM;
%------产生z------%
a1=-;
a2=;
b1=1;
b2=;
z=[1 0];
for i=3:(L)
z(i)=b1*u(i-1)+b2*u(i-2)+v(i)-a1*z(i-1)-a2*z(i-2);
end
%-----递推部分初始化----%
P=(10^6)*eye(4);
theta=*ones(4,1);
K=[10;10;10;10]; %增益
J=0;
thetaC=[2;3;1;];
D=[1 0 0 0;0 1 0 0;0 0 0 0;0 0 0 0];
na=2;
nb=2;
a_1(1)=theta(1);
a_2(1)=theta(2);
b_1(1)=theta(3);
b_2(1)=theta(4);
for i=3:(L-1)
for k=1:(na+nb) %生成h向量
if k<=na
h0(k)=-z(i-k);
else
h0(k)=u(i+na-k);
end
end
h=h0';
K=P*h*inv(h'*P*h+1);
theta=theta+K*(z(i)-h'*theta);
P=(eye(4)-K*h')*P;
J=J+(z(i-1)-h'*theta)^2/(1+h'*P*h);
g=J/((i-1)*(1+(thetaC'*D*theta)));
thetaC=theta+(i-1)*g*P*D*thetaC;
a_1(i-1)=theta(1); %存储每一次计算a1,a2,b1,b2的值
a_2(i-1)=theta(2);
b_1(i-1)=theta(3);
b_2(i-1)=theta(4);
end
%-------------%
disp(' 偏差补偿最小二乘法算法的