《自动控制原理》 - 第五章线性系统的频域分析法

啊阿狸不会拉杆

发布于 2026-01-21 11:41:34

2300

5-1 频率特性

频率特性是线性系统在正弦信号作用下的稳态响应特性。对于一个稳定的线性定常系统，当输入为正弦信号时，其稳态输出也为同频率的正弦信号，但幅值和相位会发生变化。幅值的变化倍数和相位的变化量都是输入信号频率的函数，分别称为幅频特性和相频特性，二者合称为系统的频率特性。

频率特性的数学表示

Python 实现示例

下面我们使用 Python 计算并绘制一个系统的频率特性：

import control
import matplotlib.pyplot as plt
import numpy as np

# 设置中文显示
plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号

# 定义系统传递函数
# G(s) = 1 / (s^2 + 2s + 1)
num = [1]  # 分子多项式系数
den = [1, 2, 1]  # 分母多项式系数
sys = control.TransferFunction(num, den)

# 计算频率特性
omega = np.logspace(-2, 2, 1000)  # 频率范围从0.01到100 rad/s
mag, phase, omega = control.bode(sys, omega, dB=True, deg=True, plot=False)

# 绘制幅频特性和相频特性
plt.figure(figsize=(10, 8))

# 幅频特性图
plt.subplot(2, 1, 1)
plt.semilogx(omega, 20 * np.log10(mag))  # 将幅值转换为dB
plt.title('幅频特性')
plt.xlabel('频率 (rad/s)')
plt.ylabel('幅值 (dB)')
plt.grid(True, which='both', linestyle='--', alpha=0.7)

# 相频特性图
plt.subplot(2, 1, 2)
plt.semilogx(omega, phase)  # 相位已经是角度
plt.title('相频特性')
plt.xlabel('频率 (rad/s)')
plt.ylabel('相位 (度)')
plt.grid(True, which='both', linestyle='--', alpha=0.7)

plt.tight_layout()
plt.show()

# 计算指定频率点的频率特性
omega_points = [0.1, 1, 10]
for w in omega_points:
    response = control.evalfr(sys, 1j * w)
    mag = np.abs(response)
    phase = np.angle(response, deg=True)
    print(f"当ω={w} rad/s时：")
    print(f"  幅值: {mag:.4f} ({20 * np.log10(mag):.4f} dB)")
    print(f"  相位: {phase:.4f} 度")

这个示例展示了如何使用 Python 的control库计算并绘制系统的频率特性。我们定义了一个二阶系统，然后计算了其在不同频率下的幅值和相位，并绘制了 Bode 图。

5-2 典型环节与开环系统的频率特性

控制系统通常由多个典型环节组成，了解这些典型环节的频率特性有助于我们分析和设计复杂的控制系统。

典型环节的频率特性

Python 实现示例

下面我们使用 Python 绘制几种典型环节的 Bode 图：

import control
import matplotlib.pyplot as plt
import numpy as np

# 设置中文显示
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams["axes.unicode_minus"] = False  # 解决负号显示问题

# 定义频率范围
omega = np.logspace(-2, 2, 1000)  # 频率范围从0.01到100 rad/s

# 典型环节的传递函数
# 比例环节
K = 2
sys_proportional = control.TransferFunction([K], [1])

# 积分环节
sys_integral = control.TransferFunction([1], [1, 0])

# 微分环节
sys_derivative = control.TransferFunction([1, 0], [1])

# 惯性环节
T = 1
sys_first_order = control.TransferFunction([1], [T, 1])

# 振荡环节
xi = 0.3  # 阻尼比
sys_second_order = control.TransferFunction([1], [T**2, 2*xi*T, 1])

# 绘制Bode图
plt.figure(figsize=(12, 10))

# 比例环节
plt.subplot(5, 2, 1)
mag, phase, _ = control.bode(sys_proportional, omega, dB=True, deg=True, plot=False)
plt.semilogx(omega, 20 * np.log10(mag))
plt.title('比例环节 幅频特性')
plt.grid(True, which='both')

