状态方程的离散化仿真

1.离散化模型的推导

设系统的状态空间表达式为:
$$\begin{gathered} \dot{\mathbf{x}}(t)=A\mathbf{x}+B\mathbf{u}(t) \\ \mathbf{y}(t)=C\mathbf{x}(t)+D\mathbf{u}(t) \end{gathered}$$
对上面两个式子分别取拉普拉斯变换得到:
$$X(s)=(sI-A{)}^{-1}X(0)+(sI-A{)}^{-1}BU(s)$$
再对上面的式子求拉普拉斯反变换得到:
$$\begin{gathered} \mathbf{X}(t)=e^{At}\mathbf{X}(0)+{\int }_0^t{e}^{A(t-\tau )}B\mathbf{u}(\tau )d\tau \\ =\mathbf{F}(t)\mathbf{X}(0)+{\int }_0^t\mathbf{F}(t-\tau )B\mathbf{u}(\tau )d\tau \\ \mathbf{F}(t)=e^{At} \end{gathered}$$
对上面的结果进行离散化处理,设采样周期是$T$，考察$t=(k+1)T,t=kT$时$\mathbf{X}(t)$的值。
$$\begin{gathered} \mathbf{X}((k+1)T)=e^{A(k+1)T}\mathbf{X}(0)+{\int }_0^{(k+1)T}{e}^{A((k+1)T-\tau )}B\mathbf{u}(\tau )d\tau \\ =e^{AT}(e^{AkT}\mathbf{X}(0)+{\int }_0^{kT}{e}^{A(kT-\tau )}B\mathbf{u}(\tau )d\tau +{\int }_{kT}^{k(T+1)}{e}^{A(kT-\tau )}B\mathbf{u}(\tau )d\tau ) \\ =e^{AT}(\mathbf{X}(kT)+{\int }_{kT}^{k(T+1)}{e}^{A(kT-\tau )}B\mathbf{u}(\tau )d\tau ) \end{gathered}$$
若离散化时对输入$u$采用零阶保持器,则有:
$$u(t)=u(kT),(kT\le t<(k+1)T)$$
于是上式可以变化为:
$$\begin{gathered} \mathbf{X}((k+1)T)=e^{AT}\mathbf{X}(kT)+{\int }_{kT}^{(k+1)T}{e}^{A((k+1)T-\tau )}B\mathbf{u}(\tau )d\tau \\ =e^{AT}\mathbf{X}(kT)+{\int }_{kT}^{(k+1)T}{e}^{A((k+1)T-\tau )}Bd\tau \mathbf{u}(kT) \\ =e^{AT}\mathbf{X}(kT)+{\int }_0^T{e}^{A(T-\tau )}Bd\tau \mathbf{u}(kT) \end{gathered}$$
写成一般形式如下:
$$\mathbf{X}(k+1)=\mathbf{F}(T)\mathbf{X}(k)+\mathbf{G}(T)\mathbf{u}(k)$$
其中$\mathbf{F}(T)=e^{AT}$,$\mathbf{G}(T)={\int }_0^T{e}^{A(T-\tau )}Bd\tau ={\int }_0^{T}F(T-\tau )Bd\tau$。

这样就得到采用零阶保持器的离散化模型:
$$\begin{gathered} \mathbf{X}(k+1)=\mathbf{F}(T)\mathbf{X}(k)+\mathbf{G}(T)\mathbf{u}(k) \\ \mathbf{Y}(k)=\mathbf{C}\mathbf{X}(k)+\mathbf{D}\mathbf{u}(k) \end{gathered}$$
若采用一阶保持器的离散化模型:
$$u(t)=u(kT)+\dot{u}(kT)(t-kT),(kT\le t<(k+1)T)$$
带入之前的模型得到:
$$\mathbf{X}(k+1)=\mathbf{F}(T)\mathbf{X}(k)+\mathbf{G}(T)\mathbf{u}(k)+{\mathbf{G}}_1(T)\dot{\mathbf{u}}(k)$$
其中$\mathbf{F}(T)=e^{AT}$,$\mathbf{G}(T)={\int }_0^T{e}^{A(T-\tau )}Bd\tau ={\int }_0^{T}F(T-\tau )Bd\tau$,${\mathbf{G}}_1(T)={\int }_0^{T}F(T-\tau )B\tau d\tau$。

