引言

在控制工程、信号处理和系统分析中，传递函数是描述线性时不变（LTI）系统输入输出关系的核心工具。它将复杂的微分方程转化为简洁的代数形式，极大地简化了系统分析和设计过程。本文将详细介绍系统建模求传递函数的实用方法，并深入解析常见问题，帮助读者掌握从物理模型到数学模型的完整流程。

一、系统建模的基本概念

1.1 什么是传递函数？

传递函数定义为系统输出拉普拉斯变换与输入拉普拉斯变换之比（初始条件为零）： $\( G(s) = \frac{Y(s)}{X(s)} \)$ 其中 ( Y(s) ) 是输出信号的拉普拉斯变换，( X(s) ) 是输入信号的拉普拉斯变换，( s ) 是复频率变量。

1.2 系统建模的步骤

  1. 物理建模：根据系统的工作原理，建立物理方程（如牛顿定律、基尔霍夫定律等）
  2. 数学建模：将物理方程转化为微分方程
  3. 拉普拉斯变换：对微分方程进行拉普拉斯变换，得到代数方程
  4. 传递函数求解：整理代数方程，得到传递函数表达式

二、实用建模方法详解

2.1 机械系统建模

示例：弹簧-质量-阻尼系统

考虑一个简单的机械系统：质量为 ( m ) 的物体通过弹簧（刚度 ( k )）和阻尼器（阻尼系数 ( c )）连接到固定点。

步骤1：物理建模 根据牛顿第二定律： $\( m\ddot{x} = F_{\text{spring}} + F_{\text{damper}} + F_{\text{external}} \)$ 其中：

  • 弹簧力：( F_{\text{spring}} = -kx )
  • 阻尼力：( F_{\text{damper}} = -c\dot{x} )
  • 外力：( F_{\text{external}} = f(t) )

步骤2：微分方程 $\( m\ddot{x} + c\dot{x} + kx = f(t) \)$

步骤3：拉普拉斯变换（零初始条件） $\( ms^2X(s) + csX(s) + kX(s) = F(s) \)\( \)\( (ms^2 + cs + k)X(s) = F(s) \)$

步骤4：传递函数 $\( G(s) = \frac{X(s)}{F(s)} = \frac{1}{ms^2 + cs + k} \)$

Python代码示例（系统分析）

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

# 定义系统参数
m = 1.0  # 质量 (kg)
c = 0.5  # 阻尼系数 (N·s/m)
k = 2.0  # 弹簧刚度 (N/m)

# 创建传递函数
num = [1]  # 分子多项式系数
den = [m, c, k]  # 分母多项式系数
sys = signal.TransferFunction(num, den)

# 系统分析
print("传递函数：", sys)
print("极点：", np.roots(den))

# 时域响应
t, y = signal.step(sys)
plt.figure(figsize=(10, 6))
plt.plot(t, y)
plt.title('单位阶跃响应')
plt.xlabel('时间 (s)')
plt.ylabel('位移 (m)')
plt.grid(True)
plt.show()

2.2 电路系统建模

示例：RLC电路

考虑一个串联RLC电路，输入电压 ( v_i(t) )，输出电压 ( v_o(t) ) 在电容两端。

步骤1：基尔霍夫电压定律 $\( v_i(t) = v_R(t) + v_L(t) + v_C(t) \)$ 其中：

  • ( v_R(t) = Ri(t) )
  • ( v_L(t) = L\frac{di(t)}{dt} )
  • ( v_C(t) = \frac{1}{C}\int i(t)dt )

步骤2：微分方程 对电容电压求导： $\( \frac{dv_C}{dt} = \frac{i(t)}{C} \)\( 代入KVL方程： \)\( v_i(t) = Ri(t) + L\frac{di(t)}{dt} + v_C(t) \)\( 对 \( i(t) \) 求导： \)\( \frac{di}{dt} = \frac{1}{L}(v_i - Ri - v_C) \)\( 代入电容方程： \)\( \frac{dv_C}{dt} = \frac{1}{LC}(v_i - Ri - v_C) \)\( 又因为 \( i = C\frac{dv_C}{dt} \)，代入得： \)\( \frac{dv_C}{dt} = \frac{1}{LC}(v_i - RC\frac{dv_C}{dt} - v_C) \)\( 整理得二阶微分方程： \)\( LC\frac{d^2v_C}{dt^2} + RC\frac{dv_C}{dt} + v_C = v_i \)$