plt.subplot(5, 2, 2)
plt.semilogx(omega, phase)
plt.title('比例环节 相频特性')
plt.grid(True, which='both')

# 积分环节
plt.subplot(5, 2, 3)
mag, phase, _ = control.bode(sys_integral, omega, dB=True, deg=True, plot=False)
plt.semilogx(omega, 20 * np.log10(mag))
plt.title('积分环节 幅频特性')
plt.grid(True, which='both')

plt.subplot(5, 2, 4)
plt.semilogx(omega, phase)
plt.title('积分环节 相频特性')
plt.grid(True, which='both')

# 微分环节
plt.subplot(5, 2, 5)
mag, phase, _ = control.bode(sys_derivative, omega, dB=True, deg=True, plot=False)
plt.semilogx(omega, 20 * np.log10(mag))
plt.title('微分环节 幅频特性')
plt.grid(True, which='both')

plt.subplot(5, 2, 6)
plt.semilogx(omega, phase)
plt.title('微分环节 相频特性')
plt.grid(True, which='both')

# 惯性环节
plt.subplot(5, 2, 7)
mag, phase, _ = control.bode(sys_first_order, omega, dB=True, deg=True, plot=False)
plt.semilogx(omega, 20 * np.log10(mag))
plt.title('惯性环节 幅频特性')
plt.grid(True, which='both')

plt.subplot(5, 2, 8)
plt.semilogx(omega, phase)
plt.title('惯性环节 相频特性')
plt.grid(True, which='both')

# 振荡环节
plt.subplot(5, 2, 9)
mag, phase, _ = control.bode(sys_second_order, omega, dB=True, deg=True, plot=False)
plt.semilogx(omega, 20 * np.log10(mag))
plt.title('振荡环节 幅频特性')
plt.grid(True, which='both')

plt.subplot(5, 2, 10)
plt.semilogx(omega, phase)
plt.title('振荡环节 相频特性')
plt.grid(True, which='both')

plt.tight_layout()
plt.show()

这个示例展示了如何使用 Python 绘制比例、积分、微分、惯性和振荡等典型环节的 Bode 图。通过观察这些典型环节的频率特性，我们可以更好地理解复杂系统的频率响应。

5-3 频率稳定判据

频率稳定判据是利用系统的开环频率特性来判断闭环系统稳定性的方法。最常用的频率稳定判据是奈奎斯特稳定判据（Nyquist Stability Criterion）。

奈奎斯特稳定判据

奈奎斯特稳定判据的基本思想是：通过研究系统开环频率特性 (G(jw)H(jw)) 的奈奎斯特图（Nyquist Plot），判断闭环系统的稳定性。具体判据如下：

Python 实现示例

下面我们使用 Python 绘制系统的奈奎斯特图，并判断系统的稳定性：

import control
import matplotlib.pyplot as plt
import numpy as np

# 设置中文显示
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams["axes.unicode_minus"] = False  # 解决负号显示问题

# 定义系统开环传递函数
# 示例1: 稳定系统 G(s) = 1 / (s^2 + s + 1)
num1 = [1]
den1 = [1, 1, 1]
sys1 = control.TransferFunction(num1, den1)

# 示例2: 不稳定系统 G(s) = 1 / (s^2 - s - 1)
num2 = [1]
den2 = [1, -1, -1]
sys2 = control.TransferFunction(num2, den2)

# 示例3: 临界稳定系统 G(s) = 1 / (s^2 + 1)
num3 = [1]
den3 = [1, 0, 1]
sys3 = control.TransferFunction(num3, den3)

# 绘制奈奎斯特图
plt.figure(figsize=(15, 5))

# 定义频率范围
omega_low = 0.01  # 避免使用零频率
omega_high = 10

# 稳定系统
plt.subplot(1, 3, 1)
control.nyquist_plot(sys1, omega_limits=(omega_low, omega_high))
plt.title('稳定系统 奈奎斯特图')
plt.grid(True)
plt.plot(-1, 0, 'ro')  # 标记(-1,0)点

