【Matlab】极点配置控制(PPC)

前言

前段时间一直在做系统辨识相关的东西，但知道了传递函数，也要会分析传递函数，设计控制器，或者根据传递函数来计算PID参数嘛。所以对于传递函数的分析是比较重要的，因此，后续我学习了一下这个极点配置控制(PPC)，然后随便举了一个测试传递函数的例子，自己按照书上的步骤计算了一下，然后仿真效果，看起来还行吧，但是还是有很多疑点的，本节文章，只是记录一下极点配置控制的步骤，细节部分还需要持续的对与控制理论的学习。

设计极点配置控制的控制律为：

F(z^-1) * u(k) = R(z^-1) * yr(k+d) - G(z^-1) * y(k)

其中F(z^-1),R(z^-1),G(z^-1)是待定多项式，yr(k)为参考输入

极点配置控制的设计任务是：

根据系统的固有特性和控制要求，确定期望的闭环特征多项式 T(z^-1)，进而通过公式：

A*F + z^-d * B * G = T 来确定 F(z^-1) ，G(z^-1)，最后由公式：

F(z^-1) * u(k) = R(z^-1) * yr(k+d) - G(z^-1) * y(k) 来计算出控制量

极点配置控制系统框图：

使用一个传递函数测试，步骤如下：

假设传递函数为：

sys =
 
     s + 1
  -----------
  s^2 + s + 1

按照采样周期 Ts = 0.5 转换成z域为：

dsys = c2d(sys,0.5,'zoh')

dsys =
 
    0.4818 z - 0.2891
  ----------------------
  z^2 - 1.414 z + 0.6065
 
Sample time: 0.5 seconds
Discrete-time transfer function.

系统有一个零点 z1 = 0.6

假设期望的闭环系统为自然频率 ω = 1 rad/s 和阻尼系数 ξ = 0.7 的二阶系统，则其闭环特征多项式为：T(s) = s^2 + 2 * ω * ξ + ω^2 = s^2 + 1.4*s + 1 （注：特征多项式是分母）

>> sys = tf([1],[1,1.4,1])

sys =
 
         1
  ---------------
  s^2 + 1.4 s + 1
 
Continuous-time transfer function.

同样按照采样周期 Ts = 0.5 进行离散化：

>> dsys = c2d(sys,0.5,'zoh')

dsys =
 
   0.09833 z + 0.07778
  ---------------------
  z^2 - 1.32 z + 0.4966
 
Sample time: 0.5 seconds
Discrete-time transfer function.

之后可以转化成 G(z^-1)的形式，如下：

可以得到阶次 d 的值：分母最高阶次为2 分子最高阶次为1 d = 2 - 1 = 1，不过也可以不计算，因为后续的计算只是用来判断观测器的最大阶数范围。

接下来开始计算期望的闭环系统的参数，期望的闭环系统设计为 自然频率 ω = 1 rad/s 和阻尼系数 ξ = 0.7 的二阶系统

则期望的闭环特征多项式为：T(s) = s^2 + 2*ω*ξ*s + ω^2

期望闭环多项式的极点为：p1 = 0.6602 + j0.2463 和 p2 = 0.6602 - j0.2463

期望闭环多项式为：

Am(z) = z^2 - 1.32z + 0.4966

目标矩阵有一个零点，此时，B矩阵可以选择是否对消零点，因此需要降B矩阵分解成 B- * B+的形式，B+即为需要对消的零点，B-为不需要对消的零点，如果选择不对消零点，则B+ = 1：

这里我选择了，不对消零点，则B+ = 1，B-即为B矩阵本身。B矩阵不需要分解，因此 B(z) = 0.4818 z - 0.2891

B = dnum(2:na+1);    %取出分子系数中大于等于2阶的项  因为期望闭环系统为2阶

B =

    0.4818   -0.2891

为使稳态输出无误差，取Bm‘ = Am(1) / B-(1) = 0.1766 / 0.1927 = 0.9140

Bm1 = sum(Am) / sum(b0);  %求Bm'= Am(1) / B-(1)    Am(1)即为各个系数之和

Bm1 =

    0.9140

Bm(z) = Bm'(z) * B(z) = 0.9140* (0.4818 z - 0.2891) = 0.4403 z - 0.2642

Bm = Bm1 * B; %确定多项式Bm = Bm' * B-   即各个系数之和乘以可对消零点


Bm =

    0.4403   -0.2642

