⭐ 欢迎来到虫虫下载站！ | 📦 资源下载 📁 资源专辑 ℹ️ 关于我们
⭐ 虫虫下载站
📤 上传资源

📄 work.m

📁 matlab系统辨识实现
💻 M
📖 第 1 页 / 共 2 页
字号:
🔍
12 
                    

% --- Executes on button press in pushbutton6.function pushbutton6_Callback(hObject, eventdata, handles)
%%%%%%%%%%%%%%%这些是辨识系统所需要的信号因为需要观察时刻前的数值，故往前再取50个点%%%%%%%%%%%%%
A=[0 1.0 0.5];
B=[1 -1.5 0.7];
C=[1 0 0];
D=B;
global RandomSignal;
global SampleNum;
global scale;
global Gauss;
 
    ww=zeros(1,SampleNum+50);
    ww=scale*Gauss;
 
    ee=zeros(1,SampleNum+50);
global ModelType;
    ModelType=get(handles.popupmenu3,'Value');
   if (ModelType==1)
       ee=dlsim(C,D,ww);
   end
   if (ModelType==2)
       ee=ww;
   end
 
    yy=zeros(1,SampleNum+50);
    yy=dlsim(A,B,RandomSignal);    
 
    zz=zeros(1,SampleNum+50);
    for i=1:(SampleNum+50)
    zz(i)=yy(i)+ee(i);
    end
uu=zeros(1,SampleNum+50);  
uu=RandomSignal;
%%%%%完毕%%%%%%%%这些是辨识系统所需要的信号因为需要观察时刻前的数值，故往前再取50个点%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%递推最小二乘算法计算过程%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
X=zeros(SampleNum,5);
Z=zeros(SampleNum,1);
for i=1:SampleNum
    X(i,1)=-zz(i+49);
    X(i,2)=-zz(i+48);
    X(i,3)=uu(i+50);
    X(i,4)=uu(i+49);
    X(i,5)=uu(i+48);
    Z(i,1)=zz(i+50);
end

