4 仿真分析
广义最小方差间接自校正控制
参数初始化：
a=[1 -1.4 0.5]; b=[0.001 0.001]; c=[1 0.3]; d=2; %对象参数
na=length(a)-1; nb=length(b)-1; nc=length(c)-1; %na、nb、nc为多项式A、B、C阶次
nf=nb+d-1; %nf为多项式F的阶次控制步数和初值设置
L=400; %控制步数
uk=zeros(d+nb,1); %输入初值：uk(i)表示u(k-i);
yk=zeros(na,1); %输出初值
yrk=zeros(nc,1); %期望输出初值
xik=zeros(nc,1); %白噪声初值
xiek=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); %白噪声序列RELS初值设置
thetae_1=0.001*ones(na+nb+1+nc,1); %非常小的正数（这里不能为0）
P=10^6*eye(na+nb+1+nc);主循环
for k=1:Ltime(k)=k;y(k)=-a(2:na+1)*yk+b*uk(d:d+nb)+c*[xi(k);xik]; %采集输出数据%递推增广最小二乘法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); %白噪声的估计值%提取辨识参数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(be(2))>0.9be(2)=sign(ce(2))*0.9; %MVC算法要求B稳定endif abs(ce(2))>0.9ce(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); %求控制量%更新数据thetae_1=thetae(:,k);for i=d+nb:-1:2uk(i)=uk(i-1);enduk(1)=u(k);for i=na:-1:2yk(i)=yk(i-1);endyk(1)=y(k);for i=nc:-1:2yrk(i)=yrk(i-1);xik(i)=xik(i-1);xiek(i)=xiek(i-1);endif nc>0yrk(1)=yr(k);xik(1)=xi(k);xiek(1)=xie;end
end绘图：
figure(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 -40 40]);
figure(2)
subplot(211)
plot([1:L],thetae(1:na,:));
xlabel('k'); ylabel('参数估计a');
legend('a_1','a_2'); axis([0 L -3 2]);
subplot(212)
plot([1:L],thetae(na+1:na+nb+1+nc,:));
xlabel('k'); ylabel('参数估计b、c');
legend('b_0','b_1','c_1'); axis([0 L 0 1.5]);辅助函数：
单步Diophanine方程的求解！！！！这居然是个常用的自定义函数
function [d,f,e]=sindiophantine(a,b,c,k)
%*********************************************************%功能：单步Diophanine方程的求解%调用格式：[d,f,e]=sindiophantine(a,b,c,t)%输入参数：多项式A、B、C系数（行向量）及纯滞后（共4个）%输出参数：Diophanine方程的解d,f,e（共3个）% C = AD + zdE% F = BD
%*********************************************************
na=length(a)-1; nb=length(b)-1; nc=length(c)-1; %A、B、C的阶次
nd=k-1; ne=na-1; %D、E的阶次
ad=[a,zeros(1,ne+nd+1-na)]; cd=[c,zeros(1,ne+nd+1-nc)]; %ad、cd长度均设为ne+nd+2%短除法求解D(Z),E(Z)
for i=1:nd+1d(i)=0;for j=1:id(i)=d(i)+d(i+1-j)*ad(j);endd(i)=cd(i)-d(i); %计算ei
endfor i=1:ne+1e(i)=0;for j=1:nd+1e(i)=e(i)+d(nd+2-j)*ad(i+j);ende(i)=cd(i+nd+1)-e(i); %计算gi
end%得到F(Z)
f=conv(b,d); %计算F
end直接和间接控制反映在代码上的区别：
参数设置：
间接控制（需要进行参数设置）：
% 递推增广最小二乘法
phie = [-yk; uk(d:d+length(b)); xiek];
K = P * phie / (1 + phie' * P * phie);
thetae(:, k) = thetae_1 + K * (y(k) - phie' * thetae_1);
P = (eye(length(a)+length(b)+1+length(c)) - K * phie') * P;% 提取辨识参数
ae = [1 thetae(1:length(a), k)'];
be = thetae(length(a)+1:length(a)+length(b)+1, k)';
ce = [1 thetae(length(a)+length(b)+2:end, k)'];直接控制（没有参数估计过程，直接使用已知的系统参数）