2.重要矩阵的计算

目的在于用数值计算$\mathbf{G}(T),\mathbf{F}(T)$。

2.1 $\mathbf{F}(T)$的计算

F ( T ) = e A T ≈ ∑ i = 0 L ( A T ) i i ! = I + A T { I + A T 2 { I + . . . + A T L − 1 ( I + A T L ) } } \mathbf{F}(T)=e^{AT}\approx\sum_{i=0}^L\frac{(AT)^i}{i!}=I+AT\{I+\frac{AT}{2}\{I+...+\frac{AT}{L-1}(I+\frac{AT}{L})\}\} F(T)=eAT≈i=0∑L​i!(AT)i​=I+AT{I+2AT​{I+...+L−1AT​(I+LAT​)}}

2.2 $\mathbf{G}(T)$的计算

$$\begin{gathered} \mathbf{G}(T)={\int }_0^T{e}^{\mathbf{A}\lambda }\mathbf{B}d\lambda ={\int }_0^T\sum _{i=0}^L\frac{(A\lambda {)}^i}{i!}\mathbf{B}d\lambda =\sum _{i=0}^L{\int }_0^T\frac{(A\lambda {)}^i}{i!}\mathbf{B}d\lambda =\sum _{i=0}^L\frac{A^i{T}^{i+1}}{(i+1)!}\mathbf{B} \\ ={\mathbf{A}}^{-1}(\sum _{j=1}^{L+1}\frac{(AT{)}^j}{j!})\mathbf{B}={\mathbf{A}}^{-1}(\mathbf{F}(t)-\mathbf{I})\mathbf{B} \end{gathered}$$

3.仿真实验

例如仿真以下的控制系统(用Simulink搭建的):

将系统转化为以下形式:

那么就有以下的状态空间方程:
$$\begin{gathered} \dot{\mathbf{x}}=A\mathbf{x}+B\mathbf{u} \\ y=(01)\mathbf{x} \end{gathered}$$
其中$\mathbf{x}=(x_1{x}_2{)}^T$,$A=\left(\begin{array}{cc}0& 0\\ 1& -1\end{array}\right)$,$B=\left(\begin{array}{c}8\\ 0\end{array}\right)$,$r=sint$。

若设采样时间$T=0.1$,可以得到状态空间方程的零阶保持器的离散化模型:
$$\mathbf{F}(T)=\left(\begin{array}{cc}1& 0\\ 1-e^{-T}& e^{-T}\end{array}\right),\mathbf{G}(T)=\left(\begin{array}{c}8T\\ 8(T-1+e^{-T})\end{array}\right)$$
于是离散化模型为:
$$\left(\begin{array}{c}{x}_1(k+1)\\ x_2(k+1)\end{array}\right)=\left(\begin{array}{cc}1& 0\\ 1-e^{-T}& e^{-T}\end{array}\right)\left(\begin{array}{c}{x}_1(k)\\ x_2(k)\end{array}\right)+\left(\begin{array}{c}8T\\ 8(T-1+e^{-T})\end{array}\right)u(k)$$

$$y(k)=x_2(k)$$

clc,clear;close all;
%设置初始参数值
K = 8;a = 5;
%设置采样时间T
T = 0.1; t0 = 0;tf = 10;
t = t0:T:tf;x0 = [0 0]';
y0 = 0;
%设置输入r
r = sin(t);
%设置矩阵
F = [1 0;1-exp(-T) exp(-T)];
G = [K*T;K*(T - 1 + exp(-T))];
C = [0 1];
%设置状态变量和输出变量
x = x0;y = y0;
time = t0;num = size(t,2);
y_output = [];t_output = [];
%下面开始循环操作
for i = 1:num
   e = r(i) - y;
   if abs(e)>=a
       u = a;
   else
       u = e;
   end
   x = F*x + G*u;
   y = C*x;
   time = time + T;
   t_output = [t_output,time];
   y_output = [y_output,y];
end
%下面开始画图
hold on;box on;grid on;
plot(t_output,y_output,'bo-');
xlabel('t/s');ylabel('y');
title(['采样时间T:',num2str(T)]);

可以预见到离散系统仿真的效果是很好的。