自適應(yīng)控制程序

1、% M 序列及其逆序列的產(chǎn)生設(shè)M序列M(k)由如下4位移位寄存器產(chǎn)生：S(k)為方波序列，逆M序列IM(k)= M(k)ÅS(k)clear all; close all;L=60; %序列長(zhǎng)度x1=1;x2=1;x3=1;x4=0; %移位寄存器初值S=1; %方波初值for k=1:L IM=xor(S,x4); %進(jìn)行異或運(yùn)算，產(chǎn)生逆M序列 if IM=0 u(k)=-1; else u(k)=1; end S=not(S); %產(chǎn)生方波 M(k)=xor(x3,x4); %產(chǎn)生M序列 x4=x3;x3=x2;x2=x1;x1=M(k); %寄存器移位endsubplot(2,

2、1,1);stairs(M);grid;axis(0 L/2 -0.5 1.5);xlabel('k');ylabel('M序列幅值');title('M序列');subplot(2,1,2);stairs(u);grid;axis(0 L -1.5 1.5);xlabel('k');ylabel('逆M序列幅值');title('逆M序列');%白噪聲及有色噪聲序列的產(chǎn)生設(shè)x(k) 為均值為0，方差為1的高斯白噪聲序列，e(k)為有色噪聲序列： 高斯白噪聲序列 x(k)在Matlab中由rand(

3、)函數(shù)產(chǎn)生，程序如下：clear all; close all;L=500; %仿真長(zhǎng)度d=1 -1.5 0.7 0.1; c=1 0.5 0.2; % 分子分母多項(xiàng)式系數(shù)nd=length(d)-1 ;nc=length(c)-1; %階次xik=zeros(nc,1); %白噪聲初值ek=zeros(nd,1);xi=randn(L,1); %產(chǎn)生均值為0，方差為1的高斯白噪聲序列for k=1:L e(k)=-d(2:nd+1)*ek+c*xi(k);xik; %產(chǎn)生有色噪聲 %數(shù)據(jù)更新 for i=nd:-1:2 ek(i)=ek(i-1); end ek(1)=e(k); for i

4、=nc:-1:2 xik(i)=xik(i-1); end xik(1)=xi(k);endsubplot(2,1,1);plot(xi);xlabel('k');ylabel('噪聲幅值');title('白噪聲序列');subplot(2,1,2);plot(e);xlabel('k');ylabel('噪聲幅值');title('有色噪聲序列');%批處理最小二乘參數(shù)估計(jì)(LS)考慮如下系統(tǒng)：式中x(k)為方差為1的白噪聲。clear all;a=1 -1.5 0.7'b=1 0.5&

5、#39;d=3; %對(duì)象參數(shù)na=length(a)-1;nb=length(b)-1; %計(jì)算階次L=500; %數(shù)據(jù)長(zhǎng)度uk=zeros(d+nb,1);yk=zeros(na,1); %輸入初值x1=1;x2=1;x3=1;x4=0;S=1;%移位寄存器初值，方波初值xi=rand(L,1);%白噪聲序列theta=a(2:na+1);b; %對(duì)象參數(shù)真值for k=1:L phi(k,:)=-yk;uk(d:d+nb)' %phi(k,:)為行向量，便于組成phi矩陣 y(k)=phi(k,:)*theta+xi(k); %采集輸出數(shù)據(jù) IM=xor(S,x4); if IM=

