最小方差控制

最小方差自校正调节器自校正调节器(Self-TuningRegulator,STR)最早是由Astrom和Wittenmark于1973年首先提出来的,其结构如图1所示.•STR是以RLS参数估计方法在线估计最优预报模型,并在此基础上以输出方差最小为调节指标的一种可以适应参数未知或慢时变的自适应控制系统.•欲讨论参数未知时能调节系统输出方差至最小的STR,需先引入参数已知时调节系统输出方差最小的最小方差调节器.•最小方差调节的基本思想是:–由于系统中信道存在着d步时滞,这就使得当前的控制作用u(k)要到d个采样周期后才能对输出产生影响.–因此,要获得输出方差最小,就必须对输出量提前d步进行预报,然后根据预报值来计算适当的调节作用u(k).–这样,通过不断的预报和调节,就能始终保持输出量的稳态方差为最小.在最小方差调节器的研究中,所讨论的被控系统的模型为A(q-1)y(k)=B(q-1)u(k-d)+C(q-1)e(k)(1)nnnnnnqcqcqCqbqbbqBqaqaqA11111011111)()(1)(其中：对该系统,有如下假设:1.被控系统时滞时间d以及时滞算子q-1的多项式A、B和C的阶次及系数都已知;2.被控系统为最小相位系统,即多项式B(q-1)的所有零点都在单位圆内;3.A(q-1)、C(q-1)所有零点都为稳定的,即所有零点都在单位圆内;4.{e(k)}为零均值白色噪声序列,且E{e2(k)}=2.1最小方差预测设在k时刻已观测到输出值y(k),y(k-1),…等，希望由此得到预测值。)(ˆkdky)()()()()()()(1111dkeqAqCkuqAqBdky右边为u(k),u(k-1),…；e(k+d),e(k+d-1),…,e(k+1),e(k),e(k-1),…,等变量的线性组合，u(k),u(k-1)……为系统输入，可直接测量。随机变量e(k),e(k-1)……可根据在k时刻为止的系统输入输出值计算得到。e(k+1),……,e(k+d)为系统在k时刻以后的干扰输入，与直到k时刻为止的系统输入输出观测值无关。为了有效利用直到k时刻为止的系统输入输出观测值进行预测，须将这两类变量区分开。为此可将C(q-1)/A(q-1)分成两部分：(2)由式(1)有：)1(11101)1(111111111)(1)()()()()()(nndddqgqggqGqfqfqFqAqGqqFqAqCF和G可通过长除法得到，F为商，而q-dG(q-1)为余因子。也可通过将式(3)写成(3))()()()(1111qGqqFqAqCd然后比较两边系数得到。由(3)可将(2)右边的噪声项写成：(4))()()()()()()()(11111keqAqGdkeqFdkeqAqC(5)代入(2)可得：)()()()()()()()()(11111keqAqGdkeqFkuqAqBdky而由式(1)有：)()()()()()()(1111kuqqCqBkyqCqAked(6)代入式(6)有：)()()()(])()()()()()([)()()(111111111kyqCqGkuqCqAqBqqGqAqBdkeqFdkyd(7)利用式(4)可将式(8)化简为：(8))()()()()()()()()()(111111kyqCqGkuqCqFqBdkeqFdky(9)记基于k时刻的观测值对y(k+d)的预报为：则它是k时刻及以前的输入输出的函数。若对预测的要求是使预测的误差平方即系统误差的方差为最小，则损失函数可表示为：)(ˆkdky})]()()()()()()()()()(ˆ{[})]()(ˆ{[21111112kyqCqGkuqCqFqBdkeqFkdkyEdkykdkyEJ上式中F(q-1)e(k+d)与其它项均不相关，且由于{e(k)}为零均值白噪声序列，式(10)可写为})]()()()()()()()(ˆ{[})]()({[21111121kyqCqGkuqCqFqBkdkyEdkeqFEJ(10)与的选择无关)(ˆkdky因此当上式中第2项为0时，可使J最小。因此最小方差预测为：(11))()()()()()()()(ˆ11111kyqCqGkuqCqFqBkdky(12)最小方差预测估计的误差的方差为2212121)1(})]()({[)}(~var{dffdkeqFEkdky(13))(ˆ-)()(~kdkydkykdky2最小方差控制最小方差控制的目的是要确定u(k)，使得输出的方差为最小，由于u(k)最早只能影响到y(k+d)，因此选择性能指标为})({2dkyEJ上式可改写为：(14)})(~{})(ˆ{})](~)(ˆ{[})({2222kdkyEkdkyEkdkykdkyEdkyEJ(15)预测误差，e(k+1),…,e(k+d)的线性组合。显然，使式(15)中性能指标取最小值的充要条件是：0)()()()()()()()(ˆ11111kyqCqGkuqCqFqBkdky(16)因此最小方差控制律为：)()(-)()()(111kyqGkuqFqB此时系统输出的方差为：2212122)1(})(~{})({dffkdkyEdkyE(17)(18)由式(16)可见，最小方差控制律可以通过先求出输出提前d步的预测值，然后令等于理想输出值yr(这里yr=0)而得到，因此最小方差控制问题可分离成两个问题，一个是预测问题，另一个是控制问题。)(ˆkdky)(ˆkdky•最小方差自校正控制器•在自校正调节器中，通常认为参考输入为零，性能指标为输出方差最小，对控制信号没有加以限制。•按照这种目标函数设计的自校正调节器不适应用于非逆稳定的受控系统，为克服这些缺点D.W.Clarke和P.I.Gawthrop等人于1975年提出了自校正控制器的算法。•在实际的控制系统中，经常会遇到跟踪问题，希望系统的输出y(t)能很好的跟踪参考输入yr(t)。另一方面，希望对控制信号加以一定限制。•因此提出了最小方差自校正控制问题。它仍然选用二次型指标函数，但在指标函数中引入了参考输出项和控制作用的加权项。•目标函数：•其中都是时滞算子q-1的多项式，分别为系统输出y(t)，伺服输入yr(t)和控制作用u(t)的加权因子，其选取按对跟踪精度的要求和对u(t)的限制而定。•记为•预测模型与自校正调节器的一样，为：•将预测模型代入目标函数：•使得J最优的控制律•目标函数•该目标函数又称广义最小方差准则，因此这里所讨论的最小方差控制又称广义最小方差自校正控制。•最小方差自校正控制器的稳定性问题•B的零点成为闭环特征方程的根，如果B为非稳定多项式，即控制对象为非逆稳，闭环系统就不稳定了。但可以通过合理选择权因子，保证系统稳定。•仿真程序清单：%*************系统初始化************%a=[1-1.70.7];b=[10.5];c=[10.2];k=4;%对象参数na=length(a)-1;nb=length(b)-1;nc=length(c)-1;%na、nb、nc为多项式A、B、C阶次nd=k-1;nf=nb+nd;ne=na-1;%nd、nf、ng为多项式D、F、G的阶次K=400;%控制步数ut=zeros(k+nf,1);%输入历史的初值yt=zeros(k+ne,1);%输出历史的初值yet=zeros(nc,1);%输出估计历史的初值yrt=zeros(nc,1);%期望输出历史初值wt=zeros(nc,1);%白噪声历史的初值wet=zeros(nc+k-1,1);%白噪声估计历史的初值yr=10*[ones(K/2,1);-ones(K/2+k,1)];%期望输出w=sqrt(0.1)*randn(K,1);%白噪声序列%递推估计初值thetae_1=0.001*ones(na+nb+k+nc+nc+nd,1);P=10^8*eye(na+nb+k+nc+nc+nd);%**********递推辨识和控制*********%fort=1:Ky(t)=-a(2:na+1)*yt(1:na)+b*ut(k:k+nb)+c*[w(t);wt];%采集输出数据phie=[yt(k:k+ne);ut(k:k+nf);-yet(1:nc);wet];ye=phie'*thetae_1;%输出的估计值，用之前的参数向量计算%递推增广最小二乘法C(Z)*y(t)=E(Z)*y(t-k)+F(Z)*u(t-k)+C(Z)*D(Z)*w(t)%实际用此式进行计算式，使用t-k时刻的数据去计算t时刻的输出，故y(t-1)...y(t-nc)%是不知道的，故用估计输出代替L=P*phie/(1+phie'*P*phie);thetae(:,t)=thetae_1(:,1)+L*(y(t)-phie'*thetae_1(:,1));P=(eye(na+nb+k+nc+nc+nd)-L*phie')*P;thetae_1=thetae(:,t);%参数向量更新we=y(t)-phie'*thetae_1;%白噪声的估计值%提取辨识参数ge=thetae_1(1:ne+1)';fe=thetae_1(ne+2:ne+nf+2)';ce=[1thetae_1(ne+nf+3:ne+nf+2+nc)'];ifabs(ce(2))0.9%稳定边界为1，取0.9ce(2)=sign(ce(2))*0.9;endiffe(1)0.1%设f0的下界为0.1fe(1)=0.1;end%求控制量u(k)=[C(Z)*yr(k+d)-E(Z)*y(k)]/[B(Z)*D(Z)]u(t)=(-fe(2:nf+1)*ut(1:nf)+ce*[yr(t+k:-1:t+k-min(k,nc));yrt(1:nc-k)]-ge*[y(t);yt(1:na-1)])/fe(1);%更新数据fori=k+nf:-1:2ut(i)=ut(i-1);endut(1)=u(t);fori=k+ne:-1:2yt(i)=yt(i-1);endyt(1)=y(t);fori=nc+nd:-1:2wet(i)=wet(i-1);endwet(1)=we;fori=nc:-1:2yet(i)=yet(i-1);yrt(i)=yrt(i-1);wt(i)=wt(i-1);endifnc0yet(1)=ye;yrt(1)=yr(t);wt(1)=w(t);endend%************显示图像***************%figure(1);subplot(2,1,1);plot([1:K],yr(1:K),'r:',[1:K],y);xlabel('k');ylabel('y_r(k),y(k)');legend('y_r(k)','y(k)');axis([0K-2020]);subplot(2,1,2);plot([1:K],u);xlabel('k');ylabel('u(k)');axis([0K-4040]);figure(2)subplot(2,1,1)plot([1:K],thetae(1:ne+1,:),[1:K],thetae(ne+nf+3:ne+2+nf+nc,:));xlabel('k');ylabel('参数估计e,c');legend('e_0','e_1','c_1');axis([0K-34]);subplot(2,1,2)plot([1:K],thetae(ne+2:ne+2+nf,:));xla