状态方程的离散化仿真

1.离散化模型的推导

设系统的状态空间表达式为:
x ˙ ( t ) = A x + B u ( t ) y ( t ) = C x ( t ) + D u ( t ) \dot{\mathbf{x}}(t)=A\mathbf{x}+B\mathbf{u}(t)\\ \mathbf{y}(t)=C\mathbf{x}(t)+D\mathbf{u}(t) x˙(t)=Ax+Bu(t)y(t)=Cx(t)+Du(t)
对上面两个式子分别取拉普拉斯变换得到:
X ( s ) = ( s I − A ) − 1 X ( 0 ) + ( s I − A ) − 1 B U ( s ) X(s)=(sI-A)^{-1}X(0)+(sI-A)^{-1}BU(s) X(s)=(sI−A)−1X(0)+(sI−A)−1BU(s)
再对上面的式子求拉普拉斯反变换得到:
X ( t ) = e A t X ( 0 ) + ∫ 0 t e A ( t − τ ) B u ( τ ) d τ = F ( t ) X ( 0 ) + ∫ 0 t F ( t − τ ) B u ( τ ) d τ F ( t ) = e A t \mathbf{X}(t)=e^{At}\mathbf{X}(0)+\int_0^te^{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} X(t)=eAtX(0)+∫0teA(t−τ)Bu(τ)dτ=F(t)X(0)+∫0tF(t−τ)Bu(τ)dτF(t)=eAt
对上面的结果进行离散化处理,设采样周期是 T T T，考察 t = ( k + 1 ) T , t = k T t=(k+1)T,t=kT t=(k+1)T,t=kT时 X ( t ) \mathbf{X}(t) X(t)的值。
X ( ( k + 1 ) T ) = e A ( k + 1 ) T X ( 0 ) + ∫ 0 ( k + 1 ) T e A ( ( k + 1 ) T − τ ) B u ( τ ) d τ = e A T ( e A k T X ( 0 ) + ∫ 0 k T e A ( k T − τ ) B u ( τ ) d τ + ∫ k T k ( T + 1 ) e A ( k T − τ ) B u ( τ ) d τ ) = e A T ( X ( k T ) + ∫ k T k ( T + 1 ) e A ( k T − τ ) B u ( τ ) d τ ) \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) X((k+1)T)=eA(k+1)TX(0)+∫0(k+1)TeA((k+1)T−τ)Bu(τ)dτ=eAT(eAkTX(0)+∫0kTeA(kT−τ)Bu(τ)dτ+∫kTk(T+1)eA(kT−τ)Bu(τ)dτ)=eAT(X(kT)+∫kTk(T+1)eA(kT−τ)Bu(τ)dτ)
若离散化时对输入 u u u采用零阶保持器,则有:
u ( t ) = u ( k T ) , ( k T ≤ t < ( k + 1 ) T ) u(t)=u(kT),(kT\leq t <(k+1)T) u(t)=u(kT),(kT≤t<(k+1)T)
于是上式可以变化为:
X ( ( k + 1 ) T ) = e A T X ( k T ) + ∫ k T ( k + 1 ) T e A ( ( k + 1 ) T − τ ) B u ( τ ) d τ = e A T X ( k T ) + ∫ k T ( k + 1 ) T e A ( ( k + 1 ) T − τ ) B d τ u ( k T ) = e A T X ( k T ) + ∫ 0 T e A ( T − τ ) B d τ u ( k T ) \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) X((k+1)T)=eATX(kT)+∫kT(k+1)TeA((k+1)T−τ)Bu(τ)dτ=eATX(kT)+∫kT(k+1)TeA((k+1)T−τ)Bdτu(kT)=eATX(kT)+∫0TeA(T−τ)Bdτu(kT)
写成一般形式如下:
X ( k + 1 ) = F ( T ) X ( k ) + G ( T ) u ( k ) \mathbf{X}(k+1)=\mathbf{F}(T)\mathbf{X}(k)+\mathbf{G}(T)\mathbf{u}(k) X(k+1)=F(T)X(k)+G(T)u(k)
其中 F ( T ) = e A T \mathbf{F}(T)=e^{AT} F(T)=eAT, G ( T ) = ∫ 0 T e A ( T − τ ) B d τ = ∫ 0 T F ( T − τ ) B d τ \mathbf{G}(T)=\int_{0}^{T}e^{A(T-\tau)}Bd\tau=\int_{0}^{T}F(T-\tau)Bd\tau G(T)=∫0TeA(T−τ)Bdτ=∫0TF(T−τ)Bdτ。