步骤3：拉普拉斯变换 $\( LCS^2V_C(s) + RCSV_C(s) + V_C(s) = V_i(s) \)\( \)\( (LCS^2 + RCS + 1)V_C(s) = V_i(s) \)$

步骤4：传递函数 $\( G(s) = \frac{V_C(s)}{V_i(s)} = \frac{1}{LCS^2 + RCS + 1} \)$

MATLAB代码示例（系统仿真）

% 定义电路参数
R = 1000;  % 电阻 (Ω)
L = 0.1;   % 电感 (H)
C = 1e-6;  % 电容 (F)

% 创建传递函数
num = [1];
den = [L*C, R*C, 1];
sys = tf(num, den);

% 系统分析
disp('传递函数：');
sys
disp('极点：');
roots(den)

% 频率响应
figure;
bode(sys);
grid on;
title('RLC电路频率响应');

% 时域响应
figure;
step(sys);
grid on;
title('单位阶跃响应');

2.3 电气-机械系统建模（机电系统）

示例：直流电动机

考虑一个直流电动机，输入电压 ( v_a(t) )，输出转速 ( \omega(t) )。

步骤1：电气方程 $\( v_a(t) = R_a i_a(t) + L_a \frac{di_a(t)}{dt} + e_b(t) \)$ 其中反电动势 ( e_b(t) = K_b \omega(t) )

步骤2：机械方程 $\( J \frac{d\omega(t)}{dt} = T_m(t) - T_L(t) \)$ 其中电磁转矩 ( T_m(t) = K_t i_a(t) )

步骤3：微分方程组 $\( \begin{cases} v_a = R_a i_a + L_a \frac{di_a}{dt} + K_b \omega \\ J \frac{d\omega}{dt} = K_t i_a - T_L \end{cases} \)$

步骤4：拉普拉斯变换（忽略负载转矩 ( T_L )） $\( \begin{cases} V_a(s) = R_a I_a(s) + L_a s I_a(s) + K_b \Omega(s) \\ J s \Omega(s) = K_t I_a(s) \end{cases} \)$

步骤5：传递函数 从第二个方程：( I_a(s) = \frac{J s}{K_t} \Omega(s) ) 代入第一个方程： $\( V_a(s) = (R_a + L_a s) \frac{J s}{K_t} \Omega(s) + K_b \Omega(s) \)\( \)\( V_a(s) = \left[ \frac{J s (R_a + L_a s)}{K_t} + K_b \right] \Omega(s) \)\( \)\( G(s) = \frac{\Omega(s)}{V_a(s)} = \frac{K_t}{J L_a s^2 + J R_a s + K_t K_b} \)$

Python代码示例（机电系统建模）

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

# 直流电动机参数
R_a = 1.0  # 电枢电阻 (Ω)
L_a = 0.01  # 电枢电感 (H)
J = 0.01  # 转动惯量 (kg·m²)
K_t = 0.1  # 转矩常数 (N·m/A)
K_b = 0.1  # 反电动势常数 (V/(rad/s))

# 传递函数
num = [K_t]
den = [J*L_a, J*R_a, K_t*K_b]
sys = signal.TransferFunction(num, den)

# 系统分析
print("直流电动机传递函数：")
print(f"G(s) = {K_t:.3f} / ({J*L_a:.3f}s² + {J*R_a:.3f}s + {K_t*K_b:.3f})")

# 伯德图
w, mag, phase = signal.bode(sys)
plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
plt.semilogx(w, mag)
plt.title('幅频特性')
plt.xlabel('频率 (rad/s)')
plt.ylabel('幅值 (dB)')
plt.grid(True)

plt.subplot(1, 2, 2)
plt.semilogx(w, phase)
plt.title('相频特性')
plt.xlabel('频率 (rad/s)')
plt.ylabel('相位 (°)')
plt.grid(True)
plt.tight_layout()
plt.show()

