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

啊阿狸不会拉杆

发布于 2026-01-21 11:40:31

1730

3-1 系统时间响应的性能指标

时域性能指标定义

线性系统的时域性能指标是衡量系统动态响应特性的重要参数，主要包括：

上升时间 (Rise Time)：响应从终值的 10% 上升到 90%（或 0 到 100%）所需的时间
峰值时间 (Peak Time)：响应达到第一个峰值所需的时间
超调量 (Overshoot)：响应峰值与终值的差值百分比
调整时间 (Settling Time)：响应进入并保持在终值 ±2% 或 ±5% 误差带内的时间
稳态误差 (Steady-State Error)：响应稳定后与期望值的偏差

性能指标计算与可视化

以下代码演示了如何使用 Python 计算二阶系统的时域性能指标：

import numpy as np
import matplotlib.pyplot as plt
from control import tf, step_response
from control.matlab import stepinfo  # 用于提取性能指标

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 定义二阶系统传递函数
def second_order_system(wn, zeta):
    num = [wn**2]
    den = [1, 2*zeta*wn, wn**2]
    return tf(num, den)

# 系统参数
wn = 10  # 自然频率
zeta = 0.5  # 阻尼比
sys = second_order_system(wn, zeta)

# 计算阶跃响应
t, y = step_response(sys, T=np.linspace(0, 2, 1000))

# 提取性能指标
info = stepinfo(sys)

# 绘制阶跃响应曲线
plt.figure(figsize=(10, 6))
plt.plot(t, y, 'b-', linewidth=2)
plt.axhline(y=1, color='r', linestyle='--', label='期望值')
plt.axhline(y=1+info['Overshoot']/100, color='g', linestyle=':', alpha=0.5)
plt.axhline(y=1-info['Overshoot']/100, color='g', linestyle=':', alpha=0.5)
plt.axhline(y=1.02, color='m', linestyle=':', alpha=0.5, label='±2%误差带')
plt.axhline(y=0.98, color='m', linestyle=':', alpha=0.5)
plt.axvline(x=info['SettlingTime'], color='k', linestyle='--', label=f'调整时间: {info["SettlingTime"]:.2f}s')
plt.scatter(info['PeakTime'], 1+info['Overshoot']/100, color='red', marker='o', label=f'峰值: {info["PeakTime"]:.2f}s')

plt.xlabel('时间(秒)')
plt.ylabel('系统输出')
plt.title('二阶系统阶跃响应及时域性能指标')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

# 打印性能指标
print("时域性能指标:")
print(f"上升时间: {info['RiseTime']:.4f} 秒")
print(f"峰值时间: {info['PeakTime']:.4f} 秒")
print(f"超调量: {info['Overshoot']:.2f} %")
print(f"调整时间: {info['SettlingTime']:.4f} 秒")
print(f"稳态值: {info['SteadyStateValue']:.4f}")
print(f"稳态误差: {1 - info['SteadyStateValue']:.4f}")

性能指标思维导图

图 3-1 时域性能指标思维导图

3-2 一阶系统的时域分析

一阶系统数学模型

一阶系统的标准传递函数形式为：

其中 T 为时间常数，决定了系统的响应速度。

一阶系统阶跃响应特性

单位阶跃响应：

时间常数 T 越大，响应越慢
当 t = T 时，(c(t) = 0.632)，即响应达到终值的 63.2%
调整时间

不同时间常数的响应对比

import numpy as np
import matplotlib.pyplot as plt
from control import tf, step_response

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 定义一阶系统传递函数
def first_order_system(T):
    num = [1]
    den = [T, 1]
    return tf(num, den)

# 不同时间常数
T_values = [0.1, 0.5, 1.0, 2.0]

# 绘制不同时间常数的阶跃响应
plt.figure(figsize=(10, 6))

for T in T_values:
    sys = first_order_system(T)
    t, y = step_response(sys, T=np.linspace(0, 5, 1000))
    plt.plot(t, y, label=f'T = {T}s')
    # 标记t=T时的响应值
    idx = np.argmin(np.abs(t - T))
    plt.scatter(T, y[idx], color='red', marker='o', s=100)
    plt.text(T+0.1, y[idx]+0.05, f'63.2%', fontsize=9)

