自动控制：控制系统的稳定性

在自动控制领域，控制系统的稳定性是一个至关重要的问题。稳定性决定了系统在受到扰动后是否能够恢复到平衡状态。本文将介绍控制系统稳定性的基本概念、如何利用特征值分析稳定性，并通过具体示例和Python代码展示如何判断系统的稳定性。

控制系统的基本概念

平衡点及其稳定性

平衡点是系统在没有外部扰动时所处的静止状态。根据平衡点的性质，可以将其分为稳定平衡点和不稳定平衡点。

稳定平衡点：系统在初始扰动后最终能够返回其原始位置并保持在那里。例如，一个球位于两座山之间的谷底，当球稍微受到推力上山后会滚回谷底。
不稳定平衡点：任何方向的扰动都会导致系统远离平衡点。例如，一个球位于山顶时，任何方向的轻微扰动都会导致球滚下山。

特征值和系统稳定性

特征值和特征向量是分析系统稳定性的重要工具。通过求解系统矩阵的特征值，我们可以判断系统在平衡点附近的稳定性。

特征值的影响

实部为负：系统稳定，扰动会衰减。
实部为正：系统不稳定，扰动会放大。
实部为零，虚部不为零：系统可能是边界稳定，需要进一步分析。

示例分析：弹簧-阻尼系统

考虑一个简单的弹簧-阻尼系统，其运动方程为：

$\ddot{x} + c \dot{x} + k x = 0$

转换为状态空间模型：

$\mathbf{x}'(t) = \begin{pmatrix} 0 & 1 \\ -\frac{k}{m} & -\frac{c}{m} \end{pmatrix} \mathbf{x}(t)$

通过计算系统矩阵 $A$ 的特征值来分析其稳定性。

Python实现

以下代码演示了如何使用Python计算特征值并判断系统稳定性：

import numpy as np
import scipy.linalg as la
import matplotlib.pyplot as plt# 参数
m = 1.0  # 质量
c = 0.5  # 阻尼系数
k = 2.0  # 弹簧系数# 系统矩阵A
A = np.array([[0, 1], [-k/m, -c/m]])# 计算特征值
eigvals, eigvecs = la.eig(A)print("特征值:")
print(eigvals)# 判断稳定性
stable = all(eig.real < 0 for eig in eigvals)
print("系统是否渐近稳定:", stable)# 可视化
time = np.linspace(0, 10, 1000)
x0 = [1, 0]  # 初始条件def state_space_solution(A, eigvals, eigvecs, x0, t):c = la.solve(eigvecs, x0)x_t = np.zeros((len(t), len(x0)), dtype=complex)for i in range(len(eigvals)):x_t += c[i] * np.exp(eigvals[i] * t)[:, np.newaxis] * eigvecs[:, i]return np.real(x_t)solution = state_space_solution(A, eigvals, eigvecs, x0, time)plt.plot(time, solution[:, 0], label='x1(t)')
plt.plot(time, solution[:, 1], label='x2(t)')
plt.xlabel('Time t')
plt.ylabel('Solution x(t)')
plt.legend()
plt.title('Solution of ODE using Eigenvalue Method')
plt.grid(True)
plt.show()

在这里插入图片描述

代码说明

定义系统参数：定义弹簧-阻尼系统的参数 $m$ 、 $c$ 、 $k$ 。
系统矩阵：构造状态空间矩阵 $A$ 。
计算特征值和特征向量：使用 scipy.linalg.eig 计算特征值和特征向量。
判断稳定性：根据特征值的实部判断系统是否渐近稳定。
可视化解：通过初始条件计算系统的时间响应，并绘制解的曲线。

虚数特征值的影响

虚数或复数特征值会导致系统的响应表现为振荡。具体而言：

正实部：系统不稳定，振荡的幅度会随时间增加。
负实部：系统稳定，振荡的幅度会随时间减小。

示例分析：正实部的复特征值

考虑一个系统，其特征值为 $\lambda = \alpha \pm j\beta$ ，其中 $\alpha > 0$ 。这种情况下，系统会表现出不稳定振荡，其响应随时间增加，振荡的幅度不断扩大。

import numpy as np
import matplotlib.pyplot as plt# 参数
alpha = 0.1
beta = 1.0# 时间序列
time = np.linspace(0, 50, 1000)# 解的计算
x1 = np.exp(alpha * time) * np.cos(beta * time)
x2 = np.exp(alpha * time) * np.sin(beta * time)plt.plot(time, x1, label='x1(t)')
plt.plot(time, x2, label='x2(t)')
plt.xlabel('Time t')
plt.ylabel('Solution x(t)')
plt.legend()
plt.title('Response of System with Complex Eigenvalues (Unstable)')
plt.grid(True)
plt.show()