# 阶跃响应
t, y = signal.step(sys)
plt.figure(figsize=(10, 6))
plt.plot(t, y)
plt.title('直流电动机单位阶跃响应')
plt.xlabel('时间 (s)')
plt.ylabel('转速 (rad/s)')
plt.grid(True)
plt.show()

三、常见问题解析

3.1 问题1：非线性系统如何处理？

问题描述：实际系统往往存在非线性（如饱和、死区、摩擦等），如何用传递函数近似？

解决方案：

  1. 线性化：在工作点附近进行泰勒展开
  2. 分段线性化：将非线性特性分段近似为线性
  3. 描述函数法：用于分析非线性系统的稳定性

示例：带饱和特性的放大器

import numpy as np
import matplotlib.pyplot as plt

# 非线性饱和特性
def saturate(x, limit=1.0):
    return np.clip(x, -limit, limit)

# 线性化近似（在工作点x0=0.5处）
x0 = 0.5
limit = 1.0
# 在x0处的导数（斜率）
if abs(x0) < limit:
    slope = 1.0  # 线性区斜率为1
else:
    slope = 0.0  # 饱和区斜率为0

# 绘制对比图
x = np.linspace(-2, 2, 1000)
y_nonlinear = saturate(x)
y_linear = slope * (x - x0) + saturate(x0)  # 线性化近似

plt.figure(figsize=(10, 6))
plt.plot(x, y_nonlinear, 'b-', linewidth=2, label='非线性饱和特性')
plt.plot(x, y_linear, 'r--', linewidth=2, label=f'线性化近似 (x₀={x0})')
plt.axvline(x0, color='g', linestyle=':', label=f'工作点 x₀={x0}')
plt.axhline(limit, color='k', linestyle=':', alpha=0.5)
plt.axhline(-limit, color='k', linestyle=':', alpha=0.5)
plt.title('非线性饱和特性的线性化近似')
plt.xlabel('输入')
plt.ylabel('输出')
plt.legend()
plt.grid(True)
plt.show()

print(f"在x₀={x0}处的线性化斜率: {slope}")

3.2 问题2：如何处理时变系统？

问题描述：系统参数随时间变化（如温度变化导致电阻变化），传递函数不再恒定。

解决方案：

  1. 准静态假设：假设参数变化缓慢，使用时变传递函数
  2. 自适应控制：实时调整控制器参数
  3. 鲁棒控制：设计对参数变化不敏感的控制器

示例：时变参数系统

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

# 时变参数：电阻随时间变化
def R_t(t):
    return 1000 + 200 * np.sin(2*np.pi*t/10)  # 周期性变化

# 时变传递函数
def time_varying_tf(t, R):
    L = 0.1
    C = 1e-6
    num = [1]
    den = [L*C, R*C, 1]
    return signal.TransferFunction(num, den)

# 仿真
t = np.linspace(0, 20, 1000)
R_values = R_t(t)
outputs = []

# 简化仿真：假设每个时间点系统瞬时响应
for i, t_i in enumerate(t[:100]):  # 简化计算，只取前100点
    sys = time_varying_tf(t_i, R_values[i])
    # 计算阶跃响应（简化）
    t_resp, y_resp = signal.step(sys, T=np.linspace(0, 0.1, 100))
    outputs.append(y_resp[-1])  # 取最终值

plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
plt.plot(t[:100], R_values[:100])
plt.title('时变电阻 R(t)')
plt.xlabel('时间 (s)')
plt.ylabel('电阻 (Ω)')
plt.grid(True)

plt.subplot(1, 2, 2)
plt.plot(t[:100], outputs)
plt.title('系统响应（时变参数）')
plt.xlabel('时间 (s)')
plt.ylabel('输出')
plt.grid(True)
plt.tight_layout()
plt.show()

3.3 问题3：如何处理高阶系统？

问题描述：实际系统可能非常复杂，包含多个状态变量，导致高阶微分方程。

解决方案：

  1. 降阶处理：忽略次要动态，保留主要动态
  2. 主导极点法：保留主导极点，忽略非主导极点
  3. 状态空间模型：使用状态空间表示法，便于数值计算

示例：高阶系统降阶

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

# 高阶系统：5阶系统
num = [1]
den = [1, 3, 6, 6, 3, 1]  # (s+1)^5 的展开

