控制系统/线性系统解决方案

状态方程解

本章的解决方案深深扎根于对 常微分方程 的先验知识。读者在阅读本章之前应该先了解该主题。

状态方程是一个一阶线性微分方程，或更准确地说，是一个线性微分方程组。由于这是一个一阶方程，我们可以使用 常微分方程 的结果来找到关于状态变量x 的方程的通解。一旦状态方程对x 求解，该解就可以代入输出方程。得到的方程将显示系统输入和系统输出之间的直接关系，无需明确考虑系统的内部状态。本章中的各节将讨论状态空间方程的解，从最简单的案例（时不变、无输入）开始，到最复杂的案例（时变系统）结束。

求解无输入时的 x(t)

再次看看状态方程

$x'=Ax(t)+Bu(t)$

我们可以看到，这个方程是一个一阶微分方程，只是变量是向量，系数是矩阵。然而，由于矩阵微积分的规则，这些区别并不重要。我们可以忽略输入项（现在），并将该方程改写为以下形式

${\frac {dx(t)}{dt}}=Ax(t)$

我们可以像这样分离变量

${\frac {dx(t)}{x(t)}}=Adt$

对两边积分，并将两边都乘以e 的幂，我们得到结果

$x(t)=e^{At+C}$

其中C 是一个常数。我们可以将D = e^C 赋值以使方程更容易，但我们也知道，D 然后将成为系统的初始条件。如果我们将变量t 中的值零代入，这将变得很明显。该方程的最终解如下所示

$x(t)=e^{A(t-t_{0})}x(t_{0})$

我们将矩阵指数e^At 称为状态转移矩阵，计算它虽然有时很困难，但对于分析和操纵系统至关重要。我们将在下面详细讨论计算矩阵指数。

求解有非零输入时的 x(t)

然而，如果我们的输入是非零的（通常情况下，任何有趣的系统都是这样），我们的解决方案会稍微复杂一些。请注意，现在我们有了输入项，我们不能再轻松地分离变量并对两边积分。

$x'(t)=Ax(t)+Bu(t)$

我们减去以获得 $Ax(t)$ 在左侧，然后我们做一些奇怪的事情；我们用逆状态转移矩阵在两边预乘

$e^{-At}x'(t)-e^{-At}Ax(t)=e^{-At}Bu(t)$

这最后一步的理由可能看起来很模糊，所以我们将用一个例子来说明这一点

示例

对以下表达式关于时间求导

$e^{-At}x(t)$

微分中的乘积法则提醒我们，如果我们有两个函数相乘

$f(t)g(t)$

然后我们对 t 求导，结果为

$f(t)g'(t)+f'(t)g(t)$

如果我们相应地设置我们的函数

$f(t)=e^{-At}\qquad f'(t)=-Ae^{-At}$

$g(t)=x(t)\qquad g'(t)=x'(t)$

则输出结果为

$e^{-At}x'(t)-e^{-At}Ax(t)$

如果我们看一下这个结果，它与我们上面方程的结果相同。

使用我们例子中的结果，我们可以将我们方程左侧浓缩成一个导数

${\frac {d(e^{-At}x(t))}{dt}}=e^{-At}Bu(t)$

现在我们可以对两边积分，从初始时间 (t₀) 到当前时间 (t)，使用一个哑变量 τ，我们将更接近我们的结果。最后，如果我们预乘以 e^At，我们将得到最终结果

[一般状态方程解]

$x(t)=e^{A(t-t_{0})}x(t_{0})+\int _{t_{0}}^{t}e^{A(t-\tau )}Bu(\tau )d\tau$

如果我们将这个解代入输出方程，我们将得到

[一般输出方程解]

$y(t)=Ce^{A(t-t_{0})}x(t_{0})+C\int _{t_{0}}^{t}e^{A(t-\tau )}Bu(\tau )d\tau +Du(t)$

这是状态空间方程的非零输入的广义时不变解。这些方程是重要的结果，对控制系统有进一步研究兴趣的学生应该记住这些方程。

状态转移矩阵

有关矩阵指数的更多信息，请参阅
工程分析