plt.axhline(y=1, color='r', linestyle='--', label='期望值')
plt.xlabel('时间(秒)')
plt.ylabel('系统输出')
plt.title('一阶系统不同时间常数的阶跃响应')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

# 计算调整时间
print("一阶系统调整时间(±2%误差带):")
for T in T_values:
    sys = first_order_system(T)
    t, y = step_response(sys, T=np.linspace(0, 5*T, 1000))
    # 找到进入±2%误差带的时间
    idx = np.where(np.abs(y-1) < 0.02)[0]
    if len(idx) > 0:
        ts = t[idx[0]]
        print(f"T = {T}s, 调整时间: {ts:.2f}s (理论值: {4*T:.2f}s)")

一阶系统流程图

图 3-2 一阶系统分析流程图

3-3 二阶系统的时域分析

二阶系统标准模型

二阶系统的标准传递函数形式为：

其中：

w为自然频率，决定响应速度
ζ为阻尼比，决定响应特性

阻尼比分类

欠阻尼 (0 <ζ < 1)：振荡响应，超调量存在
临界阻尼 (ζ = 1)：非振荡响应，最快调整时间
过阻尼 (ζ> 1)：非振荡响应，调整时间较长
无阻尼 (ζ = 0)：等幅振荡

不同阻尼比的响应对比

import numpy as np
import matplotlib.pyplot as plt
from control import tf, step_response

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 定义二阶系统传递函数
def second_order_system(wn, zeta):
    num = [wn**2]
    den = [1, 2*zeta*wn, wn**2]
    return tf(num, den)

# 固定自然频率，改变阻尼比
wn = 5  # 自然频率
zeta_values = [0.2, 0.5, 1.0, 1.5, 0.0]  # 阻尼比(包括无阻尼)

# 绘制不同阻尼比的阶跃响应
plt.figure(figsize=(12, 8))

for zeta in zeta_values:
    sys = second_order_system(wn, zeta)
    t, y = step_response(sys, T=np.linspace(0, 5, 1000))
    
    # 标记阻尼比类型
    if zeta == 0:
        zeta_text = "无阻尼(ζ=0)"
    elif zeta < 1:
        zeta_text = f"欠阻尼(ζ={zeta})"
    elif zeta == 1:
        zeta_text = "临界阻尼(ζ=1)"
    else:
        zeta_text = f"过阻尼(ζ={zeta})"
    
    plt.plot(t, y, label=zeta_text)

plt.axhline(y=1, color='r', linestyle='--', label='期望值')
plt.xlabel('时间(秒)')
plt.ylabel('系统输出')
plt.title('二阶系统不同阻尼比的阶跃响应')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

# 计算欠阻尼系统性能指标
print("欠阻尼二阶系统性能指标:")
zeta = 0.5
sys = second_order_system(wn, zeta)
from control.matlab import stepinfo
info = stepinfo(sys)

print(f"阻尼比 ζ = {zeta}")
print(f"自然频率 ωn = {wn} rad/s")
print(f"上升时间: {info['RiseTime']:.4f} 秒")
print(f"峰值时间: {info['PeakTime']:.4f} 秒")
print(f"超调量: {info['Overshoot']:.2f} %")
print(f"调整时间: {info['SettlingTime']:.4f} 秒")

二阶系统性能指标公式

图 3-3 二阶系统性能指标公式

3-4 高阶系统的时域分析

高阶系统模型

高阶系统的传递函数一般形式为：

主导极点概念

距离虚轴最近的极点，且附近无零点
其他极点距离虚轴的距离是主导极点的 5 倍以上
高阶系统可近似为以主导极点为特征的低阶系统

高阶系统与主导极点近似对比

import numpy as np
import matplotlib.pyplot as plt
from control import tf, step_response

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 定义高阶系统传递函数(三阶系统)
def high_order_system():
    # 三阶系统: 主导极点(-1+1j,-1-1j)和远极点(-10)
    num = [1]
    den = [1, 12, 38, 40]  # (s+10)(s^2+2s+2)
    return tf(num, den)

