智能控制系统matlab仿真

智能控制系统实验报告ARMA模型ARMA(p,q)是一个线性时间序列预测模型，适用于平稳的时间序列，即对于任何时刻t，都有()tEZ,E(at)=0.协方差矩阵(')ttEaa,对于任意0l有(')0ttlEaa。AR模型11ttptptZZZ(0.1)当误差项t自相关时，可以被有限阶滑动平均表示11tttqtqaaa(0.2)这里ta是零均白噪声，协方差矩阵a非奇异。结合AR过程和MA误差项，得到ARMA过程：111111ttptpttptpttqtqZZZZZaaa(0.3)对于一个很大的阶数n，AR(n)接近ARMA(p,q)1()()nttitiiZnZan(0.4)从(0.4)得到残差的估计值1ˆˆ()()nttitiianZnZ(0.5)其中ˆ()in利用多变量最小二乘法求解，然后使用估计值ˆ()tan建立多变量回归模型1111ˆˆttptpttqtqZZZaaa(0.6)1111[,,:,.]()ˆ()ˆ()ttptpqtttqZZZananan(0.7)(1:)0[,]TZYA(0.8)01[,,]TAaa(0.9)000'ˆˆˆ()aAAnT(0.10)最小二乘法求解公式，以AR(p)为例。111121[,,,,]ttptpttpttpZZZaZaZ112111112111[,,][,,,]TTkTpppTpZZZZZZZZ(0.11)记11(1)101(1),:[,,]:[,,,],[1,',,']':[,,]1TkTTkTpkkptttpTkpTZ,,ZaavYZZYYZ=Y+AZ:=[]AY(0.12)()()kzIuvecY(0.13)最小化目标函数11()'[(')]'()[(')]ukTakSuuzIIzIYY(0.14)解得111ˆ'(')(+)'(')'(')ZYYYYAYYYAYYY(0.15)以下是仿真例题：例题1：使用AR模型预测德国西部经济数据。该数据是三元时间序列，包括投资(Investment)，可支配输入(income)和消费支出(consumption)三个变量。现使用AR模型预测INCOME，将1960年到1978年的75组用于训练，其余用于测试，利用最小二乘回归估计系数矩阵12[,,,,]pv代码如下：clearclcdata1=load('e1.txt');%data1=data1';data2=log(data1);data3=data2(:,2:end)-data2(:,1:end-1);%差分数据[K,num]=size(data3);data=data3;T=73;recordaic=100;forp=2:-1:2clearZ;clearZc;fort=p+1:T+pY(:,t-p)=data(:,t);%构造矩阵Ymiddle1=data(:,(t-p):(t-1));middle=fliplr(middle1);Z(:,t-p)=[1;middle(:)];%构造矩阵Zend%---------------直接求伪逆%B1=Y*pinv(Z)%------------------y=Y(:);middleb=inv(Z*Z')*Z;b=kron(middleb,eye(K))*y;%------------------最小二乘函数b1=regress(Y(:),kron(Z',eye(K)));B1=reshape(b,K,K*p+1)%权值矩阵%-------------B=reshape(b,K,K*p+1);%权值矩阵%B=[-0.017,-0.32,0.146,0.961,-0.161,0.115,0.934;0.016,0.044,-0.153,0.289,0.050,...%0.019,-0.010;0.013,-0.002,0.225,-0.264,0.034,0.355,-0.022];fort=p+1:num%fort=T+p+1:numO(:,t-p)=data(:,t);middle2=data(:,(t-p):(t-1));middlec=fliplr(middle2);Zc(:,t-p)=[1;middlec(:)];endYc=B*Zc;%----------直接伪逆Yc1=B1*Zc;data41=Yc1+data2(:,p+1:end-1);%第三个到第八十九个数dataf1=exp(data41);%反对数rmse1=(sum((data1(3,3:91)-dataf1(3,:)).^2)/(89-1))^0.5%计算均方根误差%-----------data4=Yc+data2(:,p+1:end-1);%%反差分第三个到第八十九个数dataf=exp(data4);%反对数U=dataf(:,1:T)-data1(:,p+1:T+p);sigmaU=(1/T)*U*U';aic=log(det(sigmaU))+2*p*K*K/T;ifaicrecordaicrecordaic=aic;m=p;Yc_b=Yc;dataf_b=dataf;endendsubplot(2,1,1)plot(1:92-1-m,Yc_b(3,:),'b',1:92-1-m,data3(3,m+1:91),'r');%差分数据画图subplot(2,1,2)plot(1:92-1-m,data1(3,m+1:91),'b',1:92-1-m,dataf_b(3,:),'r');%实际数据画图rmse=(sum((data1(3,m+1:91)-dataf(3,:)).^2)/(91-m-1))^0.5%计算均方根误差运行结果如下：实验分析：效果图如下图所示上面的子图是对数差分数据预测效果，下面的子图是反差分反对数的预测效果图，预测均方根误差为25.1096.可以看出，预测效果很好。例题2：课件ARMA_AR_MA_RepresentGrangerAnalysis,P34的例题AR(p)模型转换成MA(∞)模型。代码如下：function[psi]=ARchangeMA(AR_fi,num_psi)%输入AR_fi和num_psi，求MA模型中的psipsi=[];p=size(AR_fi,2);%AR_fi的列数if(p=0)errordlg('modelisnotcorrected');endif(p0)k=size(AR_fi,1);%AR_fi的行数p=p/k;%求pD=eye(k);I=eye(k);O=zeros(k,k*(p-1));J=[I,O];%求Jif(p==1)xi=AR_fi;J=I;elsexi=[AR_fi;[eye(k*(p-1)),zeros(k,k)]];%求xi,xi为kp×kp阶矩阵end%判断是否平稳lamta=eig(xi);fori=1:p*kif(abs(lamta(i))1)%判断此模型成立条件,lamta为xi的特征值，lamta=1/zreturnendend%求psipsi=eye(k);fori=1:num_psi-1psi=[psi,J*xi^i*J'];endend运行结果：例题3：卡尔曼滤波仿真中所选的模型的状态方程是AR(2)模型，基于老师所给数据。代码如下：clc;clearall;H=diag([111000]);G=zeros(6,6);v=[-0.017,0.016,0.013,-0.017,0.016,0.013]';phai_1=[-0.320,0.146,0.961;0.044,-0.153,0.289;-0.002,0.225,-0.264];phai_2=[-0.161,0.115,0.934;0.050,0.019,-0.010;0.034,0.355,-0.022];B=zeros(6,6);F=diag(ones(1,6));B(1:3,1:3)=phai_1;B(1:3,4:6)=phai_2;B(4:6,1:3)=diag(ones(1,3));Original_data=zeros(92,3);Original_data(:,1)=[180;179;185;192;211;202;207;214;231;229;234;237;206;250;259;263;264;280;282;292;286;302;304;307;317;314;306;304;292;275;273;301;280;289;303;322;315;339;364;371;375;432;453;460;475;496;494;498;526;519;516;531;573;551;538;532;558;524;525;519;526;510;519;538;549;570;559;584;611;597;603;619;635;658;675;700;692;759;782;816;844;830;853;852;833;860;870;830;801;824;831;830;];Original_data(:,2)=[451;465;485;493;509;520;521;540;548;558;574;583;591;599;610;627;642;653;660;694;709;734;751;763;766;779;808;785;794;799;799;812;837;853;876;897;922;949;979;988;1025;1063;1104;1131;1137;1178;1211;1256;1290;1314;1346;1385;1416;1436;1462;1493;1516;1557;1613;1642;1690;1759;1756;1780;1807;1831;1873;1897;1910;1943;1976;2018;2040;2070;2121;2132;2199;2253;2276;2318;2369;2423;2457;2470;2521;2545;2580;2620;2639;2618;2628;2651;];Original_data(:,3)=[415;421;434;448;459;458;479;487;497;510;516;525;529;538;546;555;574;574;586;602;617;639;653;668;679;686;697;688;704;699;709;715;724;746;758;779;798;816;837;858;881;905;934;968;983;1013;1034;1064;1101;1102;1145;1173;1216;1229;1242;1267;1295;1317;1355;1371;1402;1452;1485;1516;1549;1567;1588;1631;1650;1685;1722;1752;1774;1807;1831;1842;1890;1958;1948;1994;2061;2056;2102;2121;2145;2164;2206;2225;2235;2237;2250;2271;];Sample_Original_data=Original_data(4:92,:)';data=zeros(92,3);data(:,1)=log10(Original_data(:,1));data(:,2)=log10(Original_data(:,2));data(:,3)=log10(Or