这样就得到采用零阶保持器的离散化模型:
X ( k + 1 ) = F ( T ) X ( k ) + G ( T ) u ( k ) Y ( k ) = C X ( k ) + D u ( k ) \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) X(k+1)=F(T)X(k)+G(T)u(k)Y(k)=CX(k)+Du(k)
若采用一阶保持器的离散化模型:
u ( t ) = u ( k T ) + u ˙ ( k T ) ( t − k T ) , ( k T ≤ t < ( k + 1 ) T ) u(t)=u(kT)+\dot{u}(kT)(t-kT),(kT\leq t <(k+1)T) u(t)=u(kT)+u˙(kT)(t−kT),(kT≤t<(k+1)T)
带入之前的模型得到:
X ( k + 1 ) = F ( T ) X ( k ) + G ( T ) u ( k ) + G 1 ( T ) u ˙ ( k ) \mathbf{X}(k+1)=\mathbf{F}(T)\mathbf{X}(k)+\mathbf{G}(T)\mathbf{u}(k)+\mathbf{G_1}(T)\dot{\mathbf{u}}(k) X(k+1)=F(T)X(k)+G(T)u(k)+G1(T)u˙(k)
其中 F ( T ) = e A T \mathbf{F}(T)=e^{AT} F(T)=eAT, G ( T ) = ∫ 0 T e A ( T − τ ) B d τ = ∫ 0 T F ( T − τ ) B d τ \mathbf{G}(T)=\int_{0}^{T}e^{A(T-\tau)}Bd\tau=\int_{0}^{T}F(T-\tau)Bd\tau G(T)=∫0TeA(T−τ)Bdτ=∫0TF(T−τ)Bdτ, G 1 ( T ) = ∫ 0 T F ( T − τ ) B τ d τ \mathbf{G_1}(T)=\int_{0}^{T}F(T-\tau)B\tau d\tau G1(T)=∫0TF(T−τ)Bτdτ。

2.重要矩阵的计算

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

2.1 F ( T ) \mathbf{F}(T) 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∑Li!(AT)i=I+AT{I+2AT{I+...+L−1AT(I+LAT)}}

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

G ( T ) = ∫ 0 T e A λ B d λ = ∫ 0 T ∑ i = 0 L ( A λ ) i i ! B d λ = ∑ i = 0 L ∫ 0 T ( A λ ) i i ! B d λ = ∑ i = 0 L A i T i + 1 ( i + 1 ) ! B = A − 1 ( ∑ j = 1 L + 1 ( A T ) j j ! ) B = A − 1 ( F ( t ) − I ) B \mathbf{G}(T)=\int_0^Te^{\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^iT^{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} G(T)=∫0TeAλBdλ=∫0Ti=0∑Li!(Aλ)iBdλ=i=0∑L∫0Ti!(Aλ)iBdλ=i=0∑L(i+1)!AiTi+1B=A−1(j=1∑L+1j!(AT)j)B=A−1(F(t)−I)B

3.仿真实验

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

将系统转化为以下形式:

那么就有以下的状态空间方程:
x ˙ = A x + B u y = ( 0 1 ) x \dot{\mathbf{x}}=A\mathbf{x}+B\mathbf{u}\\ y=(0 \quad 1)\mathbf{x} x˙=Ax+Buy=(01)x
其中 x = ( x 1 x 2 ) T \mathbf{x}=(x_1\quad x_2)^T x=(x1x2)T, A = ( 0 0 1 − 1 ) A=\begin{pmatrix}0&0\\1&-1 \end{pmatrix} A=(010−1), B = ( 8 0 ) B=\begin{pmatrix}8\\0 \end{pmatrix} B=(80), r = s i n t r=sint r=sint。

若设采样时间 T = 0.1 T=0.1 T=0.1,可以得到状态空间方程的零阶保持器的离散化模型:
F ( T ) = ( 1 0 1 − e − T e − T ) , G ( T ) = ( 8 T 8 ( T − 1 + e − T ) ) \mathbf{F}(T)=\begin{pmatrix}1&0\\1-e^{-T}&e^{-T} \end{pmatrix},\mathbf{G}(T)=\begin{pmatrix}8T\\8(T-1+e^{-T}) \end{pmatrix} F(T)=(11−e−T0e−T),G(T)=(8T8(T−1+e−T))
于是离散化模型为:
( x 1 ( k + 1 ) x 2 ( k + 1 ) ) = ( 1 0 1 − e − T e − T ) ( x 1 ( k ) x 2 ( k ) ) + ( 8 T 8 ( T − 1 + e − T ) ) u ( k ) \begin{pmatrix}x_1(k+1)\\x_2(k+1)\end{pmatrix}=\begin{pmatrix}1&0\\1-e^{-T}&e^{-T} \end{pmatrix}\begin{pmatrix}x_1(k)\\x_2(k)\end{pmatrix}+\begin{pmatrix}8T\\8(T-1+e^{-T}) \end{pmatrix}u(k) (x1(k+1)x2(k+1))=(11−e−T0e−T)(x1(k)x2(k))+(8T8(T−1+e−T))u(k)

y ( k ) = x 2 ( k ) y(k)=x_2(k) y(k)=x2(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:nume = r(i) - y;if abs(e)>=au = a;elseu = e;endx = 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)]);

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