最小方差自校正控制Matlab程序

最小方差自校正控制Matlab程序最小方差自校正控制Matlab程序1.自校正控制自适应控制有很多种，例如模型参考自适应控制系统、自校正控制系统等。自校正控制(STC)最早是由R.E.Kalman在1958年提出的，他设计了基于最小二乘估计和有限拍控制的自适应控制器，并为了实现这个控制器，还建造了一台专用模拟计算机，但其发展受到了当时的硬件问题的闲扰。图1间接自校正控制系统图2直接自校正控制系统自校正控制系统也有内环和外环。内环与常规反馈系统类似，外环由对象参数递推估计器和控制器计算机构组成，其任务是由递推估计器在线估计被控对象参数，用以代替对象的未知参数，然后由设计机构按一定的规则对可调控制器的参数进行在线求解，用以修改内环的控制器。自校正控制器是在线参数估计和控制参数在线设计两者的有机结合。另外，在参数估计时，对观测数据的使用方式有两种。一种是不直接更新控制器参数，而是先估计被控对象模型本身的未知参数，然后再通过设计机构得到控制器参数，如图1所示，称为间接算法，另一种是直接估计控制器参数，这时需要将过程重新参数化，建立一个与控制器参数直接关联的估计模型，称为直接算法，如图2。2.最小方差自校正Matlab算法仿真（直接自校正和间接自校正）设被调对象为CARMA模型111()()()()()()dAzytzBzutCztξ----=+其中，1121111()11.70.7()10.5()10.2AzzzBqzCzz-------=-+=+=+式中，()kξ为方差为1的白噪声。（一）取初值6?(0)10(0)0PIθ==、，0f的下界为min0.1f=，期望输出()ryk为幅值为10的方波信号，采用最小方差直接自校正控制算法，观察不同时滞d=1、4、8时，最小方差自校正算法的控制效果。(a)时滞为1(b)时滞为4(c)时滞为8图4-2最小方差直接自校正算法仿真结果（二）取初值6?(0)10(0)0.001PIθ==、。期望输出()ryk为幅值为10的方波信号，用增广递推最小二乘法进行参数估计，采用最小方差间接自校正控制算法，仿真效果如图4-3，d=4。图4-3最小方差间接自校正算法仿真效果（d=4）Matlab程序：%程序：最小方差直接自校正控制（MVC）算法clearall;closeall;a=[1-1.70.7];b=[10.5];c=[10.2];d=4;%对象参数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);%白噪声初值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);%白噪声序列[e,f,g]=sindiophantine(a,b,c,d);%求解单步Diophantine方程fork=1:Ltime(k)=k;y(k)=-a(2:na+1)*yk+b*uk(d:d+nb)+c*[xi(k);xik];%采集输出数据u(k)=(-f(2:nf+1)*uk(1:nf)+c*[yr(k+d:-1:k+d-min(d,nc));yrk(1:nc-d)]-g*[y(k);yk(1:na-1)])/f(1);%求控制量%更新数据fori=d+nb:-1:2uk(i)=uk(i-1);enduk(1)=u(k);fori=na:-1:2yk(i)=yk(i-1);endyk(1)=y(k);fori=nc:-1:2yrk(i)=yrk(i-1);xik(i)=xik(i-1);endifnc>0yrk(1)=yr(k);xik(1)=xi(k);endendsubplot(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)');%程序：最小方差间接自校正控制clearall;closeall;a=[1-1.70.7];b=[10.5];c=[10.2];d=4;%对象参数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);fork=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=[1thetae(1:na,k)'];be=thetae(na+1:na+nb+1,k)';ce=[1thetae(na+nb+2:na+nb+1+nc,k)'];ifabs(be(2))>0.9be(2)=sign(ce(2))*0.9;%MVC算法要求B稳定endifabs(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);fori=d+nb:-1:2uk(i)=uk(i-1);enduk(1)=u(k);fori=na:-1:2yk(i)=yk(i-1);endyk(1)=y(k);fori=nc:-1:2yrk(i)=yrk(i-1);xik(i)=xik(i-1);xiek(i)=xiek(i-1);endifnc>0yrk(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(k)');legend('y_r(k)','y(k)');axis([0L-2020]);subplot(2,1,2);plot(time,u);xlabel('k');ylabel('u(k)');axis([0L-4040]);figure(2)subplot(211)plot([1:L],thetae(1:na,:));xlabel('k')