# 原始系统
sys_full = signal.TransferFunction(num, den)
print("原始5阶系统极点：", np.roots(den))

# 主导极点法降阶：保留最靠近虚轴的两个极点
poles = np.roots(den)
# 按实部绝对值排序（虚部相同）
sorted_poles = sorted(poles, key=lambda p: abs(p.real))
dominant_poles = sorted_poles[:2]  # 取两个主导极点

# 构建降阶系统（二阶）
# 假设主导极点为 -0.5 ± 0.866j
p1, p2 = dominant_poles
# 二阶系统特征方程：s² + 2ζω_n s + ω_n² = 0
# 从极点计算参数
omega_n = abs(p1)  # 自然频率
zeta = -p1.real / omega_n  # 阻尼比

# 降阶系统传递函数
num_reduced = [1]
den_reduced = [1, 2*zeta*omega_n, omega_n**2]
sys_reduced = signal.TransferFunction(num_reduced, den_reduced)

# 时域响应对比
t = np.linspace(0, 10, 1000)
t_full, y_full = signal.step(sys_full, T=t)
t_red, y_red = signal.step(sys_reduced, T=t)

plt.figure(figsize=(10, 6))
plt.plot(t_full, y_full, 'b-', linewidth=2, label='原始5阶系统')
plt.plot(t_red, y_red, 'r--', linewidth=2, label=f'降阶2阶系统 (ζ={zeta:.3f}, ω_n={omega_n:.3f})')
plt.title('高阶系统降阶对比')
plt.xlabel('时间 (s)')
plt.ylabel('输出')
plt.legend()
plt.grid(True)
plt.show()

print(f"主导极点：{dominant_poles}")
print(f"降阶系统参数：阻尼比 ζ={zeta:.3f}, 自然频率 ω_n={omega_n:.3f}")

3.4 问题4：如何验证传递函数的正确性？

问题描述：建立的传递函数是否准确反映了实际系统？

解决方案：

  1. 物理一致性检查：量纲分析、极限情况分析
  2. 实验验证：通过实际测量数据对比
  3. 仿真验证：与详细仿真结果对比

示例：传递函数验证

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

# 假设我们通过实验得到了一组数据
# 输入：单位阶跃信号
# 输出：测量数据（含噪声）
t_exp = np.linspace(0, 5, 50)
y_exp = 1 - np.exp(-0.5*t_exp) * (np.cos(0.866*t_exp) + 0.577*np.sin(0.866*t_exp))
# 添加噪声
np.random.seed(42)
noise = 0.02 * np.random.randn(len(t_exp))
y_exp += noise

# 假设我们建立的传递函数
num = [1]
den = [1, 1, 1]  # 二阶系统：s² + s + 1
sys = signal.TransferFunction(num, den)

# 仿真响应
t_sim, y_sim = signal.step(sys, T=t_exp)

# 计算误差
error = y_exp - y_sim
mse = np.mean(error**2)

# 绘制对比
plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
plt.plot(t_exp, y_exp, 'bo', label='实验数据', markersize=4)
plt.plot(t_sim, y_sim, 'r-', linewidth=2, label='传递函数仿真')
plt.title('实验数据与仿真对比')
plt.xlabel('时间 (s)')
plt.ylabel('输出')
plt.legend()
plt.grid(True)

plt.subplot(1, 2, 2)
plt.plot(t_exp, error, 'g-', linewidth=2)
plt.title(f'误差 (MSE={mse:.6f})')
plt.xlabel('时间 (s)')
plt.ylabel('误差')
plt.grid(True)
plt.tight_layout()
plt.show()

print(f"均方误差：{mse:.6f}")
print(f"传递函数：G(s) = 1/(s² + s + 1)")

四、高级建模技巧

4.1 状态空间模型与传递函数的转换

状态空间模型： $\( \begin{cases} \dot{x} = Ax + Bu \\ y = Cx + Du \end{cases} \)$

转换为传递函数： $\( G(s) = C(sI - A)^{-1}B + D \)$

Python实现：

import numpy as np
import control as ct

# 定义状态空间模型
A = np.array([[0, 1],
              [-1, -1]])
B = np.array([[0],
              [1]])