# 不稳定系统
plt.subplot(1, 3, 2)
control.nyquist_plot(sys2, omega_limits=(omega_low, omega_high))
plt.title('不稳定系统 奈奎斯特图')
plt.grid(True)
plt.plot(-1, 0, 'ro')  # 标记(-1,0)点

# 临界稳定系统
plt.subplot(1, 3, 3)
control.nyquist_plot(sys3, omega_limits=(omega_low, omega_high))
plt.title('临界稳定系统 奈奎斯特图')
plt.grid(True)
plt.plot(-1, 0, 'ro')  # 标记(-1,0)点

plt.tight_layout()
plt.show()

# 判断系统稳定性
def check_stability(sys):
    poles = control.pole(sys)
    P = sum(np.real(poles) > 0)  # 右半平面的极点数
    # 计算奈奎斯特曲线包围(-1,0)点的圈数
    # 这里简化处理，实际应用中需要更复杂的算法来计算圈数
    # 对于示例系统，我们可以直接分析
    # 示例1: P=0，不包围(-1,0)点，稳定
    # 示例2: P=1，不包围(-1,0)点，不稳定
    # 示例3: P=0，通过(-1,0)点，临界稳定
    return P

print("系统稳定性分析:")
print(f"示例1: P = {check_stability(sys1)}，闭环系统稳定")
print(f"示例2: P = {check_stability(sys2)}，闭环系统不稳定")
print(f"示例3: P = {check_stability(sys3)}，闭环系统临界稳定")

这个示例展示了如何使用 Python 绘制系统的奈奎斯特图，并判断系统的稳定性。我们定义了三个系统：稳定系统、不稳定系统和临界稳定系统，分别绘制了它们的奈奎斯特图，并分析了稳定性。

5-4 稳定裕度

稳定裕度是衡量系统稳定性的一个重要指标，它表示系统距离不稳定状态的 “距离”。常用的稳定裕度有相位裕度（Phase Margin, PM）和增益裕度（Gain Margin, GM）。

Python 实现示例

下面我们使用 Python 计算并绘制系统的稳定裕度：

import control
import matplotlib.pyplot as plt
import numpy as np

# 设置中文显示
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams["axes.unicode_minus"] = False  # 解决负号显示问题

# 定义系统开环传递函数
# G(s) = 10 / (s(s+1)(s+2))
num = [10]
den = [1, 3, 2, 0]
sys = control.TransferFunction(num, den)

# 计算稳定裕度
gm, pm, wg, wc = control.margin(sys)

print(f"增益裕度: {gm:.4f} (dB: {20 * np.log10(gm):.4f} dB)")
print(f"相位裕度: {pm:.4f} 度")
print(f"相位穿越频率: {wg:.4f} rad/s")
print(f"截止频率: {wc:.4f} rad/s")

# 绘制Bode图并标注稳定裕度
plt.figure(figsize=(10, 8))
mag, phase, omega = control.bode(sys, dB=True, deg=True, plot=True)

# 标注增益裕度
plt.subplot(2, 1, 1)
plt.axvline(x=wg, color='r', linestyle='--')
plt.axhline(y=20 * np.log10(gm), color='g', linestyle='--')
plt.text(wg, 20 * np.log10(gm), f'GM={20 * np.log10(gm):.2f} dB', ha='right', va='bottom')

# 标注相位裕度
plt.subplot(2, 1, 2)
plt.axvline(x=wc, color='r', linestyle='--')
plt.axhline(y=-180 + pm, color='g', linestyle='--')
plt.text(wc, -180 + pm, f'PM={pm:.2f}°', ha='right', va='bottom')

plt.tight_layout()
plt.show()

# 分析不同增益下的稳定裕度
gains = [1, 5, 10, 20]
plt.figure(figsize=(12, 10))