之后是求观测器多项式A0，根据各个矩阵的阶次，可以得出观测多项式A0的阶数

degA0 ≤ degA + degB + d - 1 - degAm = 2 + 1 + 1 - 1 - 2 = 1

degA0 ≥ 2 * degA - degAm - degB - 1= 2 * 2 - 2 - 1 - 1 = 0

取degA0 = 1 取A0 = z + 0.5

%****************************************确定观测器多项式A0******************************************************%
%  观测器A0阶次需满足：
%     degA0 ≤ degA + degB- + d - 1 - degAm    (d为传递函数最高阶次)
%     degA0 ≥ 2 * degA - degAm - degB+ - 1
%***************************************************************************************************************%
na0 = 2 * na - 1 - nam - (length(b1)-1); %观测器直接取最低阶次  为1 
A0 = 1; %取A0 = 1
for i=1:1
    A0 = conv(A0,[1 0.5]); %生成观测器
end

之后根据方程组最高次幂相等的原理，求出 F1,G,R：

A*F1 + z^-d * B * G = A0 * Am

R = Bm' * A0

%******************************计算Diophantine方程，得到F、G、R********************************************%
[F1,G] = diophantine(A,B,d,A0,Am); %注意，此处为b0
F = conv(F1,1); 
R = Bm1 * A0; 

F =

    1.0000   -0.0018
    
G =

    1.2353   -0.8627 

R矩阵很好求得：

R = 0.9165 * (z + 0.5) = 0.9140 z + 0.4570

R =

    0.9140    0.4570

我们可以按照公式来验证计算出来的F，G，R矩阵是否正确，即：

A * F1 + z^-d * B * G = A0 * Am

A0 * Am = (z + 0.5) * (z^2 - 1.32 z + 0.4966) = z^3 + 1.32 z^2 + 0.4966 z + 0.5 z^2 - 0.66 z + 0.2483

= z^3 - 0.8205 z^2 - 0.1634z + 0.2483

A0Am = conv(A0,Am);  %求多项式相乘系数

A0Am =

    1.0000   -0.8205   -0.1637    0.2483

A * F1 + B * G = (z^2 - 1.414 z + 0.6065) * (z - 0.0018) + (0.4818 z - 0.2891) * (1.2353 z - 0.8627)

= z^3 + (-1.4156 + 0.5951) z^2 + (0.6090 - 0.7727) z^1 + (-0.0011 + 0.2494)

= z^3 - 0.8205 z^2 - 0.1637 z^1 + 0.2483

AF = conv(A,F);
BG = conv(B,G);

AF =

    1.0000   -1.4156    0.6090   -0.0011
    
 BG =

    0.5951   -0.7727    0.2494   

验证结果无误。

之后便可以讲控制器的实际参数带入差分方程，输出控制器效果波形：

%*****************************************输出控制器效果波形***********************************************%
L=400; %控制步数
uk = zeros(d + nb,1); %输入初值：uk(i)表示u(k-i)
yk = zeros(na,1); %输出初值
yrk = zeros(na,1); %期望输出初值
yr = 10 * [ones(L/4,1);-ones(L/4,1);ones(L/4,1);-ones(L/4+d,1)]; %期望输出

for k=1:L
    time(k)=k*Ts;
    y(k) = -A(2:na+1) * yk + B * uk(d:d+nb); %采集输出数据
    
    u(k)=(-F(2:nf+1) * uk(1:nf) + R * [yr(k+d:-1:k+d-min(d,nr)) ; yrk(1:nr-d)]-G*[y(k) ; yk(1:ng)]) /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);
        yrk(i)=yrk(i-1);
    end
    yk(1)=y(k);
    yrk(1)=yr(k);
end
subplot(2,1,1);
plot(time,yr(1:L),'r:',time,y);
xlabel('t'); ylabel('y_r(t)、y(t)');
legend('y_r(t)','y(t)');
subplot(2,1,2);
plot(time,u);
xlabel('t'); ylabel('u(t)');

以上只是PPC控制的基本步骤，关于PPC控制以及整个控制理论的学习还需要继续深入，比如：期望闭环系统的自然频率与阻尼系数如何选取，是自行设计还是测量得到？如何自动判定哪些零点是可以对消的，哪些是不可以对消的？零极点对于整个控制的影响到底是怎样的？如果期望的闭环传递函数更高阶改怎样设计参数？观测器A多项式的参数如何选取？

这些问题都要更加深入的去学习控制理论，学无止境！希望能给我点个赞！