状态转移矩阵，e^At，是上述时不变情况下广义状态空间解的重要组成部分。计算这个矩阵指数函数是分析新系统时首先要做的几件事之一，计算结果将揭示有关所讨论系统的关键信息。

可以使用泰勒级数展开直接计算矩阵指数

$e^{At}=\sum _{n=0}^{\infty }{\frac {(At)^{n}}{n!}}$

有关对角矩阵和约旦型矩阵的更多信息，请参阅
工程分析

此外，我们可以尝试将矩阵 A 对角化为对角矩阵或约旦标准型矩阵。对角矩阵的指数只是对角元素分别乘以该指数。约旦标准型矩阵的指数稍微复杂一些，但存在一个有用的模式可以利用来快速找到解。有兴趣的读者应该阅读工程分析中的相关段落。

状态转移矩阵，以及一般意义上的矩阵指数，是控制工程中非常重要的工具。

对角矩阵

如果矩阵是对角矩阵，则状态转移矩阵可以通过将矩阵的每个对角元素作为e的幂次方来计算。

约旦标准型

如果A矩阵是约旦标准型，则可以使用以下公式快速生成矩阵指数

$e^{Jt}=e^{\lambda t}{\begin{bmatrix}1&t&{\frac {1}{2!}}t^{2}&\cdots &{\frac {1}{n!}}t^{n}\\0&1&t&\cdots &{\frac {1}{(n-1)!}}t^{n-1}\\\vdots &\vdots &\vdots &\ddots &\vdots \\0&0&0&\cdots &1\end{bmatrix}}$

其中λ是约旦标准型矩阵的特征值（对角线上的值）。

拉普拉斯反变换方法

我们可以通过对以下拉普拉斯反变换求解来计算状态转移矩阵（或任何矩阵指数函数）

$e^{At}={\mathcal {L}}^{-1}[(sI-A)^{-1}]$

如果A是高阶矩阵，则求逆运算可能难以解决。

如果A矩阵是约旦标准型，则可以使用以下公式快速生成矩阵指数

   $e^{Jt}=e^{\lambda t}{\begin{bmatrix}1&t&{\frac {1}{2!}}t^{2}&\cdots &{\frac {1}{n!}}t^{n}\\0&1&t&\cdots &{\frac {1}{(n-1)!}}t^{n-1}\\\vdots &\vdots &\vdots &\ddots &\vdots \\0&0&0&\cdots &1\end{bmatrix}}$

其中λ是约旦标准型矩阵的特征值（对角线上的值）。

谱分解

如果我们知道A的所有特征值，我们可以创建我们的转移矩阵T和逆转移矩阵T^-1。这些矩阵将分别是右特征向量和左特征向量的矩阵。如果我们同时拥有左特征向量和右特征向量，我们可以计算出状态转移矩阵为

[谱分解]

$e^{At}=\sum _{i=1}^{n}e^{\lambda _{i}t}v_{i}w_{i}'$

请注意，w_i' 是第i个左特征向量的转置，而不是它的导数。我们将在后面的章节中详细讨论“特征值”、“特征向量”以及谱分解技术。

凯莱-哈密顿定理

有关凯莱-哈密顿定理的更多信息，请参阅
工程分析

凯莱-哈密顿定理也可以用来找到矩阵指数的解。对于系统矩阵A的任何特征值λ，我们可以证明这两个方程是等价的

$e^{\lambda t}=a_{0}+a_{1}\lambda t+a_{2}\lambda ^{2}t^{2}+\cdots +a_{n-1}\lambda ^{n-1}t^{n-1}$

一旦我们解出方程的系数a，就可以将这些系数代入以下方程

$e^{At}=a_{0}I+a_{1}At+a_{2}A^{2}t^{2}+\cdots +a_{n-1}A^{n-1}t^{n-1}$

示例：非对角矩阵

给定以下矩阵A，求状态转移矩阵

$A={\begin{bmatrix}0&1\\-1&0\end{bmatrix}}$

我们可以找到该矩阵的特征值为λ = i，-i。如果我们将这些值代入我们的特征向量方程，得到

${\begin{vmatrix}i&-1\\1&i\end{vmatrix}}v_{1}=0$

${\begin{vmatrix}-i&-1\\1&-i\end{vmatrix}}v_{2}=0$