for i, K in enumerate(gains):
    sys_K = K * sys
    gm, pm, wg, wc = control.margin(sys_K)

    plt.subplot(2, 2, i + 1)
    control.bode(sys_K, dB=True, deg=True, plot=True)
    plt.title(f'K={K}, PM={pm:.2f}°, GM={20 * np.log10(gm):.2f} dB')

    # 标注截止频率和相位穿越频率
    plt.subplot(2, 2, i + 1)
    plt.subplot(2, 1, 1)
    plt.axvline(x=wc, color='r', linestyle='--')
    plt.subplot(2, 1, 2)
    plt.axvline(x=wg, color='r', linestyle='--')

plt.tight_layout()
plt.show()

这个示例展示了如何使用 Python 计算并绘制系统的稳定裕度。我们定义了一个三阶系统，计算了其增益裕度和相位裕度，并在 Bode 图上进行了标注。此外，我们还分析了不同增益对稳定裕度的影响。

5-5 闭环系统的频域性能指标

闭环系统的频域性能指标是通过系统的闭环频率特性来描述系统性能的指标。主要包括：

带宽频率（Bandwidth Frequency, ωb）：闭环幅频特性下降到零频率值的 0.707 倍（即 - 3dB）时的频率。带宽反映了系统对不同频率输入信号的跟踪能力，带宽越大，系统响应速度越快，但抗噪声能力越弱。
谐振峰值（Resonant Peak, Mr）：闭环幅频特性的最大值。谐振峰值反映了系统的相对稳定性，Mr 越大，系统超调量越大，稳定性越差。
谐振频率（Resonant Frequency, ωr）：闭环幅频特性达到最大值时的频率。谐振频率反映了系统的振荡频率。

Python 实现示例

下面我们使用 Python 计算并分析闭环系统的频域性能指标：

import control
import matplotlib.pyplot as plt
import numpy as np

# 设置中文显示
plt.rcParams["font.family"] = ["SimHei"]
plt.rcParams["axes.unicode_minus"] = False  # 解决负号显示问题

# 定义系统开环传递函数
# G(s) = 10 / (s(s+1)(s+5))
num = [10]
den = [1, 6, 5, 0]
sys_open = control.TransferFunction(num, den)

# 计算闭环传递函数
sys_closed = control.feedback(sys_open, 1)

# 计算闭环频率响应
omega = np.logspace(-2, 2, 1000)
mag, phase, omega = control.bode(sys_closed, omega, dB=True, deg=True, plot=False)

# 计算零频率值
mag_0 = 20 * np.log10(np.abs(control.evalfr(sys_closed, 0)))

# 计算带宽频率(下降到-3dB的频率)
bandwidth_idx = np.where(20 * np.log10(mag) <= mag_0 - 3)[0][0]
bandwidth = omega[bandwidth_idx]

# 计算谐振峰值和谐振频率
mag_db = 20 * np.log10(mag)
Mr_db = np.max(mag_db)
Mr = 10 ** (Mr_db / 20)
resonant_idx = np.argmax(mag_db)
resonant_freq = omega[resonant_idx]

print(f"零频率幅值: {mag_0:.4f} dB")
print(f"带宽频率: {bandwidth:.4f} rad/s")
print(f"谐振峰值: {Mr_db:.4f} dB ({Mr:.4f})")
print(f"谐振频率: {resonant_freq:.4f} rad/s")

# 绘制闭环频率响应
plt.figure(figsize=(10, 8))

# 幅频特性
plt.subplot(2, 1, 1)
plt.semilogx(omega, 20 * np.log10(mag))
plt.axhline(y=mag_0, color='r', linestyle='--', label=f'零频率幅值: {mag_0:.2f} dB')
plt.axhline(y=mag_0 - 3, color='g', linestyle='--', label='-3dB线')
plt.axvline(x=bandwidth, color='g', linestyle='--', label=f'带宽频率: {bandwidth:.2f} rad/s')
plt.axhline(y=Mr_db, color='m', linestyle='--', label=f'谐振峰值: {Mr_db:.2f} dB')
plt.axvline(x=resonant_freq, color='m', linestyle='--', label=f'谐振频率: {resonant_freq:.2f} rad/s')
plt.title('闭环系统幅频特性')
plt.xlabel('频率 (rad/s)')
plt.ylabel('幅值 (dB)')
plt.grid(True, which='both')
plt.legend()

