Goodwin模型参考自适应控制方案

本篇博客简要介绍Goodwin控制方案，并进行了Matlab实现。

      博主在完成自适应控制小课程设计时需要用到Goodwin MRAC控制方案，但是在国内外网站上均没有找到简单易懂的参考资料，而学校用到的教材中也只是简要介绍甚至还有错误。因此，仔细拜读了Goodwin大牛编写的教材《自适应滤波预测与控制》明白了一二，故写下本篇博客进行记录。

      首先，先介绍我们需要完成的控制任务

        在Goodwin教授的书上模型参考控制的结构图如下，我们这里只要求跟踪输入信号，故G=H/E为1.

        系统对象参数未知，因此要用到自适应预测器估计参数值。Goodwin教授介绍了两种估计方法，直接估计和间接估计，而博主课上只提到了直接估计，因此选择直接估计法里不受限的情况。

        直接算法构建自适应预测器，当预测输出是参数θ的线性函数时，对于不受限的情况，预测器的输出可描述为：

        这个式子是通过离散系统推导定义出来的，具体过程参见《自适应技术的理论及应用(第二版)》3.5节离散时域模型参考自适应控制系统。

        此处注意：在本书此节有一处错误和一处未点清的地方

        ①下图中a0的位置应该是a1

         ②在求gi时，如果i+d-j>n,ai+d-j认为是零

        其中信号向量：

         基于投影算法自适应d步超前预测器为：

        对于给定的参考序列，要求

        定义输出误差为：

        控制规律如下，可导致输出误差为零：

        完整代码如下：

        （可任意更改参考模型阶数及延迟时间）

Author:GBH
Date:2022.5.8
clear;clc;
%仿真步长及仿真次数
h=0.1; L=100/h;
m=4; n=4; d=2;
%测试信号-参考输出
r=1;     %幅值为r
t=1:1000+d;
%y_ref=r*square(t/250*pi);  %周期为250 4pi/(pi/125)
y_ref=sin(t/250*pi);  %周期为250 2pi/(pi/125)
%初始化，移迟部分设置为0
for i=1:d
    u(i)=0;
    y(i)=0;
    Ceta(:,i)=zeros(m+n+d,1);
    Fai_turn(i,:)=zeros(1,m+n+d);
end
%求解Fai(k)的转置
%Fai_turn(k)=[y(k),y(k-1),```,y(k-n+1),u(k),u(k-1),```,u(k-m-d+1)]
for k=d+1:L
    %Fain_turn(k);
    for i=1:n
        if((k-i+1<1)||(i==1))   %k-i+1<0即y(k)中k<0;i==1即y(k).两者分别为不存在或未知的因此设为零
             Fai_turn(k,i)=0;
        else
            Fai_turn(k,i)=y(k-i+1);
        end
    end
    for i=1:m+d
         if((k-i+1<1)||(i==1))
            Fai_turn(k,n+i)=0;
         else
            Fai_turn(k,n+i)=u(k-i+1);
         end
    end
    %根据Fai_turn(k)、Ceta(k)、y(k)求解Ceta(k)，设a(k)=1，c=1
    Ceta(:,k)=Ceta(:,k-1)+Fai_turn(k-d,:)'*(y_ref(k)-Fai_turn(k-d,:)*Ceta(:,k-1))/(1+Fai_turn(k-d,:)*Fai_turn(k-d,:)');
    %输入差值为：参考数值-预测输出，设Beta(k)=1
    u(k)=y_ref(k+d)-[Fai_turn(k,1:n-1),Fai_turn(k,n+1:m+n+d)]*[Ceta(1:n-1,k);Ceta(n+1:m+n+d,k)];
    %将缺失的beta0填上
    Fai_turn(k,n)=u(k);
    y(k+d)=Fai_turn(k,:)*Ceta(:,k);    %预测器
end
%图像显示
%合并显示
figure(1);
plot(t,y_ref,'b',t,y,'r');set(gca,'YLim',[-1.5 1.5]);
%分开显示
figure(2);
subplot(2,1,1);
plot(t,y_ref,'b');set(gca,'YLim',[-1.5 1.5]);
subplot(2,1,2);
plot(t,y,'r');set(gca,'YLim',[-1.5 1.5]);

结果（蓝色的为输入信号，红色的为跟踪信号）：