6、0 u(k)=-1; else u(k)=1; end S=not(S);M=xor(x3,x4); %產(chǎn)生M序列 %更新數(shù)據(jù) x4=x3;x3=x2;x2=x1;x1=M; for i=nb+d:-1:2 uk(i)=uk(i-1); end uk(1)=u(k); for i=na:-1:2 yk(i)=yk(i-1); end yk(1)=y(k);endthetaevaluation=inv(phi'*phi)*phi'*y' %計(jì)算參數(shù)估計(jì)值thetaevaluation = -1.5362 0.6802 1.00680.4864%遺忘因子遞推最小二乘參數(shù)估計(jì)

7、(FFRLS)考慮如下系統(tǒng)：式中x(k)為均值為0、方差為0.1的白噪聲，對(duì)象時(shí)變參數(shù)為：取遺忘因子l=0.98，clear all; close all;a=1 -1.5 0.7'b=1 0.5'd=3; %對(duì)象參數(shù)na=length(a)-1;nb=length(b)-1; %計(jì)算階次L=1000;%數(shù)據(jù)長(zhǎng)度uk=zeros(d+nb,1);yk=zeros(na,1); %輸入輸出初值u=randn(L,1); %輸入采用方差為1的白噪聲序列xi=sqrt(0.1)*randn(L,1); % 方差為0.1的白噪聲干擾序列%theta=a(2:na+1);b; %對(duì)象參數(shù)

8、真值thetae_1=zeros(na+nb+1,1); %參數(shù)初值P=106*eye(na+nb+1);lambda=0.98; %遺忘因子范圍0.9 1for k=1:L if k=501 a=1 -1 0.4'b=1.5 0.2' %對(duì)象參數(shù)突變 end theta(:,k)=a(2:na+1);b; %對(duì)象參數(shù)真值 phi=-yk;uk(d:d+nb); y(k)=phi'*theta(:,k)+xi(k); %采樣輸出數(shù)據(jù) %遺忘因子遞推最小二乘公式 K=P*phi/(lambda+phi'*P*phi); thetae(:,k)=thetae_1+K