# 相频特性
plt.subplot(2, 1, 2)
plt.semilogx(omega, phase)
plt.title('闭环系统相频特性')
plt.xlabel('频率 (rad/s)')
plt.ylabel('相位 (度)')
plt.grid(True, which='both')

plt.tight_layout()
plt.show()

# 分析不同增益对闭环频域性能的影响
gains = [1, 5, 10, 20]
plt.figure(figsize=(12, 10))

for i, K in enumerate(gains):
    sys_open_K = K * sys_open
    sys_closed_K = control.feedback(sys_open_K, 1)

    # 计算闭环频率响应
    mag, phase, omega = control.bode(sys_closed_K, omega, dB=True, deg=True, plot=False)

    # 计算带宽频率
    mag_0 = 20 * np.log10(np.abs(control.evalfr(sys_closed_K, 0)))
    bandwidth_idx = np.where(20 * np.log10(mag) <= mag_0 - 3)[0][0]
    bandwidth = omega[bandwidth_idx]

    # 计算谐振峰值和谐振频率
    mag_db = 20 * np.log10(mag)
    Mr_db = np.max(mag_db)
    Mr = 10 ** (Mr_db / 20)
    resonant_idx = np.argmax(mag_db)
    resonant_freq = omega[resonant_idx]

    plt.subplot(2, 2, i + 1)
    plt.semilogx(omega, 20 * np.log10(mag))
    plt.axhline(y=mag_0 - 3, color='g', linestyle='--')
    plt.axvline(x=bandwidth, color='g', linestyle='--')
    plt.title(f'K={K}, 带宽={bandwidth:.2f} rad/s, Mr={Mr_db:.2f} dB')
    plt.xlabel('频率 (rad/s)')
    plt.ylabel('幅值 (dB)')
    plt.grid(True, which='both')

plt.tight_layout()
plt.show()

这个示例展示了如何使用 Python 计算并分析闭环系统的频域性能指标。我们定义了一个三阶系统，计算了其带宽频率、谐振峰值和谐振频率，并在幅频特性图上进行了标注。此外，我们还分析了不同增益对闭环频域性能的影响。

5-6 控制系统频域设计

控制系统的频域设计是指利用系统的频率特性来设计控制器，使系统满足给定的性能指标。常见的频域设计方法包括：

增益调整：通过调整系统的开环增益来满足稳态性能和稳定性要求。
超前校正：通过引入超前校正网络来提高系统的相位裕度，改善系统的动态性能。
滞后校正：通过引入滞后校正网络来提高系统的稳态精度，同时保持系统的稳定性。
滞后 - 超前校正：综合利用滞后校正和超前校正的优点，同时改善系统的稳态性能和动态性能。

Python 实现示例

下面我们使用 Python 实现一个控制系统的频域设计：

import control
import matplotlib.pyplot as plt
import numpy as np

# 设置中文显示
plt.rcParams["font.family"] = ["SimHei"]
plt.rcParams["axes.unicode_minus"] = False  # 解决负号显示问题

# 原系统开环传递函数
# G(s) = 1 / (s(s+1)(s+5))
num_original = [1]
den_original = [1, 6, 5, 0]
sys_original = control.TransferFunction(num_original, den_original)

# 设计要求: 相位裕度 ≥ 45°, 增益裕度 ≥ 10dB
# 计算原系统的稳定裕度
gm_original, pm_original, wg_original, wc_original = control.margin(sys_original)
print(f"原系统相位裕度: {pm_original:.2f}°, 增益裕度: {20*np.log10(gm_original):.2f} dB")

# 设计超前校正器
# 计算需要增加的相位裕度
pm_required = 45
pm_additional = pm_required - pm_original + 5  # 增加5°的余量

# 计算超前校正器参数
alpha = (1 + np.sin(np.radians(pm_additional))) / (1 - np.sin(np.radians(pm_additional)))
print(f"超前校正器参数 α = {alpha:.4f}")

