简单二阶系统

先研究最简单的二阶积分器串联型系统
$$\begin{array}{rl}& {\dot{x}}_1=x_2\\ & {\dot{x}}_2=u\end{array}$$
使用零阶保持法离散化（见附录），
$$A=\left[\begin{array}{cc}0& 1\\ 0& 0\end{array}\right],B=\left[\begin{array}{c}0\\ 1\end{array}\right]$$
设控制器周期（即离散化周期）为 $T$ 已知，则矩阵指数
$$\begin{array}{rl}& (sI-A{)}^{-1}={\left[\begin{array}{cc}s& -1\\ 0& s\end{array}\right]}^{-1}=\frac{1}{s^2}\left[\begin{array}{cc}s& 1\\ 0& s\end{array}\right]\\ & {\text{e}}^{At}={\mathcal{L}}^{-1}(sI-A{)}^{-1}={\mathcal{L}}^{-1}\left[\begin{array}{cc}\frac{1}{s}& \frac{1}{s^2}\\ 0& \frac{1}{s}\end{array}\right]=\left[\begin{array}{cc}1& t\\ 0& 1\end{array}\right]\\ & F={\text{e}}^{AT}=\left[\begin{array}{cc}1& T\\ 0& 1\end{array}\right]\\ & G={\int }_0^T{\text{e}}^{At}\text{d}tB=\left[\begin{array}{cc}T& \frac{T^2}{2}\\ 0& T\end{array}\right]\left[\begin{array}{c}0\\ 1\end{array}\right]=\left[\begin{array}{c}\frac{T^2}{2}\\ T\end{array}\right]\end{array}$$
于是
$$\begin{array}{rl}\vec{x}(n+1)=& F\vec{x}(n)+Gu(n)\\ =& \left[\begin{array}{cc}1& T\\ 0& 1\end{array}\right]\vec{x}(n)+\left[\begin{array}{c}\frac{T^2}{2}\\ T\end{array}\right]u(n)\end{array}$$
将
$$u(n)=-K\vec{x}(n)=-\left[\begin{array}{cc}{k}_1& k_2\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\end{array}\right]$$
代入得
$$\begin{array}{rl}\vec{x}(n+1)=& (F-GK)\vec{x}(n)\\ =& \left(\left[\begin{array}{cc}1& T\\ 0& 1\end{array}\right]+\left[\begin{array}{cc}\frac{T^2}{2}{k}_1& \frac{T^2}{2}{k}_2\\ T{k}_1& T{k}_2\end{array}\right]\right)\vec{x}(n)\end{array}$$
给定两个极点就可以使用比较系数法计算出 $k_1,k_2$。除此以外还有一种便于使用计算机求解的方法
$$K=\left[\begin{array}{cccc}0& \cdots & 0& 1\end{array}\right]\left[\begin{array}{ccc}G& FG& F^{n-1}G\end{array}\right]\alpha (F)$$
其中 $\alpha (x)$ 为给定极点的特征多项式，例如给定两个极点为 ${\lambda }_1,{\lambda }_2$，则 $\alpha (x)=(x-{\lambda }_1)(x-{\lambda }_2)$。设给定重极点 ${\lambda }_1={\lambda }_2=-5$，变换到离散域 $z={\text{e}}^{sT}={\text{e}}^{-0.5}$。仿真结果如图所示，图中 $u$ 为控制量，$y=x_1$ 为状态量。

参考

教材[1]里有更完整的推导过程，链接[2]里有另一个混合系统控制的例子。

1. 孙增圻. 计算机控制理论与应用. 清华大学出版社.
2. simucpp系列教程(3)例程解析(第二部分)

附录：零阶保持器离散化

对使用状态空间描述的线性系统
$$\begin{array}{l}{\vec{x}}^{\prime }=A\vec{x}+Bu\\ y=C\vec{x}+Du\end{array}$$
使用零阶保持法可以离散化成
$$\begin{array}{l}\vec{x}(k+1)=F\vec{x}(k)+Gu(k)\\ y(k)=C\vec{x}(k)+Du(k)\end{array}$$
其中
$$F={\text{e}}^{AT},G={\int }_0^T{\text{e}}^{At}\text{d}tB$$
${\text{e}}^{AT}$ 是矩阵指数。
推导
$$\begin{array}{rl}& x^{\prime }(t)=Ax(t)+Bu(t)\\ & x(t)={\text{e}}^{A(t-t_0)}x(t_0)+{\int }_{t_0}^t{\text{e}}^{A(t-\tau )}Bu(\tau )\text{d}\tau \end{array}$$
令$t_0=kT$，$t=(k+1)T$，则
$$\begin{array}{rl}& x((k+1)T)={\text{e}}^{AT}x(kT)+{\int }_{kT}^{(k+1)T}{\text{e}}^{A(kT+T-\tau )}Bu(\tau )\text{d}\tau \\ & {\int }_{kT}^{(k+1)T}{\text{e}}^{A(kT+T-\tau )}Bu(\tau )\text{d}\tau \stackrel{t=kT+T-\tau }{=}{\int }_0^T{\text{e}}^{At}Bu(kT+T-t)\text{d}t\\ & x(k+1)={\text{e}}^{AT}x(k)+{\int }_0^T{\text{e}}^{At}B\text{d}tu(k)\end{array}$$

附录：代码

计算反馈控制系数 $K$ 的代码

#include <iostream>
#include "zhnmat.hpp"
using namespace zhnmat;
using namespace std;
typedef std::vector<double>  vecdble;
Mat polynomial(Mat x, vecdble lambda) {
    double order = x.row();
    if (order != x.col()) return Mat();
    Mat ans = eye(order);
    for (uint32_t i = 0; i < order; i++)
        ans *= x - lambda[i]*eye(order);
    return x;
}
int main(void) {
    double T = 0.1;
    Mat F(2, 2, vecdble{1, T, 0, 1});
    Mat G(2, 1, vecdble{T*T*0.5, T});
    Mat K(1, 2);
    K.set(0, 0, G.at(1,0));
    K.set(0, 1, (F*G).at(1,0));
    K *= polynomial(F, vecdble{-3, -3});
    cout << K << endl;
}

二阶积分器串联型系统仿真代码

#include "simucpp.hpp"
using namespace simucpp;
int main() {
    Simulator sim1;
    FUIntegrator(uint1, &sim1);
    FUIntegrator(uint2, &sim1);
    FUConstant(cnst1, &sim1);
    FUOutput(uout1, &sim1);
    sim1.connectU(cnst1, uint1);
    sim1.connectU(uint1, uint2);
    sim1.connectU(uint2, uout1);
    sim1.Initialize();
    double ctrl;
    for (uint32_t n1 = 0; n1 < 20; n1++) {
        ctrl = (3-uint2->Get_OutValue()) * 15.24;
        ctrl -= uint1->Get_OutValue() * 7.04;
        for (uint32_t i = 0; i < 100; i++) {
            cnst1->Set_OutValue(ctrl);
            sim1.Simulate_OneStep();
        }
    }
    sim1.Plot();
    return 0;
}