# 定义主导极点对应的二阶系统
def dominant_pole_system():
    num = [2]  # 调整增益使稳态值相同
    den = [1, 2, 2]
    return tf(num, den)

# 获取系统
sys_high = high_order_system()
sys_dominant = dominant_pole_system()

# 计算阶跃响应
t_high, y_high = step_response(sys_high, T=np.linspace(0, 10, 1000))
t_dominant, y_dominant = step_response(sys_dominant, T=np.linspace(0, 10, 1000))

# 绘制响应曲线
plt.figure(figsize=(10, 6))
plt.plot(t_high, y_high, 'b-', label='三阶系统')
plt.plot(t_dominant, y_dominant, 'r--', label='主导极点二阶系统')
plt.axhline(y=1, color='g', linestyle='--', label='期望值')
plt.xlabel('时间(秒)')
plt.ylabel('系统输出')
plt.title('高阶系统与主导极点近似系统的阶跃响应对比')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

# 计算性能指标
from control.matlab import stepinfo
print("三阶系统性能指标:")
info_high = stepinfo(sys_high)
print(f"上升时间: {info_high['RiseTime']:.4f} 秒")
print(f"峰值时间: {info_high['PeakTime']:.4f} 秒")
print(f"超调量: {info_high['Overshoot']:.2f} %")
print(f"调整时间: {info_high['SettlingTime']:.4f} 秒")

print("\n主导极点二阶系统性能指标:")
info_dominant = stepinfo(sys_dominant)
print(f"上升时间: {info_dominant['RiseTime']:.4f} 秒")
print(f"峰值时间: {info_dominant['PeakTime']:.4f} 秒")
print(f"超调量: {info_dominant['Overshoot']:.2f} %")
print(f"调整时间: {info_dominant['SettlingTime']:.4f} 秒")

# 计算误差
error_rise = abs(info_high['RiseTime'] - info_dominant['RiseTime']) / info_high['RiseTime'] * 100
error_overshoot = abs(info_high['Overshoot'] - info_dominant['Overshoot']) / info_high['Overshoot'] * 100
print(f"\n上升时间近似误差: {error_rise:.2f}%")
print(f"超调量近似误差: {error_overshoot:.2f}%")

高阶系统分析流程

图 3-4 高阶系统分析流程图

3-5 线性系统的稳定性分析

稳定性基本概念

稳定系统：对有界输入产生有界输出 (BIBO 稳定)
渐近稳定：自由响应随时间衰减到零
临界稳定：自由响应为等幅振荡
不稳定：自由响应随时间增长

劳斯判据实现

import numpy as np

