《智能控制系统matlab仿真(共9页).doc》由会员分享，可在线阅读，更多相关《智能控制系统matlab仿真(共9页).doc（9页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、精选优质文档-倾情为你奉上智能控制系统实验报告ARMA模型ARMA(p,q)是一个线性时间序列预测模型，适用于平稳的时间序列，即对于任何时刻t，都有,E(at)=0.协方差矩阵,对于任意有。AR模型(0.1)当误差项自相关时，可以被有限阶滑动平均表示(0.2)这里是零均白噪声，协方差矩阵非奇异。结合AR过程和MA误差项，得到ARMA过程：(0.3)对于一个很大的阶数n，AR(n)接近ARMA(p,q)(0.4)从得到残差的估计值(0.5)其中利用多变量最小二乘法求解，然后使用估计值建立多变量回归模型(0.6)(0.7)(0.8)(0.9)(0.10)最小二乘法求解公式，以AR(p)为例。(0.
2、11)记(0.12)(0.13)最小化目标函数(0.14)解得(0.15)以下是仿真例题：例题1：使用AR模型预测德国西部经济数据。该数据是三元时间序列，包括投资(Investment)，可支配输入(income)和消费支出(consumption)三个变量。现使用AR模型预测INCOME，将1960年到1978年的75组用于训练，其余用于测试，利用最小二乘回归估计系数矩阵代码如下：专心-专注-专业clearclcdata1=load(e1.txt); %data1=data1;data2=log(data1);data3=data2(:,2:end)-data2(:,1:end-1); %差
3、分数据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)=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(:),k
4、ron(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;for t=p+1:num % for t=T+p+1:num O(:,t-p)=data(:,t);middle2=data(:,(t-p):(t-1);middlec=flipl
5、r(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
6、;aic=log(det(sigmaU)+2*p*K*K/T;if aicrecordaic recordaic=aic; m=p; Yc_b=Yc; dataf_b=dataf; end endsubplot(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
7、-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 (p0) k=size(AR_fi,1); %AR_fi的
8、行数 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;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; e
9、ndend运行结果：例题3：卡尔曼滤波仿真中所选的模型的状态方程是AR(2)模型，基于老师所给数据。代码如下：clc;clear all;H=diag(1 1 1 0 0 0);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)
10、;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;
11、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;8
12、08;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;
13、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;1
14、485;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(Origi
15、nal_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=zeros(6,89);Estimate_data_t_t_1=zeros(6,89);Estimate_data_t_t=zeros(6,89);Estimate_Smoothing=zeros(6,89);a=1.0e-003 *0.3922 0.0154 0.0236; 0.0154 0.0274 0.012
16、2; 0.0236 0.0122 0.0168;Esti_sigma=diag(ones(1,6);Esti_sigma(1:3,1:3)=a;Esti_sigma_y=zeros(6,6); Esti_sigma_y(1:3,1:3)=(0.1)*a;Esti_sigma_y(4:6,4:6)=a;Esti_sigma_out=zeros(6,6);Esti_sigma_t_1=zeros(6,6);Esti_sigma_t=0.004*diag(ones(1,6);Esti_eror=zeros(3,89);kalman_p=zeros(6,6);%相关参数for k=1:89 Study
17、_data(1:3,k)=Difference_data(k+2,:); Study_data(4:6,k)=0 0 0; %Study_data(4:6,k)=Difference_data(k+1,:); %Estimate_data_t_t_1(:,k)=B*Study_data(:,k)+F*v;endfor t=1:89 if t=1 Estimate_data_t_t_1(:,t)=v; %Esti_sigma_t_1=Esti_sigma; else Estimate_data_t_t_1(:,t)=B*Estimate_data_t_t(:,t-1)+v; %Esti_sigm
18、at_1=B*Esti_sigma_t+Esti_sigma; end Esti_sigma_t_1=B*Esti_sigma_t+Esti_sigma; Output_data(:,t)=H*Estimate_data_t_t_1(:,t); Esti_sigma_out=H*Esti_sigma_t_1*H+Esti_sigma_y+0.05*Esti_sigma; %Esti_sigma_out=H*Esti_sigma_t_1*H+Esti_sigma_y+Esti_sigma; kalman_p=Esti_sigma_t_1*(H)*(inv(Esti_sigma_out); Est
19、imate_data_t_t(:,t)=Estimate_data_t_t_1(:,t)+kalman_p*(Study_data(:,t)-Output_data(:,t); Esti_sigma_t=Esti_sigma_t_1-kalman_p*Esti_sigma*kalman_p; Esti_error(:,t)=Study_data(1:3,t)-Output_data(1:3,t);endkalman_pfiguresubplot(221)stairs(Output_data(1,:),r);hold onstairs(Study_data(1,:),b);grid on;leg
20、end(预报值,真实值);xlabel(投资);subplot(222)stairs(Output_data(2,:),r);hold on;stairs(Study_data(2,:),b);grid on;legend(预报值,真实值);xlabel(收入);subplot(223)stairs(Output_data(3,:),r);hold on;stairs(Study_data(3,:),b);grid on;legend(预报值,真实值);xlabel(消费);subplot(224)stairs(Esti_error(1,:),r);hold on;stairs(Esti_error(2,:),b);hold on;stairs(Esti_error(3,:),b:);grid on;legend(投资误差,输入误差,消费误差);运行结果：训练出的卡尔曼滤波的系数矩阵：实验分析：由于此仿真中的输出方程为随意选定，输出方程的随机干扰噪声与实际偏差较大，故此仿真的跟踪效果不是很理想。