sita=zeros(5,1);
mv=zeros(5,1);    
counter=0;      
a1=zeros(1,SampleNum);
a2=zeros(1,SampleNum);
b0=zeros(1,SampleNum);
b1=zeros(1,SampleNum);
b2=zeros(1,SampleNum);
ed=zeros(1,SampleNum);
%
P=zeros(5,5);
gama=0;
E=eye(5);
arfa=10^6;
P=(arfa^2).*E;
epsitong=0.0001;
ed=zeros(1,SampleNum);
for i=1:SampleNum
    counter=i;
    a1(i)=sita(1,1);
    a2(i)=sita(2,1);
    b0(i)=sita(3,1);
    b1(i)=sita(4,1);
    b2(i)=sita(5,1);
    gama=1/(1+X(i,:)*P*(X(i,:))');
    sita=mv+P*(X(i,:))'*gama*(Z(i,1)-X(i,:)*mv);
    ed(i)=Z(i,1)-X(i,:)*mv;
    P=P-P*(X(i,:))'*gama*X(i,:)*P;
    if (mv~=0)&(abs((sita-mv)./mv)<epsitong);
        break;
    end
    mv=sita;
end
%%%%%%%%%%%%%%%%这些是辨识系统所需要的信号因为需要观察时刻前的数值，故往前再取50个点%%%%%%%%%%%%%

%取预报误差植ed(i)
for i=1:counter
    ed1(i)=ed(i);
end

%%用所得到结果的检验误差
A1=[0,sita(1,1),sita(2,1)];B1=[sita(3,1),sita(4,1),sita(5,1)];
lt=round(counter/2);%取后一半数据
uu1=zeros(1,counter-lt); 
for i=(lt+1):counter
    uu1(i-lt)=uu(50+i);
end
yy1=dlsim(A,B,uu1);%根据预测的模型判断输出 
for i=(lt+1):counter
    jy(i)=zz(50+i)-yy1(i-lt);%计算误差检验
end
axes(handles.axes1);

t1=1:1:counter;%输出函数
plot(t1,ed1,'r',t1,jy,'b');grid;  title('预报误差曲线(可用鼠标选定制定区域放大)');legend('\fontsize{8}预报误差','参数检验');zoom on;hold off;


% --- Executes on button press in pushbutton7.function pushbutton7_Callback(hObject, eventdata, handles)%%%%%%%%%%%%%%%%这些是辨识系统所需要的信号因为需要观察时刻前的数值，故往前再取50个点%%%%%%%%%%%%%
A=[0 1.0 0.5];
B=[1 -1.5 0.7];
C=[1 0 0];
D=B;
global RandomSignal;
global SampleNum;
global scale;
global Gauss;
 
    ww=zeros(1,SampleNum+50);
    ww=scale*Gauss;
 
    ee=zeros(1,SampleNum+50);
global ModelType;
    ModelType=get(handles.popupmenu3,'Value');
   if (ModelType==1)
       ee=dlsim(C,D,ww);
   end
   if (ModelType==2)
       ee=ww;
   end
 
    yy=zeros(1,SampleNum+50);
    yy=dlsim(A,B,RandomSignal);    
 
    zz=zeros(1,SampleNum+50);
    for i=1:(SampleNum+50)
    zz(i)=yy(i)+ee(i);
    end
   
 
uu=zeros(1,SampleNum+50);  
uu=RandomSignal;
%%%%%完毕%%%%%%%%这些是辨识系统所需要的信号因为需要观察时刻前的数值，故往前再取50个点%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%递推最小二乘算法计算过程%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
X=zeros(SampleNum,5);
Z=zeros(SampleNum,1);
for i=1:SampleNum
    X(i,1)=-zz(i+49);
    X(i,2)=-zz(i+48);
    X(i,3)=uu(i+50);
    X(i,4)=uu(i+49);
    X(i,5)=uu(i+48);
    Z(i,1)=zz(i+50);
end

sita=zeros(5,1);
mv=zeros(5,1);    
counter=0;      
a1=zeros(1,SampleNum);
a2=zeros(1,SampleNum);
b0=zeros(1,SampleNum);
b1=zeros(1,SampleNum);
b2=zeros(1,SampleNum);
%
P=zeros(5,5);
gama=0;
E=eye(5);
arfa=10^6;
P=(arfa^2).*E;
epsitong=0.0001;

for i=1:SampleNum
    counter=i;
    a1(i)=sita(1,1);
    a2(i)=sita(2,1);
    b0(i)=sita(3,1);
    b1(i)=sita(4,1);
    b2(i)=sita(5,1);
    gama=1/(1+X(i,:)*P*(X(i,:))');
    sita=mv+P*(X(i,:))'*gama*(Z(i,1)-X(i,:)*mv);
    P=P-P*(X(i,:))'*gama*X(i,:)*P;
    if (mv~=0)&(abs((sita-mv)./mv)<epsitong);
        break;
    end
    mv=sita;
end
%%%%%%完毕%%%%%%%%%%%%%%%%%%%%%递推最小二乘算法计算过程%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%显示数据收敛过程　
viewa1=zeros(1,counter);
viewa2=zeros(1,counter);
viewb0=zeros(1,counter);
viewb1=zeros(1,counter);
viewb2=zeros(1,counter);
for i=1:counter
    viewa1(i)=a1(i);
    viewa2(i)=a2(i);
    viewb0(i)=b0(i);
    viewb1(i)=b1(i);
    viewb2(i)=b2(i);
end
axes(handles.axes2);
t=1:1:counter;        %输出函数

plot(t,viewa1,'b',t,viewa2,'g',t,viewb0,'r',t,viewb1,'m',t,viewb2,'k'); 
grid;title('递推最小二乘法数据收敛过程(可用鼠标选定制定区域放大观察)');legend('\fontsize{5}a1','a2','b0','b1','b2',-1);
zoom on;% --- Executes on button press in pushbutton8.function pushbutton8_Callback(hObject, eventdata, handles)%%%%%%%%%%%%%%%%这些是辨识系统所需要的信号因为需要观察时刻前的数值，故往前再取50个点%%%%%%%%%%%%%
A=[0 1.0 0.5];
B=[1 -1.5 0.7];
C=[1 0 0];
D=B;
global RandomSignal;
global SampleNum;
global scale;
global Gauss;
 
    ww=zeros(1,SampleNum+50);
    ww=scale*Gauss;
 
    ee=zeros(1,SampleNum+50);
global ModelType;
    ModelType=get(handles.popupmenu3,'Value');
   if (ModelType==1)
       ee=dlsim(C,D,ww);
   end
   if (ModelType==2)
       ee=ww;
   end
 
    yy=zeros(1,SampleNum+50);
    yy=dlsim(A,B,RandomSignal);    
 
    zz=zeros(1,SampleNum+50);
    for i=1:(SampleNum+50)
    zz(i)=yy(i)+ee(i);
    end

uu=zeros(1,SampleNum+50);  
uu=RandomSignal;
%%%%%完毕%%%%%%%%这些是辨识系统所需要的信号因为需要观察时刻前的数值，故往前再取50个点%%%%%%%%%%%%%


%用整批算法获得一初始数据
X=zeros(SampleNum,5);
Z=zeros(SampleNum,1);
for i=1:SampleNum
X(i,1)=-zz(i+47);
X(i,2)=-zz(i+46);
X(i,3)=uu(i+48);
X(i,4)=uu(i+47);
X(i,5)=uu(i+46);
Z(i,1)=zz(i+48);
end
sita=zeros(5,1);
sita=(inv(X'*X))*X'*Z;


mv=zeros(5,1);%引入中间变量middlevariable，整批算法获得的初始数据


E=zeros(SampleNum,2);
ee=zeros(SampleNum+2,1);
eee=zeros(SampleNum,1);
f=zeros(2,1);
zzz=zeros(1,SampleNum+50);%中间变量
uuu=zeros(1,SampleNum+50);%中间变量
counter=0;%计数器
epsitong=0.000000001;%循环停机条件　
circletime=10000;%设定最大循环次数
a1=zeros(1,circletime);
a2=zeros(1,circletime);
b0=zeros(1,circletime);
b1=zeros(1,circletime);
b2=zeros(1,circletime);



%%%%%%%广义最小二乘法计算过程%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%％％％
for i=1:circletime%进入循环
    counter=i;
    a1(i)=sita(1,1);
    a2(i)=sita(2,1);
    b0(i)=sita(3,1);
    b1(i)=sita(4,1);
    b2(i)=sita(5,1);
    %模型残差　
    for k=1:SampleNum+2
        ee(k,1)=zz(k+48)+sita(1,1)*zz(k+47)+sita(2,1)*zz(k+46)-sita(3,1)*uu(k+48)-sita(4,1)*uu(k+47)-sita(5,1)*uu(k+46);
    end
    %对误差序列e(k)进行估计％％％％％％％％％％％％％％％％
    for j=1:SampleNum
        E(j,1)=-ee(j+1,1);
        E(j,2)=-ee(j,1);
        eee(j,1)=ee(j+2,1);
    end
    f=inv(E'*E)*E'*eee;
    %完毕%%对误差序列e(k)进行估计％％％％％％％％％％％％％％％％
    
    %%对输入输出观测序列滤波%%%%%%%%%%%%%%
    for i=49:SampleNum+50
        zzz(i)=zz(i)+f(1,1)*zz(i-1)+f(2,1)*zz(i-2);
        uuu(i)=uu(i)+f(1,1)*uu(i-1)+f(2,1)*uu(i-2);
    end
      
    for i=49:SampleNum+50
        zz(i)=zzz(i);
        uu(i)=uuu(i);
    end
    
    
    for i=1:SampleNum
        X(i,1)=-zz(i+47);
        X(i,2)=-zz(i+46);
        X(i,3)=uu(i+48);
        X(i,4)=uu(i+47);
        X(i,5)=uu(i+46);
        Z(i,1)=zz(i+48);
    end
   sita=zeros(5,1);
   sita=(inv(X'*X))*X'*Z;
    if (mv~=0)&(abs((sita-mv)./mv)<epsitong);%满足停机条件
        break;
    end
    mv=sita;%%得到结果
end
   s=(sita-[-1.5 0.7 0 1 0.5]')'*(sita-[-1.5 0.7 0 1 0.5]');
%%%完毕%广义最小二乘法计算过程%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%输出
    set(handles.text77,'string',sita(1,1));
    set(handles.text80,'string',sita(2,1));
    set(handles.text78,'string',sita(3,1));
    set(handles.text81,'string',sita(4,1));
    set(handles.text79,'string',sita(5,1));
    set(handles.text85,'string', counter);
    set(handles.text94,'string',s);
 % --- Executes on button press in pushbutton9.function pushbutton9_Callback(hObject, eventdata, handles)%%%%%%%%%%%%%%%%这些是辨识系统所需要的信号因为需要观察时刻前的数值，故往前再取50个点%%%%%%%%%%%%%
A=[0 1.0 0.5];
B=[1 -1.5 0.7];
C=[1 0 0];
D=B;
global RandomSignal;
global SampleNum;
global scale;
global Gauss;
 
    ww=zeros(1,SampleNum+50);
    ww=scale*Gauss;
 
    ee=zeros(1,SampleNum+50);
global ModelType;
    ModelType=get(handles.popupmenu3,'Value');
   if (ModelType==1)
       ee=dlsim(C,D,ww);
   end
   if (ModelType==2)
       ee=ww;
   end
 
    yy=zeros(1,SampleNum+50);
    yy=dlsim(A,B,RandomSignal);    
 
    zz=zeros(1,SampleNum+50);
    for i=1:(SampleNum+50)
    zz(i)=yy(i)+ee(i);
    end

uu=zeros(1,SampleNum+50);  
uu=RandomSignal;
%%%%%完毕%%%%%%%%这些是辨识系统所需要的信号因为需要观察时刻前的数值，故往前再取50个点%%%%%%%%%%%%%

%%%%%%%广义最小二乘法计算过程%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%％％％
X=zeros(SampleNum,5);
Z=zeros(SampleNum,1);
for i=1:SampleNum
X(i,1)=-zz(i+47);
X(i,2)=-zz(i+46);
X(i,3)=uu(i+48);
X(i,4)=uu(i+47);
X(i,5)=uu(i+46);
Z(i,1)=zz(i+48);
end
sita=zeros(5,1);
sita=(inv(X'*X))*X'*Z;


mv=zeros(5,1);   %引入中间变量middlevariable

E=zeros(SampleNum,2);
ee=zeros(SampleNum+2,1);
eee=zeros(SampleNum,1);
f=zeros(2,1);
zzz=zeros(1,SampleNum+50);%中间变量
uuu=zeros(1,SampleNum+50);%中间变量
counter=0;%计数器
epsitong=0.000000001;
circletime=10000;
a1=zeros(1,circletime);
a2=zeros(1,circletime);
b0=zeros(1,circletime);
b1=zeros(1,circletime);
b2=zeros(1,circletime);

for i=1:circletime %进入循环
    counter=i;
    a1(i)=sita(1,1);
    a2(i)=sita(2,1);
    b0(i)=sita(3,1);
    b1(i)=sita(4,1);
    b2(i)=sita(5,1);
    for k=1:SampleNum+2
        ee(k,1)=zz(k+48)+sita(1,1)*zz(k+47)+sita(2,1)*zz(k+46)-sita(3,1)*uu(k+48)-sita(4,1)*uu(k+47)-sita(5,1)*uu(k+46);
    end
    for j=1:SampleNum
        E(j,1)=-ee(j+1,1);
        E(j,2)=-ee(j,1);
        eee(j,1)=ee(j+2,1);
    end
  f=inv(E'*E)*E'*eee;
    for i=49:SampleNum+50
        zzz(i)=zz(i)+f(1,1)*zz(i-1)+f(2,1)*zz(i-2);
        uuu(i)=uu(i)+f(1,1)*uu(i-1)+f(2,1)*uu(i-2);
    end
    for i=49:SampleNum+50
        zz(i)=zzz(i);
        uu(i)=uuu(i);
    end
    
    
    for i=1:SampleNum
        X(i,1)=-zz(i+47);
        X(i,2)=-zz(i+46);
        X(i,3)=uu(i+48);
        X(i,4)=uu(i+47);
        X(i,5)=uu(i+46);
        Z(i,1)=zz(i+48);
    end
   sita=zeros(5,1);
   sita=(inv(X'*X))*X'*Z;

    if (mv~=0)&(abs((sita-mv)./mv)<epsitong);
        break;
    end
    mv=sita;
end
%%%完毕%广义最小二乘法计算过程%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%       


%%%%数据收敛过程
viewa1=zeros(1,counter);
viewa2=zeros(1,counter);
viewb0=zeros(1,counter);
viewb1=zeros(1,counter);
viewb2=zeros(1,counter);
for i=1:counter
    viewa1(i)=a1(i);
    viewa2(i)=a2(i);
    viewb0(i)=b0(i);
    viewb1(i)=b1(i);
    viewb2(i)=b2(i);
end

axes(handles.axes2);
t=1:1:counter;        %输出函数
plot(t,viewa1,'b',t,viewa2,'g',t,viewb0,'r',t,viewb1,'m',t,viewb2,'k');hold on;
plot(t,viewa1,'b+',t,viewa2,'g+',t,viewb0,'r+',t,viewb1,'m+',t,viewb2,'k+');zoom on;title('广义最小二乘法数据收敛过程(可用鼠标选定制定区域放大观察)'); legend('\fontsize{5}a1','a2','b0','b1','b2',-1);
hold off;
function pushbutton11_Callback(hObject, eventdata, handles)
msgbox('采用直接给定递推算法的初始值，系统初始参数设置为0，采用停机准则一，当两次迭带参数差值小于0.00001则停止迭带','递推最小二乘法参数设置') ;function pushbutton12_Callback(hObject, eventdata, handles)
msgbox('广义算法开始迭带前已经有整批算法给定一个初值，当两次迭带参数差值小于0.000000001则停止迭带','广义最小二乘法参数设置') ;
12

⌨️ 快捷键说明

复制代码 Ctrl + C
搜索代码 Ctrl + F
全屏模式 F11
切换主题 Ctrl + Shift + D
显示快捷键 ?
增大字号 Ctrl + =
减小字号 Ctrl + -