我们可以解出特征向量

$v_{1}={\begin{bmatrix}1\\i\end{bmatrix}}$

$v_{2}={\begin{bmatrix}1\\-i\end{bmatrix}}$

有了特征向量，我们就可以解出左特征向量

$w_{1}={\begin{bmatrix}1\\-i\end{bmatrix}}$

$w_{2}={\begin{bmatrix}1\\i\end{bmatrix}}$

现在，利用谱分解，我们可以构建状态转移矩阵

$e^{At}=e^{it}{\begin{bmatrix}1\\i\end{bmatrix}}{\begin{bmatrix}1&-i\end{bmatrix}}+e^{-it}{\begin{bmatrix}1\\-i\end{bmatrix}}{\begin{bmatrix}1&i\end{bmatrix}}$

如果我们记得欧拉恒等式，我们可以将复指数分解成正弦函数。执行向量相乘后，所有的虚数项都会抵消，剩下的就是我们的结果

$e^{At}={\begin{bmatrix}\cos t&\sin t\\-\sin t&\cos t\end{bmatrix}}$

鼓励读者进行乘法运算，并尝试推导出该结果。

示例：Sympy 计算

使用免费的 Python 库 'sympy'，我们可以非常容易地自动计算状态转移矩阵

>>> from sympy import *
>>> t = symbols('t', positive = true)
>>> A = Matrix([[0,1],[-1,0]])
>>> exp(A*t).expand(complex=True)

⎡cos(t)   sin(t)⎤
⎢               ⎥
⎣-sin(t)  cos(t)⎦

你也可以在这个网站上尝试一下

sympy live

示例：MATLAB 计算

使用 MATLAB 中的符号工具箱，我们可以编写 MATLAB 代码来自动生成给定输入矩阵A的状态转移矩阵。以下是一个可以执行此任务的 MATLAB 代码示例

function [phi] = statetrans(A)
   t = sym('t');
   phi = expm(A * t);
end

使用此 MATLAB 函数来查找以下矩阵的状态转移矩阵（警告，计算可能需要一些时间）

1. $A_{1}={\begin{bmatrix}2&0\\0&2\end{bmatrix}}$
2. $A_{2}={\begin{bmatrix}0&1\\-1&0\end{bmatrix}}$
3. $A_{3}={\begin{bmatrix}2&1\\0&2\end{bmatrix}}$

矩阵 1 是一个对角矩阵，矩阵 2 具有复特征值，矩阵 3 是约旦标准型。这三个矩阵应该代表系统矩阵的一些常见形式。以下代码片段是输入到 MATLAB 中以生成这些矩阵的输入命令以及输出结果

矩阵 A1

>> A1 = [2 0 ; 0 2];
>> statetrans(A1)
 
ans =
 
[ exp(2*t),        0]
[        0, exp(2*t)]

矩阵 A2

>> A2 = [0 1 ; -1 0];
>> statetrans(A1)
 
ans =
 
[  cos(t),  sin(t)]
[ -sin(t),  cos(t)]

矩阵 A3

>> A1 = [2 1 ; 0 2];
>> statetrans(A1)
 
ans =
 
[   exp(2*t), t*exp(2*t)]
[          0,   exp(2*t)]

示例：MATLAB 中的多种方法

在 MATLAB 中，有多种方法可以从标量（时不变）矩阵 A 计算状态转移矩阵。以下方法都将依赖于 *符号工具箱* 来执行方程操作。在每个代码片段的末尾，变量 **eAt** 包含矩阵 A 的状态转移矩阵。

直接法

t = sym('t');
eAt = expm(A * t);

拉普拉斯变换法

s = sym('s');
[n,n] = size(A);
in = inv(s*eye(n) - A);
eAt = ilaplace(in);

谱分解

t = sym('t');
[n,n] = size(A);
[V, e] = eig(A);
W = inv(V);
sum = [0 0;0 0];
for I = 1:n
   sum = sum + expm(e(I,I)*t)*V(:,I)*W(I,:);
end;
eAt = sum;

这三种方法都应该产生相同的答案。鼓励学生验证这一点。

← 状态空间方程

控制系统

时变系统解 →