# 选择校正器的转折频率
# 通常选择在原系统相位裕度为-90°左右的频率
omega_m = wc_original  # 选择在原系统截止频率附近
T = 1 / (omega_m * np.sqrt(alpha))
print(f"超前校正器转折频率: {omega_m:.4f} rad/s")
print(f"超前校正器时间常数: {T:.4f} s")

# 超前校正器传递函数
num_c = [T, 1]
den_c = [T * alpha, 1]
sys_c = control.TransferFunction(num_c, den_c)

# 校正后系统开环传递函数
sys_compensated = sys_c * sys_original

# 计算校正后系统的稳定裕度
gm_compensated, pm_compensated, wg_compensated, wc_compensated = control.margin(sys_compensated)
print(f"校正后系统相位裕度: {pm_compensated:.2f}°, 增益裕度: {20*np.log10(gm_compensated):.2f} dB")

# 绘制Bode图比较
plt.figure(figsize=(12, 10))

# 原系统Bode图
plt.subplot(2, 2, 1)
control.bode(sys_original, dB=True, deg=True, plot=True)
plt.title('原系统Bode图')

# 校正后系统Bode图
plt.subplot(2, 2, 2)
control.bode(sys_compensated, dB=True, deg=True, plot=True)
plt.title('校正后系统Bode图')

# 比较原系统和校正后系统的阶跃响应
plt.subplot(2, 2, 3)
# 获取原系统的阶跃响应
t_original, y_original = control.step_response(control.feedback(sys_original, 1))
# 获取校正后系统的阶跃响应
t_compensated, y_compensated = control.step_response(control.feedback(sys_compensated, 1))

# 确保时间向量和响应向量长度一致
min_len_original = min(len(t_original), len(y_original))
min_len_compensated = min(len(t_compensated), len(y_compensated))

plt.plot(t_original[:min_len_original], y_original[:min_len_original], 'b-', label='原系统')
plt.plot(t_compensated[:min_len_compensated], y_compensated[:min_len_compensated], 'r-', label='校正后系统')
plt.title('阶跃响应比较')
plt.xlabel('时间 (s)')
plt.ylabel('输出')
plt.legend()
plt.grid(True)

# 比较原系统和校正后系统的奈奎斯特图
plt.subplot(2, 2, 4)
omega_low = 0.01  # 避免使用零频率
control.nyquist_plot(sys_original, omega_limits=(omega_low, 10))
control.nyquist_plot(sys_compensated, omega_limits=(omega_low, 10))
plt.title('奈奎斯特图比较')
plt.plot(-1, 0, 'ro')  # 标记(-1,0)点
plt.grid(True)

plt.tight_layout()
plt.show()

# 分析系统性能
print("\n原系统性能分析:")
poles_original = control.pole(control.feedback(sys_original, 1))
print(f"  闭环极点: {poles_original}")
print(f"  相位裕度: {pm_original:.2f}°")
print(f"  增益裕度: {20*np.log10(gm_original):.2f} dB")

print("\n校正后系统性能分析:")
poles_compensated = control.pole(control.feedback(sys_compensated, 1))
print(f"  闭环极点: {poles_compensated}")
print(f"  相位裕度: {pm_compensated:.2f}°")
print(f"  增益裕度: {20*np.log10(gm_compensated):.2f} dB")

# 计算超调量
def calculate_overshoot(y):
    max_value = np.max(y)
    steady_state = y[-1]
    return (max_value - steady_state) / steady_state * 100

print(f"  原系统超调量: {calculate_overshoot(y_original):.2f}%")
print(f"  校正后系统超调量: {calculate_overshoot(y_compensated):.2f}%")

# 计算调节时间(2%误差带)
def calculate_settling_time(t, y):
    steady_state = y[-1]
    for i in range(len(y)):
        if np.all(np.abs(y[i:] - steady_state) <= 0.02 * steady_state):
            return t[i]
    return t[-1]

print(f"  原系统调节时间: {calculate_settling_time(t_original, y_original):.2f}s")
print(f"  校正后系统调节时间: {calculate_settling_time(t_compensated, y_compensated):.2f}s")