电力系统暂态稳定计算与单机无穷大系统建模

电力系统暂态稳定计算与单机无穷大系统建模

1. 引言

1.1 电力系统暂态稳定概述

电力系统暂态稳定是指电力系统在遭受大扰动（如短路故障、线路切换、发电机退出运行等）后，各同步发电机保持同步运行并过渡到新的或恢复到原来稳态运行方式的能力。暂态稳定分析是电力系统安全运行的重要保障，对于电网规划、运行和控制具有重要意义。

在实际工程应用中，使用完整的省级电网参数进行暂态稳定计算会面临模型复杂、计算量大、耗时长等问题。特别是当需要多次仿真或进行参数灵敏度分析时，完整的详细模型会大大降低工作效率。因此，建立简化模型成为提高计算效率的有效途径。

1.2 单机无穷大系统模型的价值

单机无穷大系统（Single Machine Infinite Bus, SMIB）是电力系统稳定性分析中最基本且重要的简化模型。它将整个系统简化为一台等效发电机通过输电线路与一个电压和频率恒定的无穷大母线相连。这种简化模型具有以下优点：

1. 计算效率高：大大减少了计算复杂度，提高了仿真速度
2. 物理意义明确：便于理解电力系统稳定性的基本机理
3. 参数调整方便：易于进行参数灵敏度分析
4. 教育价值高：适合用于教学和基本原理研究

本文将详细介绍如何使用Python实现单机无穷大系统的建模与暂态稳定计算，为实际工程应用提供参考。

2. 单机无穷大系统数学模型

2.1 同步发电机模型

同步发电机是电力系统的核心元件，其动态行为可以用一组微分方程描述。对于暂态稳定分析，通常采用三阶模型，考虑发电机的转子运动方程和暂态电动势变化。

转子运动方程：

dδ/dt = ω - ω0
dω/dt = (Pm - Pe - D(ω - ω0)) / M

其中：

• δ：发电机功角（弧度）
• ω：发电机转子角速度（弧度/秒）
• ω0：同步角速度（通常为2πf0，f0=50Hz或60Hz）
• Pm：机械功率（标幺值）
• Pe：电磁功率（标幺值）
• D：阻尼系数（标幺值）
• M：惯性时间常数（秒）

暂态电动势方程：