C = np.array([[1, 0]])
D = np.array([[0]])

# 创建状态空间系统
sys_ss = ct.ss(A, B, C, D)

# 转换为传递函数
sys_tf = ct.tf(sys_ss)

print("状态空间模型：")
print(sys_ss)
print("\n传递函数：")
print(sys_tf)

# 验证转换
print("\n极点验证：")
print("状态空间极点：", np.linalg.eigvals(A))
print("传递函数极点：", np.roots([1, 1, 1]))

4.2 多输入多输出（MIMO）系统建模

示例：2×2 MIMO系统

import numpy as np
import control as ct

# 定义MIMO系统
A = np.array([[0, 1],
              [-2, -3]])
B = np.array([[0, 1],
              [1, 0]])
C = np.array([[1, 0],
              [0, 1]])
D = np.array([[0, 0],
              [0, 0]])

# 创建MIMO系统
sys_mimo = ct.ss(A, B, C, D)

# 获取传递函数矩阵
G = ct.tf(sys_mimo)
print("MIMO系统传递函数矩阵：")
print(G)

# 绘制频率响应
ct.bode(G)

4.3 离散时间系统建模

连续系统离散化： $\( G(z) = \mathcal{Z}\{G(s)\} \)$

Python实现：

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

# 连续系统
num = [1]
den = [1, 1, 1]  # s² + s + 1
sys_cont = signal.TransferFunction(num, den)

# 离散化（采样时间T=0.1s）
T = 0.1
sys_disc = signal.cont2discrete(sys_cont, T, method='zoh')

print("连续系统：", sys_cont)
print("离散系统：", sys_disc)

# 时域响应对比
t_cont = np.linspace(0, 5, 1000)
t_disc = np.arange(0, 5, T)

_, y_cont = signal.step(sys_cont, T=t_cont)
_, y_disc = signal.step(sys_disc, T=t_disc)

plt.figure(figsize=(10, 6))
plt.plot(t_cont, y_cont, 'b-', linewidth=2, label='连续系统')
plt.plot(t_disc, y_disc, 'ro', markersize=6, label='离散系统 (T=0.1s)')
plt.title('连续系统与离散系统响应对比')
plt.xlabel('时间 (s)')
plt.ylabel('输出')
plt.legend()
plt.grid(True)
plt.show()

五、实用建议与最佳实践

5.1 建模流程检查清单

  1. 明确系统边界：确定输入输出变量
  2. 简化假设：明确哪些因素可以忽略
  3. 量纲检查：确保方程两边量纲一致
  4. 极限情况分析：验证模型在极端条件下的行为
  5. 参数敏感性分析：了解关键参数对系统的影响

5.2 常见错误及避免方法

  1. 错误1：忽略初始条件

    • 问题：传递函数假设零初始条件
    • 解决：对于非零初始条件，使用拉普拉斯变换的完整形式

  2. 错误2：非线性直接线性化

    • 问题：在整个工作范围内线性化
    • 解决：在工作点附近线性化，或使用分段线性化

  3. 错误3：单位不一致

    • 问题：混合使用不同单位制
    • 解决：统一使用国际单位制（SI）

  4. 错误4：过度简化

    • 问题：忽略重要动态导致模型不准确
    • 解决：通过实验验证模型精度

5.3 软件工具推荐

  1. MATLAB/Simulink：强大的系统建模和仿真工具
  2. Python（SciPy, Control库）：开源、灵活的数值计算
  3. Modelica：多领域物理建模语言
  4. LabVIEW：图形化系统设计平台

六、总结

系统建模求传递函数是一个系统化的过程，需要结合物理原理、数学工具和工程经验。通过本文介绍的方法，读者可以：

  1. 掌握机械、电气、机电系统的建模方法
  2. 理解并解决非线性、时变、高阶系统等常见问题
  3. 使用现代工具进行仿真和验证
  4. 避免常见建模错误

记住，好的模型不是最复杂的，而是在满足精度要求下最简单的。在实际工程中，模型精度与计算复杂度需要权衡，根据具体应用场景选择合适的建模方法。

最后建议：在实际项目中，建议先建立简单模型，然后通过实验数据逐步修正和细化，最终得到既准确又实用的系统模型。