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

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

1.自校正控制

自适应控制有很多种，例如模型参考自适应控制系统、自校正控制系统等。

自校正控制(STC)最早是由R.E. Kalman 在1958年提出的，他设计了基于最小二乘估计和有限拍控制的自适应控制器，并为了实现这个控制器，还建造了一台专用模拟计算机，但其发展受到了当时的硬件问题的闲扰。

图1间接自校正控制系统 图2直接自校正控制系统 自校正控制系统也有内环和外环。内环与常规反馈系统类似，外环由对象参数递推估计器和控制器计算机构组成，其任务是由递推估计器在线估计被控对象参数，用以代替对象的未知参数，然后由设计机构按一定的规则对可调控制器的参数进行在线求解，用以修改内环的控制器。

自校正控制器是在线参数估计和控制参数在线设计两者的有机结合。另外，在参数估计时，对观测数据的使用方式有两种。一种是不直接更新控制器参数，而是先估计被控对象模型本身的未知参数，然后再通过设计机构得到控制器参数，如图 1所示，称为间接算法，另一种是直接估计控制器参数，这时需要将过程重新参数化，建立一个与控制器参数直接关联的估计模型，称为直接算法，如图2。

2.最小方差自校正Matlab 算法仿真（直接自校正和间接自校正）

设被调对象为CARMA 模型

111()()()()()()d A z y t z B z u t C z t ξ----=+

其中，

112

11

11()1 1.70.7()10.5()10.2A z z z B q z C z z -------=-+=+=+

式中，()k ξ为方差为1的白噪声。

（一）取初值6ˆ(0)10(0)0P I θ==、，0

ˆf 的下界为min 0.1f =，期望输出()r y k 为幅值为10的方波信号，采用最小方差直接自校正控制算法，观察不同时滞d=1、4、8时，最小方差自校正算法的控制效果。

(a)时滞为1 (b)时滞为4

(c)时滞为8

图4-2 最小方差直接自校正算法仿真结果

（二）取初值6ˆ(0)10(0)0.001P I θ==、。期望输出()r

y k 为幅值为10的方波信号，用增广递推最小二乘法进行参数估计，采用最小方差间接自校正控制算法，仿真效果如图4-3，d=4。

图4-3 最小方差间接自校正算法仿真效果（d=4）

Matlab 程序：

%程序：最小方差直接自校正控制（MVC ）算法

clear all; close all;

a=[1 -1.7 0.7]; b=[1 0.5]; c=[1 0.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方程

for k=1:L

time(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);%求控制量

%更新数据

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

end

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)');

subplot(2,1,2);

plot(time,u);

xlabel('k'); ylabel('u(k)');

%程序：最小方差间接自校正控制

clear all; close all;

a=[1 -1.7 0.7]; b=[1 0.5]; c=[1 0.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);

for k=1:L

time(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.9

be(2)=sign(ce(2))*0.9; %MVC算法要求B稳定

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); %求控制量

%更新数据

thetae_1=thetae(:,k);

for i=d+nb:-1:2

uk(i)=uk(i-1);

end