dE'/dt = (Ef - E' - (Xd - X'd)Id) / T'd0

其中：

• E’：q轴暂态电动势（标幺值）
• Ef：励磁电动势（标幺值）
• Xd：d轴同步电抗（标幺值）
• X’d：d轴暂态电抗（标幺值）
• Id：d轴电流（标幺值）
• T’d0：d轴开路暂态时间常数（秒）

2.2 输电系统模型

单机无穷大系统中，发电机通过输电线路与无穷大母线相连。输电系统可以用一个等效阻抗表示：

Zline = Rline + jXline

其中：

• Rline：线路电阻（标幺值）
• Xline：线路电抗（标幺值）

2.3 无穷大母线

无穷大母线假设其电压幅值和频率恒定不变，不受所连接发电机运行状态的影响。这是一个理想化的概念，在实际系统中相当于一个容量远大于所研究发电机的系统。

2.4 系统完整数学模型

将发电机模型与输电系统模型结合，可以得到单机无穷大系统的完整数学模型：

dδ/dt = ω - ω0
dω/dt = (Pm - Pe - D(ω - ω0)) / M
dE'/dt = (Ef - E' - (Xd - X'd)Id) / T'd0Pe = (E'V∞/XΣ)sinδ + (V∞²/2)(1/Xq - 1/Xd)sin2δ
Id = (E' - V∞cosδ)/X'dΣ

其中：

• V∞：无穷大母线电压（标幺值）
• XΣ = X’d + Xline：系统总电抗（标幺值）
• X’dΣ = X’d + Xline：d轴总暂态电抗（标幺值）
• Xq：q轴同步电抗（标幺值）

3. Python实现单机无穷大系统建模

3.1 环境配置与库导入

首先，我们需要配置Python环境并导入必要的库：

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
import pandas as pd
from matplotlib.font_manager import FontProperties# 设置中文字体
font = FontProperties(fname=r"c:\windows\fonts\simsun.ttc", size=12)# 忽略警告
import warnings
warnings.filterwarnings('ignore')

3.2 系统参数定义

接下来，我们定义单机无穷大系统的参数。这些参数通常来自实际系统的等效化简：

class PowerSystemParameters:"""定义电力系统参数类"""def __init__(self):# 基本参数self.fn = 50.0  # 额定频率(Hz)self.omega_0 = 2 * np.pi * self.fn  # 同步电角速度(rad/s)# 发电机参数self.M = 10.0  # 惯性时间常数(s)self.D = 2.0  # 阻尼系数(pu)self.Xd = 1.8  # d轴同步电抗(pu)self.Xq = 1.7  # q轴同步电抗(pu)self.Xd_prime = 0.3  # d轴暂态电抗(pu)self.Td0_prime = 8.0  # d轴开路暂态时间常数(s)# 励磁系统参数self.Efd = 1.5  # 励磁电动势(pu)# 线路参数self.Xline = 0.3  # 线路电抗(pu)self.Rline = 0.02  # 线路电阻(pu)# 无穷大母线参数self.V_inf = 1.0  # 无穷大母线电压(pu)# 初始运行条件self.Pm0 = 0.9  # 初始机械功率(pu)self.delta0 = None  # 初始功角(rad)，将通过计算得到self.E_prime0 = None  # 初始暂态电动势(pu)，将通过计算得到# 计算总电抗self.X_total = self.Xd_prime + self.Xline  # 系统总电抗(pu)# 计算初始条件self.calculate_initial_conditions()def calculate_initial_conditions(self):"""计算系统初始运行条件"""# 计算初始功角# Pe0 = (E'V_inf/X_total) * sin(delta0) = Pm0# 假设初始运行时E'≈V_inf+X_total*Q0/V_inf，但为简化，我们使用迭代方法# 初始估计delta_guess = np.arcsin(self.Pm0 * self.X_total / (1.0 * self.V_inf))# 使用简单迭代计算初始条件tolerance = 1e-6max_iter = 100E_prime_guess = 1.2  # 初始猜测值for i in range(max_iter):# 计算电磁功率Pe_calculated = (E_prime_guess * self.V_inf / self.X_total) * np.sin(delta_guess)# 计算d轴电流Id = (E_prime_guess - self.V_inf * np.cos(delta_guess)) / self.X_total# 计算新的E_primeE_prime_new = self.V_inf * np.cos(delta_guess) + self.X_total * Id# 检查收敛if abs(Pe_calculated - self.Pm0) < tolerance and abs(E_prime_new - E_prime_guess) < tolerance:break# 更新猜测值delta_guess = np.arcsin(self.Pm0 * self.X_total / (E_prime_guess * self.V_inf))E_prime_guess = E_prime_newself.delta0 = delta_guessself.E_prime0 = E_prime_guessprint(f"初始条件计算完成:")print(f"初始功角: {np.degrees(self.delta0):.2f}°")print(f"初始暂态电动势: {self.E_prime0:.4f} pu")print(f"初始电磁功率: {Pe_calculated:.4f} pu")# 创建系统参数实例
system_params = PowerSystemParameters()

3.3 系统动态方程实现

实现单机无穷大系统的动态方程：

def power_system_dynamics(t, y, params, fault_scenario):"""定义电力系统动态方程y: 状态变量 [delta, omega, E_prime]params: 系统参数对象fault_scenario: 故障场景描述"""delta, omega, E_prime = y# 解包参数omega_0 = params.omega_0M = params.MD = params.DXd = params.XdXd_prime = params.Xd_primeXline = params.XlineTd0_prime = params.Td0_primeEfd = params.EfdV_inf = params.V_infX_total = params.X_totalPm = params.Pm0  # 假设机械功率恒定# 处理故障场景if fault_scenario['type'] == 'none':# 正常情况X_total_effective = X_totalelif fault_scenario['type'] == 'three_phase':# 三相短路故障if t >= fault_scenario['start_time'] and t < fault_scenario['end_time']:X_total_effective = fault_scenario['fault_impedance']  # 故障期间电抗else:X_total_effective = X_total  # 故障清除后else:X_total_effective = X_total# 计算d轴电流Id = (E_prime - V_inf * np.cos(delta)) / X_total_effective# 计算电磁功率Pe = (E_prime * V_inf / X_total_effective) * np.sin(delta)# 计算状态变量的导数ddelta_dt = omega - omega_0domega_dt = (Pm - Pe - D * (omega - omega_0)) / MdE_prime_dt = (Efd - E_prime - (Xd - Xd_prime) * Id) / Td0_primereturn [ddelta_dt, domega_dt, dE_prime_dt]

3.4 故障场景定义

定义不同的故障场景来测试系统的暂态稳定性：

# 定义故障场景
fault_scenarios = {'normal': {'type': 'none','start_time': 0.0,'end_time': 0.0,'fault_impedance': 0.0},'mild_fault': {'type': 'three_phase','start_time': 1.0,'end_time': 1.1,'fault_impedance': 0.1  # 轻微故障，电抗较小},'severe_fault': {'type': 'three_phase','start_time': 1.0,'end_time': 1.2,'fault_impedance': 0.01  # 严重故障，电抗非常小}
}

3.5 仿真执行与结果分析

实现仿真执行和结果分析功能：

def run_simulation(params, fault_scenario, simulation_time=10.0):"""执行暂态稳定仿真"""# 初始状态initial_state = [params.delta0, params.omega_0, params.E_prime0]# 时间点t_eval = np.linspace(0, simulation_time, 1000)# 求解微分方程sol = solve_ivp(fun=lambda t, y: power_system_dynamics(t, y, params, fault_scenario),t_span=[0, simulation_time],y0=initial_state,t_eval=t_eval,method='RK45',rtol=1e-6,atol=1e-8)return soldef analyze_results(sol, params, fault_scenario):"""分析仿真结果"""t = sol.tdelta = sol.y[0]omega = sol.y[1]E_prime = sol.y[2]# 计算电磁功率Pe = np.zeros_like(t)for i in range(len(t)):# 创建临时故障场景用于功率计算temp_fault = fault_scenario.copy()# 确保使用正确的时间点temp_fault['current_time'] = t[i]Pe[i] = (E_prime[i] * params.V_inf / params.X_total) * np.sin(delta[i])# 计算转子相对角速度omega_relative = omega - params.omega_0# 绘制结果plt.figure(figsize=(15, 10))# 功角曲线plt.subplot(2, 2, 1)plt.plot(t, np.degrees(delta))plt.xlabel('时间 (s)')plt.ylabel('功角 (度)')plt.title('发电机功角随时间变化', fontproperties=font)plt.grid(True)# 角速度偏差曲线plt.subplot(2, 2, 2)plt.plot(t, omega_relative)plt.xlabel('时间 (s)')plt.ylabel('角速度偏差 (rad/s)')plt.title('发电机角速度偏差随时间变化', fontproperties=font)plt.grid(True)# 暂态电动势曲线plt.subplot(2, 2, 3)plt.plot(t, E_prime)plt.xlabel('时间 (s)')plt.ylabel('E\' (pu)')plt.title('暂态电动势随时间变化', fontproperties=font)plt.grid(True)# 功率曲线plt.subplot(2, 2, 4)plt.plot(t, Pe, label='电磁功率 Pe')plt.plot(t, [params.Pm0] * len(t), 'r--', label='机械功率 Pm')plt.xlabel('时间 (s)')plt.ylabel('功率 (pu)')plt.title('功率随时间变化', fontproperties=font)plt.legend()plt.grid(True)plt.tight_layout()plt.show()# 稳定性判断delta_final = delta[-100:]  # 最后100个点delta_max = np.max(delta_final)delta_min = np.min(delta_final)delta_variation = delta_max - delta_min# 如果最后阶段的功角波动小于阈值，则认为系统稳定if delta_variation < 0.1:  # 约5.7度stability = "稳定"else:stability = "失稳"print(f"系统暂态稳定性: {stability}")print(f"最大功角: {np.degrees(np.max(delta)):.2f}°")print(f"最小功角: {np.degrees(np.min(delta)):.2f}°")print(f"功角波动范围: {np.degrees(delta_variation):.2f}°")return {'time': t,'delta': delta,'omega': omega,'E_prime': E_prime,'Pe': Pe,'stability': stability}

3.6 完整仿真流程

整合以上组件，形成完整的仿真流程：

def complete_simulation_analysis(params, scenario_name, simulation_time=10.0):"""完整的仿真分析流程"""print(f"正在执行 {scenario_name} 场景仿真...")print(f"故障类型: {fault_scenarios[scenario_name]['type']}")if fault_scenarios[scenario_name]['type'] == 'three_phase':print(f"故障持续时间: {fault_scenarios[scenario_name]['end_time'] - fault_scenarios[scenario_name]['start_time']:.2f} s")# 执行仿真sol = run_simulation(params, fault_scenarios[scenario_name], simulation_time)# 分析结果results = analyze_results(sol, params, fault_scenarios[scenario_name])return results# 执行不同场景的仿真
normal_results = complete_simulation_analysis(system_params, 'normal')
mild_fault_results = complete_simulation_analysis(system_params, 'mild_fault')
severe_fault_results = complete_simulation_analysis(system_params, 'severe_fault')

4. 高级功能扩展

4.1 自动稳定性评估

实现更复杂的稳定性评估算法：

def advanced_stability_analysis(results, params):"""高级稳定性分析"""delta = results['delta']t = results['time']# 寻找第一个极值点peak_indices = []for i in range(1, len(delta)-1):if (delta[i] > delta[i-1] and delta[i] > delta[i+1]) or \(delta[i] < delta[i-1] and delta[i] < delta[i+1]):peak_indices.append(i)# 计算摆动周期和阻尼比if len(peak_indices) >= 3:periods = []dampings = []for i in range(2, len(peak_indices)):T = t[peak_indices[i]] - t[peak_indices[i-2]]periods.append(T)# 计算阻尼比A1 = delta[peak_indices[i-2]]A2 = delta[peak_indices[i]]if A1 != 0:damping_ratio = -np.log(A2/A1) / (2 * np.pi * (i-2))dampings.append(damping_ratio)avg_period = np.mean(periods) if periods else 0avg_damping = np.mean(dampings) if dampings else 0print(f"平均摆动周期: {avg_period:.3f} s")print(f"估计阻尼比: {avg_damping:.4f}")# 基于阻尼比判断稳定性if avg_damping > 0.05:  # 阻尼比大于5%认为稳定stability = "良好阻尼"elif avg_damping > 0.01:  # 阻尼比大于1%认为弱阻尼但稳定stability = "弱阻尼但稳定"else:stability = "负阻尼或零阻尼，可能失稳"print(f"阻尼特性: {stability}")return results# 对每个场景进行高级分析
advanced_stability_analysis(normal_results, system_params)
advanced_stability_analysis(mild_fault_results, system_params)
advanced_stability_analysis(severe_fault_results, system_params)

4.2 参数灵敏度分析

实现参数灵敏度分析功能：

def parameter_sensitivity_analysis(base_params, param_name, values, fault_scenario):"""参数灵敏度分析"""results = {}for value in values:# 创建参数副本modified_params = base_params# 修改参数setattr(modified_params, param_name, value)# 重新计算初始条件modified_params.calculate_initial_conditions()# 运行仿真sol = run_simulation(modified_params, fault_scenario)# 分析稳定性delta = sol.y[0]delta_final = delta[-100:]delta_variation = np.max(delta_final) - np.min(delta_final)# 判断稳定性stable = delta_variation < 0.1results[value] = {'stable': stable,'max_delta': np.max(delta),'min_delta': np.min(delta),'variation': delta_variation}print(f"{param_name} = {value}: {'稳定' if stable else '失稳'}, "f"最大功角: {np.degrees(np.max(delta)):.2f}°")return results# 示例：分析惯性时间常数对稳定性的影响
M_values = [5.0, 7.5, 10.0, 12.5, 15.0]
M_sensitivity = parameter_sensitivity_analysis(system_params, 'M', M_values, fault_scenarios['mild_fault']
)# 可视化灵敏度分析结果
plt.figure(figsize=(10, 6))
stable_M = [M for M in M_values if M_sensitivity[M]['stable']]
unstable_M = [M for M in M_values if not M_sensitivity[M]['stable']]plt.plot(stable_M, [np.degrees(M_sensitivity[M]['max_delta']) for M in stable_M], 'go-', label='稳定', linewidth=2)
plt.plot(unstable_M, [np.degrees(M_sensitivity[M]['max_delta']) for M in unstable_M], 'ro-', label='失稳', linewidth=2)plt.xlabel('惯性时间常数 M (s)')
plt.ylabel('最大功角 (度)')
plt.title('惯性时间常数对暂态稳定性的影响', fontproperties=font)
plt.legend()
plt.grid(True)
plt.show()

4.3 临界切除时间计算

实现临界切除时间（CCT）计算功能：

def calculate_critical_clearing_time(params, fault_scenario_base, max_clearing_time=1.0, tolerance=0.01):"""计算临界切除时间"""# 二分法查找临界切除时间low = 0.0high = max_clearing_timecct = (low + high) / 2.0iterations = 0max_iterations = 20while (high - low) > tolerance and iterations < max_iterations:# 设置当前切除时间current_scenario = fault_scenario_base.copy()current_scenario['end_time'] = cct# 运行仿真sol = run_simulation(params, current_scenario)delta = sol.y[0]# 判断稳定性delta_final = delta[-100:]delta_variation = np.max(delta_final) - np.min(delta_final)stable = delta_variation < 0.1print(f"尝试切除时间: {cct:.3f} s, 系统{'稳定' if stable else '失稳'}")# 调整搜索区间if stable:low = cct  # 如果稳定，尝试更长的切除时间else:high = cct  # 如果失稳，尝试更短的切除时间cct = (low + high) / 2.0iterations += 1print(f"计算完成: 临界切除时间(CCT) = {cct:.3f} s")return cct# 计算严重故障场景的临界切除时间
severe_fault_base = fault_scenarios['severe_fault'].copy()
cct = calculate_critical_clearing_time(system_params, severe_fault_base)

5. 实际应用与工程考虑

5.1 从实际系统到简化模型的参数转换

将实际省级电网参数转换为单机无穷大系统参数是一个重要的工程问题。以下是一些关键考虑：

class ProvincialGridReducer:"""省级电网简化器"""def __init__(self, grid_data):self.grid_data = grid_dataself.base_MVA = grid_data['base_MVA']self.base_kV = grid_data['base_kV']def reduce_to_smib(self, generator_of_interest):"""将省级电网简化为单机无穷大系统generator_of_interest: 所关注的发电机节点"""# 获取关注发电机的参数gen_params = self.grid_data['generators'][generator_of_interest]# 计算等效惯性常数total_inertia = sum(gen['H'] for gen in self.grid_data['generators'].values())equivalent_M = total_inertia  # 简化假设# 计算等效电抗# 这里需要根据网络拓扑进行简化，实际应用中可能需要使用网络等值方法# 如WARD等值或REI等值# 简化假设：使用关注发电机到主网架的阻抗equivalent_X = self.calculate_equivalent_impedance(generator_of_interest)# 创建单机无穷大系统参数smib_params = PowerSystemParameters()smib_params.M = equivalent_Msmib_params.Xline = equivalent_Xreturn smib_paramsdef calculate_equivalent_impedance(self, generator_node):"""计算从发电机节点到无穷大母线的等效阻抗"""# 这里实现简化算法，实际应用中可能需要使用复杂的网络等值算法# 简化假设：使用最短路径电抗# 伪代码：# 1. 找到从发电机到主网架的最短路径# 2. 计算路径上的总电抗# 3. 返回等效电抗return 0.3  # 简化返回值

5.2 模型验证与误差分析

验证简化模型的准确性：

def validate_smib_model(detailed_model_results, smib_model_results):"""验证单机无穷大系统模型的准确性"""# 提取关键指标detailed_max_delta = np.max(np.degrees(detailed_model_results['delta']))smib_max_delta = np.max(np.degrees(smib_model_results['delta']))detailed_settling_time = calculate_settling_time(detailed_model_results['delta'])smib_settling_time = calculate_settling_time(smib_model_results['delta'])# 计算误差delta_error = abs(detailed_max_delta - smib_max_delta) / detailed_max_delta * 100time_error = abs(detailed_settling_time - smib_settling_time) / detailed_settling_time * 100print(f"最大功角误差: {delta_error:.2f}%")print(f"稳定时间误差: {time_error:.2f}%")# 绘制对比图plt.figure(figsize=(12, 6))plt.plot(detailed_model_results['time'], np.degrees(detailed_model_results['delta']), 'b-', label='详细模型', linewidth=2)plt.plot(smib_model_results['time'], np.degrees(smib_model_results['delta']), 'r--', label='SMIB模型', linewidth=2)plt.xlabel('时间 (s)')plt.ylabel('功角 (度)')plt.title('详细模型与SMIB模型对比', fontproperties=font)plt.legend()plt.grid(True)plt.show()return {'delta_error': delta_error,'time_error': time_error}def calculate_settling_time(delta, tolerance=0.01):"""计算稳定时间"""final_value = np.mean(delta[-100:])for i in range(len(delta)):if np.all(np.abs(delta[i:] - final_value) < tolerance):return i / len(delta) * 10  # 假设总仿真时间为10秒return 10.0  # 如果没有稳定，返回总时间

6. 结论与展望

本文详细介绍了使用Python实现电力系统暂态稳定计算中的单机无穷大系统建模方法。通过建立简化的单机无穷大系统模型，我们能够大大提高计算效率，同时保持足够的精度用于工程分析和决策支持。

6.1 主要贡献

1. 完整的数学模型：建立了考虑发电机动态、励磁系统和输电网络的完整数学模型。
2. 灵活的仿真框架：实现了可配置的故障场景和参数设置，支持多种稳定性分析。
3. 高级分析功能：包括参数灵敏度分析、临界切除时间计算和模型验证等功能。
4. 工程应用导向：考虑了从实际系统到简化模型的参数转换问题。

6.2 实际应用价值

单机无穷大系统模型虽然简化，但在电力系统规划和运行中具有重要价值：

1. 快速稳定性评估：在电网规划初期快速评估各种方案的稳定性。
2. 保护整定参考：为继电保护系统提供整定参考值。
3. 教育培训：用于电力系统稳定性原理的教学和培训。
4. 预筛选工具：在详细仿真前预筛选可行方案，提高工作效率。

6.3 未来工作方向

1. 多机系统等值：研究更精确的多机系统等值方法，提高简化模型的准确性。
2. 考虑更详细的设备模型：加入励磁系统、调速系统等更详细的设备模型。
3. 人工智能辅助分析：结合机器学习方法进行快速稳定性评估和控制决策。
4. 云端计算平台：开发基于云计算的暂态稳定分析平台，支持大规模并行计算。

通过持续改进和优化，单机无穷大系统模型将在电力系统暂态稳定分析中发挥更加重要的作用，为电力系统的安全稳定运行提供有力支持。

参考文献

1. Kundur, P., Balu, N. J., & Lauby, M. G. (1994). Power system stability and control. McGraw-Hill.
2. Anderson, P. M., & Fouad, A. A. (2008). Power system control and stability. John Wiley & Sons.
3. Sauer, P. W., & Pai, M. A. (1998). Power system dynamics and stability. Prentice Hall.
4. IEEE Task Force on Benchmark Systems for Stability Controls. (2013). Benchmark systems for small-signal stability analysis and control. IEEE Power & Energy Society.

附录：完整代码结构

power_system_stability/
│
├── models/
│   ├── __init__.py
│   ├── power_system_params.py    # 系统参数定义
│   └── provincial_grid_reducer.py # 电网简化器
│
├── simulation/
│   ├── __init__.py
│   ├── dynamics.py               # 系统动态方程
│   ├── simulator.py              # 仿真执行
│   └── fault_scenarios.py        # 故障场景定义
│
├── analysis/
│   ├── __init__.py
│   ├── stability_analyzer.py     # 稳定性分析
│   ├── sensitivity_analyzer.py   # 灵敏度分析
│   └── cct_calculator.py         # 临界切除时间计算
│
├── validation/
│   ├── __init__.py
│   └── model_validator.py        # 模型验证
│
├── utils/
│   ├── __init__.py
│   ├── plotter.py                # 绘图工具
│   └── data_loader.py            # 数据加载工具
│
└── main.py                       # 主程序

这个完整的Python实现为电力系统暂态稳定计算提供了一个灵活、可扩展的框架，能够有效支持实际工程应用中的决策分析需求。