《智能控制系统matlab仿真》由会员分享，可在线阅读，更多相关《智能控制系统matlab仿真（9页珍藏版）》请在金锄头文库上搜索。
1、智能控制系统实验报告ARMA模型ARMA(p,q)是一个线性时间序列预测模型，适用于平稳的时间序列，即对于任何时刻t，都有,E(at)=0.协方差矩阵,对于任意有。AR模型当误差项自相关时，可以被有限阶滑动平均表示这里是零均白噪声，协方差矩阵非奇异。结合AR过程和MA误差项，得到ARMA过程：对于一个很大的阶数n，AR(n)接近ARMA(p,q)从得到残差的估计值其中利用多变量最小二乘法求解，然后使用估计值建立多变量回归模型最小二乘法求解公式，以AR(p)为例。记最小化目标函数解得以下是仿真例题：例题1：使用AR模型预测德国西部经济数据。该数据是三元时间序列，包括投资(Investment)，
2、可支配输入(income)和消费支出(consumption)三个变量。现使用AR模型预测INCOME，将1960年到1978年的75组用于训练，其余用于测试，利用最小二乘回归估计系数矩阵代码如下：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;for p=2:-1:2clear Z;clear Zc;for t=p+1:T+p Y(:,t-p)
3、=data(:,t); %构造矩阵Ymiddle1=data(:,(t-p):(t-1);middle=fliplr(middle1);Z(:,t-p)=1;middle(:); %构造矩阵Z end% -直接求伪逆%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.16
4、1,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;for t=p+1:num % for t=T+p+1:num O(:,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); %第三个到第八十九个
5、数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;if aicrecordaic recordaic=aic; m=p; Yc_b=Yc; dataf_b=dataf; end endsubp
6、lot(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.可以看出，预测效果很好。例
7、题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 (p0) k=size(AR_fi,1); %AR_fi的行数 p=p/k; %求p D=eye(k); I=eye(k); O= zeros(k,k*(p-1); J=I,O; %求J if(p=1) xi=AR_fi; J=I; else xi=AR_fi;
8、eye(k*(p-1),zeros(k,k); %求xi,xi为kpkp阶矩阵 end%判断是否平稳 lamta=eig(xi); for i=1:p*k if(abs(lamta(i)1) %判断此模型成立条件,lamta为xi的特征值，lamta=1/z return end end%求psi psi=eye(k); for i=1:num_psi-1 psi=psi,J*xii*J; endend运行结果：例题3：卡尔曼滤波仿真中所选的模型的状态方程是AR(2)模型，基于老师所给数据。代码如下：clc;clear all;H=diag(1 1 1 0 0 0);G=zeros(6,6);
9、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
10、(:,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;65
11、8;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;13
12、14;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;
13、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;210
14、2;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(Original_data(:,3);Difference_data=zeros(91,3);for i=2:92Difference_data(i-1,:)=data(i,:)-data(i-1,:);endStudy_data=zeros(6,89); %跟踪的信号Output_data=