9、*(y(k)-phi'*thetae_1); P=(eye(na+nb+1)-K*phi')*P/lambda; %更新數(shù)據(jù) thetae_1=thetae(:,k); for i=d+nb:-1:2 uk(i)=uk(i-1); enduk(1)=u(k); for i=na:-1:2 yk(i)=yk(i-1); end yk(1)=y(k);endsubplot(2,1,1);plot(1:L,thetae(1:na,:);hold on;plot(1:L,theta(1:na,:),'k:');xlabel('k');ylabel(

10、9;參數(shù)估計(jì)a');legend('a_1','a_2');axis(0 L -2 2);subplot(2,1,2);plot(1:L,thetae(na+1:na+nb+1,:);hold on;plot(1:L,theta(na+1:na+nb+1,:),'k:');xlabel('k');ylabel('參數(shù)估計(jì)b');legend('b_0','b_1');axis(0 L -0.5 2); %增廣遞推最小二乘參數(shù)估計(jì)(ERLS)考慮如下系統(tǒng)：式中x(k)為方差為0

11、.1的白噪聲。選擇方差為1的白噪聲作為輸入信號(hào)u(k).clear all; close all;a=1 -1.5 0.7'b=1 0.5'c=1 -1 0.2'd=3; %對(duì)象參數(shù)na=length(a)-1;nb=length(b)-1;nc=length(c)-1; %計(jì)算階次L=1000; %數(shù)據(jù)長(zhǎng)度uk=zeros(d+nb,1);yk=zeros(na,1); %輸入輸出初值xik=zeros(nc,1); %噪聲初值xiek=zeros(nc,1); %噪聲估計(jì)初值u=randn(L,1); %輸入采用方差為1的白噪聲序列xi=sqrt(0.1)*rand

12、n(L,1); % 方差為0.1的白噪聲干擾序列theta=a(2:na+1);b;c(2:nc+1); %對(duì)象參數(shù)真值thetae_1=zeros(na+nb+1+nc,1); %參數(shù)初值,na+nb+1+nc為辨識(shí)參數(shù)個(gè)數(shù)P=106*eye(na+nb+1+nc);lambda=0.98; %遺忘因子范圍0.9 1for k=1:L phi=-yk;uk(d:d+nb);xik; %測(cè)量向量 y(k)=phi'*theta+xi(k); %采樣輸出數(shù)據(jù)phie=-yk;uk(d:d+nb);xiek; %估計(jì)的測(cè)量向量 %增廣遞推最小二乘公式 K=P*phie/(1+phie

13、9;*P*phie); thetae(:,k)=thetae_1+K*(y(k)-phie'*thetae_1); P=(eye(na+nb+1+nc)-K*phie')*P; xie=y(k)-phie'*thetae(:,k); %白噪聲估計(jì)值 %更新數(shù)據(jù) thetae_1=thetae(:,k); for i=d+nb:-1:2 uk(i)=uk(i-1); end uk(1)=u(k); for i=na:-1:2 yk(i)=yk(i-1); end yk(1)=y(k); for i=nc:-1:2 xik(i)=xik(i-1); xiek(i)=xiek

14、(i-1); end xik(1)=xi(k); xiek(1)=xie; endfigure(1)plot(1:L,thetae(1:na,:);xlabel('k');ylabel('參數(shù)估計(jì)a');legend('a_1','a_2');axis(0 L -2 2);figure(2)plot(1:L,thetae(na+1:na+nb+1,:);xlabel('k');ylabel('參數(shù)估計(jì)b');legend('b_0','b_1');axis(0 L 0

15、 1.5);figure(3)plot(1:L,thetae(na+nb+2:na+nb+nc+1,:);xlabel('k');ylabel('參數(shù)估計(jì)c');legend('c_1','c_2');axis(0 L -2 2);遞推最小二乘參數(shù)估計(jì)(RLS)考慮如下系統(tǒng)：式中x(k)為方差為0.1的白噪聲。clear all; close all;a=1 -1.5 0.7'b=1 0.5'd=3; %對(duì)象參數(shù)na=length(a)-1;nb=length(b)-1; %計(jì)算階次,na=2,nb=1L=500;

16、 %數(shù)據(jù)長(zhǎng)度(仿真長(zhǎng)度)uk=zeros(d+nb,1);yk=zeros(na,1); %輸入輸出初值uk:4x1,ykx1u=randn(L,1); %輸入采用方差為1的白噪聲序列xi=sqrt(0.1)*randn(L,1); %方差為0.1的白噪聲干擾序列theta=a(2:na+1);b; %對(duì)象參數(shù)真值theta=-1.5,0.7;1,0.5thetae_1=zeros(na+nb+1,1); %參數(shù)初值q為4x1的全零矩陣P=106*eye(na+nb+1);for k=1:L phi=-yk;uk(d:d+nb); %此處phi為列向量4x1 y(k)=phi'*the

17、ta+xi(k); %采集輸出數(shù)據(jù) %遞推公式 K=P*phi/(1+phi'*P*phi); thetae(:,k)=thetae_1+K*(y(k)-phi'*thetae_1); P=(eye(na+nb+1)-K*phi')*P; %更新數(shù)據(jù) thetae_1=thetae(:,k); for i=d+nb:-1:2 uk(i)=uk(i-1); end uk(1)=u(k); for i=na:-1:2 yk(i)=yk(i-1); end yk(1)=y(k);endplot(1:L,thetae); %line(1:L,theta,theta);xlabe

18、l('k');ylabel('參數(shù)估計(jì)a,b');legend('a_1','a_2','b_0','b_1');axis(0 L -2 2);%最小方差控制(MVC)考慮如下系統(tǒng)：式中x(k)為方差為0.1的白噪聲。取期望輸出yr(k)為幅值為10的方波信號(hào)。clear all;close all;a=1 -1.7 0.7;b=1 0.5;c=1 0.2;d=4;%對(duì)象參數(shù)na=length(a)-1;nb=length(b)-1;nc=length(c)-1;%計(jì)算階次nh=nb+d-1;%nh

19、為多項(xiàng)式H的階次L=400;uk=zeros(d+nb,1);yk=zeros(na,1);yrk=zeros(nc,1);xik=zeros(nc,1);yr=10*ones(L/4,1);-ones(L/4,1);ones(L/4,1);-ones(L/4+d,1);%期望輸出xi=sqrt(0.1)*randn(L,1);%方差為0.1的白噪聲序列h,f,g=singlediophantine(a,b,c,d);%求解單步Diophantine方程for k=1:L time(k)=k; y(k)=-a(2:na+1)*yk+b*uk(d:d+nb)+c*xi(k);xik;%采集輸出數(shù)

20、據(jù) u(k)=(-h(2:nh+1)*uk(1:nh)+c*yr(k+d:-1:k+d-min(d,nc);yrk(1:nc-d)-g*y(k);yk(1:na-1)/h(1);%求控制量 %更新數(shù)據(jù) for i=d+nb:-1:2 uk(i)=uk(i-1); end uk(1)=u(k); for i=na:-1:2 yk(i)=yk(i-1); end yk(1)=y(k); for i=nc:-1:2 yrk(i)=yrk(i-1); xik(i)=xik(i-1); end if nc>0 yrk(1)=yr(k); xik(1)=xi(k); endendsubplot(2,

21、1,1);plot(time,yr(1:L),'r:',time,y);xlabel('k');ylabel('y_r(k)、y(k)');legend('y_r(k)','y(k)');subplot(2,1,2);plot(time,u);xlabel('k');ylabel('u(k)'); 單步Diophantine方程求解求解下列系統(tǒng)的單步Diophantine方程：（1）（2）（3）%單步Diophantine方程的求解clear all;a=1 -1.5 0.7; b=

22、1 0.5; c=1; d=3; %例4.1（1）%a=1 -0.95; b=1 2; c=1 -0.7; d=2; %例4.1（2）%a=1 -1.7 0.7; b=0.9 1; c=1 2; d=4; %例4.1（3）e,f,g=sindiophantine(a,b,c,d) %調(diào)用函數(shù)sindiophantinefunction e,f,g=singlediophantine(a,b,c,d) %單步Diophantine方程求解na=length(a)-1;nb=length(b)-1;nc=length(c)-1;%計(jì)算階次ne=d-1;ng=na-1;%E,G的階次ad=a,zer

23、os(1,ng+ne+1-na);cd=c,zeros(1,ng+d-nc);%令a(na+2)=a(na+3)=.=0e(1)=1;for i=2:ne+1 e(i)=0; for j=2:i e(i)=e(i)+e(i+1-j)*ad(j); end e(i)=cd(i)-e(i);%計(jì)算eiendfor i=1:ng+1 g(i)=0; for j=1:ne+1 g(i)=g(i)+e(ne+2-1)*ad(i+j); end g(i)=cd(i+d)-g(i);%計(jì)算giendf=conv(b,e);%計(jì)算Fe = 1.0000 1.5000 1.5500f = 1.0000 2.00

24、00 2.3000 0.7750g =1.2750 -1.0850多步Diophantine方程求解求解如下系統(tǒng)的多步Diophantine方程，并去預(yù)測(cè)長(zhǎng)度N=3%多步Diophantine方程的求解clear all;a=1 -2.7 2.4 -0.7; b=0.9 1; c=1 2;na=length(a)-1; nb=length(b)-1; nc=length(c)-1; %A、B、C的階次N=3; %預(yù)測(cè)步數(shù)E,F,G=multidiophantine(a,b,c,N) %調(diào)用函數(shù)multidiophantinefunction E,F,G=multidiophantine(a,b

25、,c,N)%* %功能：多步Diophanine方程的求解 %調(diào)用格式：E,F,G=sindiophantine(a,b,c,N)（注：d=1） %輸入?yún)?shù)：多項(xiàng)式A、B、C系數(shù)向量及預(yù)測(cè)步數(shù)（共4個(gè)） %輸出參數(shù)：Diophanine方程的解E、F、G（共3個(gè)）%*na=length(a)-1; nb=length(b)-1; nc=length(c)-1; %A、B、C的階次%E、F、G的初值E=zeros(N); E(1,1)=1; F(1,:)=conv(b,E(1,:); if na>=nc G(1,:)=c(2:nc+1) zeros(1,na-nc)-a(2:na+1);

26、%令c(nc+2)=c(nc+3)=.=0else G(1,:)=c(2:nc+1)-a(2:na+1) zeros(1,nc-na); %令a(na+2)=a(na+3)=.=0end%求E、G、Ffor j=1:N-1 for i=1:j E(j+1,i)=E(j,i); end E(j+1,j+1)=G(j,1); for i=2:na G(j+1,i-1)=G(j,i)-G(j,1)*a(i); end G(j+1,na)=-G(j,1)*a(na+1); F(j+1,:)=conv(b,E(j+1,:);end%最小方差自校正控制用最小方差自校正控制算法對(duì)以下系統(tǒng)進(jìn)行閉環(huán)控制：式中x

27、(k)為方差為0.1的白噪聲。取期望輸出yr(k)為幅值為10的方波信號(hào)。%最小方差間接自校正控制clear all; close all;a=1 -1.7 0.7; b=1 0.5; c=1 0.2; d=4; %對(duì)象參數(shù)na=length(a)-1; nb=length(b)-1; nc=length(c)-1; %na、nb、nc為多項(xiàng)式A、B、C階次nf=nb+d-1; %nf為多項(xiàng)式F的階次L=400; %控制步數(shù)uk=zeros(d+nb,1); %輸入初值：uk(i)表示u(k-i);yk=zeros(na,1); %輸出初值yrk=zeros(nc,1); %期望輸出初值xik

28、=zeros(nc,1); %白噪聲初值xiek=zeros(nc,1); %白噪聲估計(jì)初值yr=10*ones(L/4,1);-ones(L/4,1);ones(L/4,1);-ones(L/4+d,1); %期望輸出xi=sqrt(0.1)*randn(L,1); %白噪聲序列%RELS初值設(shè)置thetae_1=0.001*ones(na+nb+1+nc,1); %非常小的正數(shù)（這里不能為0）P=106*eye(na+nb+1+nc);for k=1:L time(k)=k; y(k)=-a(2:na+1)*yk+b*uk(d:d+nb)+c*xi(k);xik; %采集輸出數(shù)據(jù) %遞推增

29、廣最小二乘法 phie=-yk;uk(d:d+nb);xiek; K=P*phie/(1+phie'*P*phie); thetae(:,k)=thetae_1+K*(y(k)-phie'*thetae_1); P=(eye(na+nb+1+nc)-K*phie')*P; xie=y(k)-phie'*thetae(:,k); %白噪聲的估計(jì)值 %提取辨識(shí)參數(shù) ae=1 thetae(1:na,k)' be=thetae(na+1:na+nb+1,k)' ce=1 thetae(na+nb+2:na+nb+1+nc,k)' if abs(

30、be(2)>0.9 be(2)=sign(ce(2)*0.9; %MVC算法要求B穩(wěn)定 end if abs(ce(2)>0.9 ce(2)=sign(ce(2)*0.9; end e,f,g=sindiophantine(ae,be,ce,d); %求解單步Diophantine方程 u(k)=(-f(2:nf+1)*uk(1:nf)+ce*yr(k+d:-1:k+d-min(d,nc);yrk(1:nc-d)-g*y(k);yk(1:na-1)/f(1); %求控制量 %更新數(shù)據(jù) thetae_1=thetae(:,k); for i=d+nb:-1:2 uk(i)=uk(i-

31、1); end uk(1)=u(k); for i=na:-1:2 yk(i)=yk(i-1); end yk(1)=y(k); for i=nc:-1:2 yrk(i)=yrk(i-1); xik(i)=xik(i-1); xiek(i)=xiek(i-1); end if nc>0 yrk(1)=yr(k); xik(1)=xi(k); xiek(1)=xie; endendfigure(1);subplot(2,1,1);plot(time,yr(1:L),'r:',time,y);xlabel('k'); ylabel('y_r(k)、y(

32、k)');legend('y_r(k)','y(k)'); axis(0 L -20 20);subplot(2,1,2);plot(time,u);xlabel('k'); ylabel('u(k)'); axis(0 L -40 40);figure(2)subplot(211)plot(1:L,thetae(1:na,:);xlabel('k'); ylabel('參數(shù)估計(jì)a');legend('a_1','a_2'); axis(0 L -3 2);s

33、ubplot(212)plot(1:L,thetae(na+1:na+nb+1+nc,:);xlabel('k'); ylabel('參數(shù)估計(jì)b、c');legend('b_0','b_1','c_1'); axis(0 L 0 1.5);%最小方差直接自校正控制設(shè)被控對(duì)象為：式中x(k)為方差為0.1的白噪聲。%最小方差直接自校正控制clear all; close all;a=1 -1.7 0.7; b=1 0.5; c=1 0.2; d=4; %對(duì)象參數(shù)na=length(a)-1; nb=length(b)

34、-1; nc=length(c)-1; %na、nb、nc為多項(xiàng)式A、B、C階次nf=nb+d-1; ng=na-1; %nf、ng為多項(xiàng)式F、G的階次L=400; %控制步數(shù)uk=zeros(d+nf,1); %輸入初值：uk(i)表示u(k-i);yk=zeros(d+ng,1); %輸出初值yek=zeros(nc,1); %最優(yōu)輸出預(yù)測(cè)估計(jì)初值yrk=zeros(nc,1); %期望輸出初值xik=zeros(nc,1); %白噪聲初值yr=10*ones(L/4,1);-ones(L/4,1);ones(L/4,1);-ones(L/4+d,1); %期望輸出xi=sqrt(0.1)

35、*randn(L,1); %白噪聲序列%遞推估計(jì)初值thetaek=zeros(na+nb+d+nc,d);P=106*eye(na+nb+d+nc);for k=1:L time(k)=k; y(k)=-a(2:na+1)*yk(1:na)+b*uk(d:d+nb)+c*xi(k);xik; %采集輸出數(shù)據(jù) %遞推增廣最小二乘法 phie=yk(d:d+ng);uk(d:d+nf);-yek(1:nc); K=P*phie/(1+phie'*P*phie); thetae(:,k)=thetaek(:,1)+K*(y(k)-phie'*thetaek(:,1); P=(eye

36、(na+nb+d+nc)-K*phie')*P; ye=phie'*thetaek(:,d); %預(yù)測(cè)輸出的估計(jì)值（必須為thetae(:,k-d)） %ye=yr(k); %預(yù)測(cè)輸出的估計(jì)值可取yr(k) %提取辨識(shí)參數(shù) ge=thetae(1:ng+1,k)' fe=thetae(ng+2:ng+nf+2,k)' ce=1 thetae(ng+nf+3:ng+nf+2+nc,k)' if abs(ce(2)>0.9 ce(2)=sign(ce(2)*0.9; end if fe(1)<0.1 %設(shè)f0的下界為0.1 fe(1)=0.1;

37、end u(k)=(-fe(2:nf+1)*uk(1:nf)+ce*yr(k+d:-1:k+d-min(d,nc);yrk(1:nc-d)-ge*y(k);yk(1:na-1)/fe(1); %控制量 %更新數(shù)據(jù) for i=d:-1:2 thetaek(:,i)=thetaek(:,i-1); end thetaek(:,1)=thetae(:,k); for i=d+nf:-1:2 uk(i)=uk(i-1); end uk(1)=u(k); for i=d+ng:-1:2 yk(i)=yk(i-1); end yk(1)=y(k); for i=nc:-1:2 yek(i)=yek(i-

38、1); yrk(i)=yrk(i-1); xik(i)=xik(i-1); end if nc>0 yek(1)=ye; yrk(1)=yr(k); xik(1)=xi(k); endendfigure(1);subplot(2,1,1);plot(time,yr(1:L),'r:',time,y);xlabel('k'); ylabel('y_r(k)、y(k)');legend('y_r(k)','y(k)'); axis(0 L -20 20);subplot(2,1,2);plot(time,u);x

39、label('k'); ylabel('u(k)'); axis(0 L -40 40);figure(2)subplot(211)plot(1:L,thetae(1:ng+1,:),1:L,thetae(ng+nf+3:ng+2+nf+nc,:);xlabel('k'); ylabel('參數(shù)估計(jì)g、c');legend('g_0','g_1','c_1'); axis(0 L -3 4);subplot(212)plot(1:L,thetae(ng+2:ng+2+nf,:);xl

40、abel('k'); ylabel('參數(shù)估計(jì)f');legend('f_0','f_1','f_2','f_3','f_4'); axis(0 L 0 4);%廣義最小方差控制（顯示控制）設(shè)被控對(duì)象為如下開(kāi)環(huán)不穩(wěn)定非最小相位系統(tǒng)：式中x(k)為方差為0.1的白噪聲。%廣義最小方差控制（顯式控制律）clear all; close all;a=1 -1.7 0.7; b=1 2; c=1 0.2; d=4; %對(duì)象參數(shù)na=length(a)-1; nb=length(b)-1; n

41、c=length(c)-1; %na、nb、nc為多項(xiàng)式A、B、C階次nf=nb+d-1; ng=na-1; %nf、ng為多項(xiàng)式F、G的階次P=1; R=1; Q=2; %加權(quán)多項(xiàng)式np=length(P)-1; nr=length(R)-1; nq=length(Q)-1;L=400; %控制步數(shù)uk=zeros(d+nb,1); %輸入初值：uk(i)表示u(k-i);yk=zeros(na,1); %輸出初值yrk=zeros(nc,1); %期望輸出初值xik=zeros(nc,1); %白噪聲初值yr=10*ones(L/4,1);-ones(L/4,1);ones(L/4,1);

42、-ones(L/4+d,1); %期望輸出xi=sqrt(0.1)*randn(L,1); %白噪聲序列e,f,g=sindiophantine(a,b,c,d); %求解單步Diophantine方程CQ=conv(c,Q); FP=conv(f,P); CR=conv(c,R); GP=conv(g,P); %CQ=C*Qfor k=1:L time(k)=k; y(k)=-a(2:na+1)*yk+b*uk(d:d+nb)+c*xi(k);xik; %采集輸出數(shù)據(jù) u(k)=(-Q(1)*CQ(2:nc+nq+1)*uk(1:nc+nq)/b(1)-FP(2:np+nf+1)*uk(1:

43、np+nf). +CR*yr(k+d:-1:k+d-min(d,nr+nc); yrk(1:nr+nc-d). -GP*y(k); yk(1:np+ng)/(Q(1)*CQ(1)/b(1)+FP(1); %求控制量 %更新數(shù)據(jù) for i=d+nb:-1:2 uk(i)=uk(i-1); end uk(1)=u(k); for i=na:-1:2 yk(i)=yk(i-1); end yk(1)=y(k); for i=nc:-1:2 yrk(i)=yrk(i-1); xik(i)=xik(i-1); end if nc>0 yrk(1)=yr(k); xik(1)=xi(k); end

44、endsubplot(2,1,1);plot(time,yr(1:L),'r:',time,y);xlabel('k'); ylabel('y_r(k)、y(k)');legend('y_r(k)','y(k)');subplot(2,1,2);plot(time,u);xlabel('k'); ylabel('u(k)'); %廣義最小方差自校正控制（間接算法）設(shè)被控對(duì)象為如下開(kāi)環(huán)不穩(wěn)定非最小相位系統(tǒng)：式中x(k)為方差為0.1的白噪聲。%廣義最小方差自校正控制（間接算法）clea

45、r all; close all;a=1 -1.7 0.7; b=1 2; c=1 0.2; d=4; %對(duì)象參數(shù)na=length(a)-1; nb=length(b)-1; nc=length(c)-1; %na、nb、nc為多項(xiàng)式A、B、C階次nf=nb+d-1; ng=na-1; %nf、ng為多項(xiàng)式F、G的階次Pw=1; R=1; Q=2; %加權(quán)多項(xiàng)式P、R、Qnp=length(Pw)-1; nr=length(R)-1; nq=length(Q)-1;L=400; %控制步數(shù)uk=zeros(d+nb,1); %輸入初值：uk(i)表示u(k-i);yk=zeros(na,1)

46、; %輸出初值yrk=zeros(nc,1); %期望輸出初值xik=zeros(nc,1); %白噪聲初值xiek=zeros(nc,1); %白噪聲估計(jì)初值yr=10*ones(L/4,1);-ones(L/4,1);ones(L/4,1);-ones(L/4+d,1); %期望輸出xi=sqrt(0.1)*randn(L,1); %白噪聲序列%RELS 初值設(shè)置thetae_1=0.001*ones(na+nb+1+nc,1);%非常小的正數(shù)（此處不能為0）P=106*eye(na+nb+1+nc);for k=1:L time(k)=k; y(k)=-a(2:na+1)*yk+b*uk

47、(d:d+nb)+c*xi(k);xik; %采集輸出數(shù)據(jù) %遞推增廣最小二乘法 phie=-yk;uk(d:d+nb);xiek; K=P*phie/(1+phie'*P*phie); thetae(:,k)=thetae_1+K*(y(k)-phie'*thetae_1); P=(eye(na+nb+1+nc)-K*phie')*P; xie=y(k)-phie'*thetae(:,k);%白噪聲的估計(jì)值 %提取辨識(shí)參數(shù) ae=1 thetae(1:na,k)' be=thetae(na+1:na+nb+1,k)' ce=1 thetae(n

48、a+nb+2:na+nb+1+nc,k)' if abs(ce(2)>0.9 ce(2)=sign(ce(2)*0.9; end e,f,g=sindiophantine(ae,be,ce,d); %求解單步Diophantine方程 CQ=conv(ce,Q); FP=conv(f,Pw); CR=conv(ce,R); GP=conv(g,Pw); %CQ=Ce*Q u(k)=(-Q(1)*CQ(2:nc+nq+1)*uk(1:nc+nq)/be(1)-FP(2:np+nf+1)*uk(1:np+nf). +CR*yr(k+d:-1:k+d-min(d,nr+nc); yrk

49、(1:nr+nc-d). -GP*y(k); yk(1:np+ng)/(Q(1)*CQ(1)/be(1)+FP(1);%求控制量 %更新數(shù)據(jù) thetae_1=thetae(:,k); for i=d+nb:-1:2 uk(i)=uk(i-1); end uk(1)=u(k); for i=na:-1:2 yk(i)=yk(i-1); end yk(1)=y(k); for i=nc:-1:2 yrk(i)=yrk(i-1); xik(i)=xik(i-1); xiek(i)=xiek(i-1); end if nc>0 yrk(1)=yr(k); xik(1)=xi(k); xiek(

50、1)=xie; endendfigure(1);subplot(2,1,1);plot(time,yr(1:L),'r:',time,y);xlabel('k'); ylabel('y_r(k)、y(k)');legend('y_r(k)','y(k)'); axis(0 L -20 20);subplot(2,1,2);plot(time,u);xlabel('k'); ylabel('u(k)'); axis(0 L -10 10);figure(2)plot(1:L,theta

51、e);xlabel('k'); ylabel('辨識(shí)參數(shù)a、b、c');legend('a_1','a_2','b_0','b_1','c_1'); axis(0 L -2 2.5); %廣義最小方差自校正控制（直接算法）設(shè)被控對(duì)象為如下開(kāi)環(huán)不穩(wěn)定非最小相位系統(tǒng)：式中x(k)為方差為0.1的白噪聲。%廣義最小方差自校正控制（直接算法）clear all; close all;a=1 -1.7 0.7; b=1 2; c=1 0.2; d=4; %對(duì)象參數(shù)na=length(a)-1;

52、 nb=length(b)-1; nc=length(c)-1; %na、nb、nc為多項(xiàng)式A、B、C階次nf=nb+d-1; ng=na-1; %nf、ng為多項(xiàng)式F、G的階次Pw=1; R=1; Q=2; %加權(quán)多項(xiàng)式P、R、Qnp=length(Pw)-1; nr=length(R)-1; nq=length(Q)-1;L=400; %控制步數(shù)uk=zeros(d+nf,1); %輸入初值：uk(i)表示u(k-i);yk=zeros(d+ng,1); %輸出初值yek=zeros(nc,1); %最優(yōu)輸出預(yù)測(cè)估計(jì)初值yrk=zeros(nc,1); %期望輸出初值xik=zeros(nc,1); %白噪聲初值yr=10*ones(L/4,1);-ones(L/4,1);ones(L/4,1);-ones(L/4+d,1); %期望輸出xi=sqrt(0.1)*randn(L,1); %白噪聲序列%遞推估計(jì)初值thetaek=ze