def routh_criterion(coefficients):
    """
    实现劳斯稳定性判据
    coefficients: 特征多项式系数[an, an-1, ..., a0]
    返回: (是否稳定, 正实部根的个数)
    """
    n = len(coefficients) - 1
    if n < 1:
        return True, 0

    # 构建劳斯表
    routh_table = np.zeros((n + 1, (n + 1) // 2 + 1))

    # 第一行和第二行
    routh_table[0, :len(coefficients[0::2])] = coefficients[0::2]
    routh_table[1, :len(coefficients[1::2])] = coefficients[1::2]

    # 填充劳斯表
    for i in range(2, n + 1):
        for j in range(routh_table.shape[1] - 1):
            divisor = routh_table[i - 1, 0]
            if divisor == 0:
                print("劳斯表中出现零，使用特殊处理")
                return False, "需要使用赫尔维茨判据或其他方法"
            routh_table[i, j] = -np.linalg.det([[routh_table[i - 2, j], routh_table[i - 2, j + 1]],
                                                 [routh_table[i - 1, j], routh_table[i - 1, j + 1]]]) / divisor

    # 检查第一列符号变化
    sign_changes = 0
    for i in range(1, n + 1):
        if routh_table[i, 0] * routh_table[i - 1, 0] < 0:
            sign_changes += 1

    is_stable = sign_changes == 0
    return is_stable, sign_changes

# 测试案例
case1 = [1, 3, 3, 1]  # s^3+3s^2+3s+1=0, 根:-1,-1,-1(稳定)
case2 = [1, 2, 1, 0]  # s^3+2s^2+s=0, 根:0,-1,-1(临界稳定)
case3 = [1, -5, 6, 0]  # s^3-5s^2+6s=0, 根:0,2,3(不稳定)
case4 = [1, 2, 2, 4]  # s^3+2s^2+2s+4=0, 根:-2,±j√2(稳定)

for case in [case1, case2, case3, case4]:
    is_stable, num_unstable = routh_criterion(case)
    print(f"特征多项式: {' + '.join([f'{c}s^{i}' for i, c in enumerate(reversed(case)) if c != 0])}")
    if is_stable:
        print("系统状态: 稳定")
    else:
        print(f"系统状态: 不稳定, 正实部根的个数: {num_unstable}")
    print("-" * 50)

稳定性分析思维导图

图 3-5 稳定性分析方法思维导图

3-6 线性系统的稳态误差计算

稳态误差基本公式

稳态误差e_ss}可通过终值定理计算：

不同输入下的稳态误差

import numpy as np
from control import tf, evalfr


def 计算稳态误差(G, R_type, params=None):
    """
    计算线性系统的稳态误差
    G: 开环传递函数
    R_type: 输入类型('step', 'ramp', 'parabola')
    params: 输入参数，如幅值
    """
    if params is None:
        params = {'amplitude': 1}

    A = params['amplitude']
    s = 1e-5  # 使用一个小的s值代替0，避免除零错误

    if R_type == 'step':
        # 阶跃输入 R(s) = A/s
        Kp = evalfr(G, s)
        if np.abs(1 + Kp) < 1e-10:
            e_ss = np.inf
        else:
            e_ss = A / (1 + Kp)
    elif R_type == 'ramp':
        # 斜坡输入 R(s) = A/s^2
        # 使用s*G(s)在s处的值
        sG_value = s * evalfr(G, s)
        if np.abs(sG_value) < 1e-10:
            e_ss = np.inf
        else:
            e_ss = A / sG_value
    elif R_type == 'parabola':
        # 抛物线输入 R(s) = A/s^3
        # 使用s²*G(s)在s处的值
        s2G_value = s ** 2 * evalfr(G, s)
        if np.abs(s2G_value) < 1e-10:
            e_ss = np.inf
        else:
            e_ss = A / s2G_value
    else:
        raise ValueError("不支持的输入类型")

    # 如果结果是复数，返回实部（虚部通常很小）
    if isinstance(e_ss, complex):
        return np.real(e_ss)
    return e_ss


# 定义系统类型
# 0型系统: G(s) = K/(Ts+1)
G0 = tf([5], [1, 1])  # K=5, T=1
# 1型系统: G(s) = K/(s(Ts+1))
G1 = tf([5], [1, 1, 0])  # K=5, T=1
# 2型系统: G(s) = K/(s²(Ts+1))
G2 = tf([5], [1, 1, 0, 0])  # K=5, T=1

# 计算不同系统在不同输入下的稳态误差
输入类型 = ['step', 'ramp', 'parabola']
幅值 = 1.0

print("不同类型系统的稳态误差:")
print("-" * 60)
print(f"{'系统类型':<10}{'输入类型':<10}{'阶跃输入误差':<15}{'斜坡输入误差':<15}{'抛物线输入误差':<15}")
print("-" * 60)

for sys_type, G in zip(['0型系统', '1型系统', '2型系统'], [G0, G1, G2]):
    # 为每种系统类型计算所有输入类型的误差
    errors = []
    for r_type in 输入类型:
        error = 计算稳态误差(G, r_type, {'amplitude': 幅值})
        errors.append(f"{error:.4f}")

    # 每个系统类型打印一行
    print(f"{sys_type:<10}{'所有':<10}{errors[0]:<15}{errors[1]:<15}{errors[2]:<15}")

# 案例：单位反馈系统G(s)=10/(s(s+1))
print("\n" + "-" * 60)
G_case = tf([10], [1, 1, 0])
print("案例: 单位反馈系统 G(s)=10/(s(s+1))")
print(f"阶跃输入稳态误差: {计算稳态误差(G_case, 'step'):.4f}")
print(f"斜坡输入稳态误差: {计算稳态误差(G_case, 'ramp'):.4f}")
print(f"抛物线输入稳态误差: {计算稳态误差(G_case, 'parabola'):.4f}")

稳态误差规律表

系统类型	阶跃输入 (r (t)=R・1 (t))	斜坡输入 (r (t)=R・t)	抛物线输入 (r (t)=R・t²/2)
0 型系统	R/(1+Kp)	∞	∞
1 型系统	0	R/Kv	∞
2 型系统	0	0	R/Ka

表 3-1 不同类型系统的稳态误差

3-7 控制系统时域设计

时域设计基本概念

基于时域性能指标进行控制器设计
常用方法：极点配置法、PID 控制
设计目标：满足上升时间、超调量、调整时间等指标

PID 控制器设计案例

import numpy as np
import matplotlib.pyplot as plt
from control import tf, step_response, feedback

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 被控对象传递函数
G_plant = tf([1], [1, 2, 0])  # G(s) = 1/(s(s+2))

# 设计PID控制器
def 设计PID控制器(Kp, Ki, Kd):
    """设计PID控制器并返回闭环系统"""
    # PID控制器传递函数: Kp + Ki/s + Kd*s = (Kd*s^2 + Kp*s + Ki)/s
    num_pid = [Kd, Kp, Ki]
    den_pid = [1, 0]
    G_pid = tf(num_pid, den_pid)

    # 闭环传递函数
    G_closed = feedback(G_pid * G_plant, 1)
    return G_closed

# 不同PID参数测试
Kp_values = [5, 10, 15]
Ki_values = [1, 2, 3]
Kd_values = [0.5, 1.0, 1.5]

# 绘制不同Kp的响应
plt.figure(figsize=(12, 8))

for i, Kp in enumerate(Kp_values):
    plt.subplot(3, 1, i + 1)
    G_closed = 设计PID控制器(Kp, Ki_values[i], Kd_values[i])
    t, y = step_response(G_closed, T=np.linspace(0, 5, 1000))
    plt.plot(t, y, 'b-', label=f'Kp={Kp}, Ki={Ki_values[i]}, Kd={Kd_values[i]}')
    plt.axhline(y=1, color='r', linestyle='--', label='期望值')
    plt.xlabel('时间(秒)')
    plt.ylabel('系统输出')
    plt.title(f'PID参数对系统响应的影响 (Kp={Kp})')
    plt.legend()
    plt.grid(True)

plt.tight_layout()
plt.show()

# 优化PID参数设计
# 目标: 超调量<10%, 调整时间<1.5秒
Kp_opt = 12
Ki_opt = 2.5
Kd_opt = 1.2

G_closed_opt = 设计PID控制器(Kp_opt, Ki_opt, Kd_opt)
t_opt, y_opt = step_response(G_closed_opt, T=np.linspace(0, 5, 1000))

# 绘制优化后的响应
plt.figure(figsize=(10, 6))
plt.plot(t_opt, y_opt, 'g-', linewidth=2, label='优化PID控制')
plt.axhline(y=1, color='r', linestyle='--', label='期望值')
plt.axhline(y=1.1, color='m', linestyle=':', alpha=0.5, label='10%超调量')
plt.axhline(y=0.9, color='m', linestyle=':', alpha=0.5)
plt.axvline(x=1.5, color='k', linestyle='--', label='调整时间限制')
plt.xlabel('时间(秒)')
plt.ylabel('系统输出')
plt.title('优化PID控制器的系统阶跃响应')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

# 计算优化后性能指标
from control.matlab import stepinfo

info_opt = stepinfo(G_closed_opt)
print("优化PID控制器性能指标:")
print(f"上升时间: {info_opt['RiseTime']:.4f} 秒")
print(f"峰值时间: {info_opt['PeakTime']:.4f} 秒")
print(f"超调量: {info_opt['Overshoot']:.2f} %")
print(f"调整时间: {info_opt['SettlingTime']:.4f} 秒")
print(f"稳态误差: {1 - info_opt['SteadyStateValue']:.4f}")