sitrun.m

基于MATLAB系统辨识工具箱的离线辨识程序

function  [sys0,sys1]=SITrun(mode)
clc;
ts=0.001;                                           %数据采样周期
N=1000;                                             %数据采样点数
if mode==0
    load SIdata.mat                                 %导入用于模型辨识的数据
    servo=iddata(z,u,ts);
    servo.InputName='position';
    servo.OutputName='distance';
    ze=servo(1:N);
    %figure, plot(ze(1:1000))
    %ze=dtrend(ze);   
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %ARX模型辨识及最优模型选择
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    nns=[1 1 0];nnt=[10 10 0];                      %设置ARX模型结构参数的范围
    %获取模型结构参数范围
    nas=nns(1);nat=nnt(1);
    nbs=nns(2);nbt=nnt(2);
    nk=nns(3);
    if (nat-nas+1)*(nbt-nbs+1)>150
        return 
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %ARX模型检验及最优模型选择
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    temp=0;
    for k1=nas:nat
        for k2=nbs:nbt
            m11=arx(ze,[k1 k2 nk]);                 %采用基于最小二乘法估计的ARX模型辨识
            [yh,fit]=compare(ze,m11);               %计算不同结构模型的的适应度
            if fit>temp
               temp=fit;m1=m11;                     %选择适应度最高的模型
               nn=[k1 k2 nk];
            end
            clc;
        end
    end
    %figure, compare(ze,m);                         %模型检验
    [yh1,fit1]=compare(ze,m1);
    save SIdata.mat 
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %ARMAX模型辨识及最优模型选择
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    disp('按任一键继续！')
    pause
    nns=[6 6 6 0];nnt=[10 10 10 0];                 %设置ARMAX模型结构参数的范围  
    %获取模型结构参数范围
    nas=nns(1);nat=nnt(1);
    nbs=nns(2);nbt=nnt(2);
    ncs=nns(3);nct=nnt(3);
    nk=nns(4);
    if (nat-nas+1)*(nbt-nbs+1)*(nct-ncs+1)>150
        return 
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %ARMAX模型检验及最优模型选择
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    temp=0;
    for k1=nas:nat
        for k2=nbs:nbt
            for k3=ncs:nct
                m22=armax(ze,[k1 k2 k3 nk]);        %采用基于最小二乘法估计的ARMAX模型辨识
                [yh,fit]=compare(ze,m22);           %计算不同结构模型的的适应度
                if fit>temp
                    temp=fit;m2=m22;                %选择适应度最高的模型
                    nn=[k1 k2 k3 nk];
                end
                clc;
            end
        end
    end
    %figure, compare(ze,m);                         %模型检验
    [yh2,fit2]=compare(ze,m2);
    save SIdata.mat                                 %模型及吻合程度存盘
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %BJ模型辨识及最优模型选择
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    disp('按任一键继续！')
    pause
    nns=[7 8 8 7 0];nnt=[10 10 10 10 0];            %设置BJ模型结构参数的范围
    %获取模型结构参数范围
    nbs=nns(1);nbt=nnt(1);
    ncs=nns(2);nct=nnt(2);
    nds=nns(3);ndt=nnt(3);
    nfs=nns(4);nft=nnt(4);
    nk=nns(5);
    if (nbt-nbs+1)*(nct-ncs+1)*(ndt-nds+1)*(nft-nfs+1)>150
            return 
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %BJ模型检验及最优模型选择
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    temp=0;
    for k1=nbs:nbt
        for k2=ncs:nct
            for k3=nds:ndt
                for k4=nfs:nft
                    m33=bj(ze,[k1 k2 k3 k4 nk]);    %采用基于最小二乘法估计的ARMAX模型辨识
                    [yh,fit]=compare(ze,m33);       %计算不同结构模型的的适应度
                    if fit>temp
                        temp=fit;m3=m33;            %选择适应度最高的模型
                        nn=[k1 k2 k3 k4 nk];
                    end
                    clc;
                end
            end
        end
    end
    %figure, compare(ze,m);                         %模型检验
    [yh3,fit3]=compare(ze,m3);
    save SIdata.mat                                 %模型及吻合程度存盘
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%最优模型及降阶模型分析
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load SIdata.mat;                                    %导入已辨识的模型和吻合程度及检验用的数据                             
svo=iddata(z,u,ts);svo.InputName='position';
svo.OutputName='distance';ze=svo(N+1:N*2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%最优模型选择，并求取其传递函数
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fit=max([fit1 fit2 fit3]);                          %选择吻合程度最高的模型
switch fit                         
    case fit1
        m=m1;disp('最优模型为ARX模型')
    case fit2
        m=m2;disp('最优模型为ARMAX模型')
    case fit3
        m=m3;disp('最优模型为Box-Jenkins模型')
end
%求最优模型的传递函数
my=noisecnv(m);
s=tf(m);                                            %将模型转化为传递函数
sys0=s(1);                                          %最优模型的传递函数
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%最优模型降阶，并转化为传递函数
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
mj=idmodred(m);                                     %模型降阶
mjj=noisecnv(mj);
sy=tf(mjj);                                         %降阶模型转化为传递函数
sys1=sy(1);                                         %降阶模型的传递函数

disp('降阶后模型的传递函数：')
sys1
[num,den]=tfdata(sys1);
disp('传递函数的系数向量：')
num=num{1}
den=den{1}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%最优模型及降阶后的模型与实际系统的吻合程度分析
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('检验数据下最优模型与实际系统的吻合程度：')
[yh1,f1]=compare(ze,m);
fit=f1   
disp('按任意键继续！')
pause

disp('检验数据下降阶模型与实际系统的吻合程度：')
[yh1,f1]=compare(ze,mj);
fit=f1  
figure,compare(ze,m,mj);