多场耦合优化-主题090-能源系统与智能电网
·
主题090：能源系统与智能电网
1. 引言
能源系统是现代社会的基础设施，涵盖电力生产、输配、存储和消费的完整链条。随着可再生能源的大规模接入、电动汽车的普及以及分布式能源的发展，传统能源系统正面临着前所未有的挑战。智能电网作为新一代电力系统，通过先进的信息通信技术和自动化控制手段，实现能源的高效、可靠、经济和环保利用。
能源系统仿真是理解和优化复杂能源网络的关键工具。通过建立精确的数学模型和高效的数值算法，我们可以分析系统的稳态和动态行为，预测可再生能源的波动性影响，优化发电调度和负荷分配，评估系统的可靠性和经济性。
本教程将系统介绍能源系统与智能电网的仿真方法，从基础的电力系统潮流计算到复杂的多能源系统优化，从传统的确定性分析到现代的不确定性量化，为读者提供全面的理论基础和实用的仿真技能。
2. 电力系统基础理论
2.1 电力系统组成与结构
电力系统是一个复杂的工程系统，主要由发电、输电、变电、配电和用电五个环节组成。各环节通过电力网络紧密连接，形成统一的电能生产和消费体系。
发电系统包括各种类型的发电机组：
- 火力发电：燃煤、燃气、燃油发电机组，具有调节灵活、响应快速的特点
- 水力发电：径流式、水库式、抽水蓄能电站，运行成本低，调节能力强
- 核能发电：核裂变反应堆，基荷电源，运行稳定
- 可再生能源：风力发电、光伏发电，清洁环保但具有间歇性和波动性
输电系统负责将电能从发电厂远距离输送到负荷中心：
- 交流输电：高压交流（HVAC）输电网络，技术成熟，成本较低
- 直流输电：高压直流（HVDC）输电，适合远距离大容量输电和异步联网
- 输电电压等级：通常包括220kV、330kV、500kV、750kV、1000kV等
变电系统通过变压器实现电压等级的变换：
- 升压变电站：将发电机电压（10-20kV）升至输电电压
- 降压变电站：将输电电压逐级降至配电电压（35kV、10kV）
- 枢纽变电站：连接多个电压等级和电源点的关键节点
配电系统将电能分配到终端用户：
- 高压配电：35kV、10kV配电网，向工商业用户供电
- 低压配电：380V/220V配电网，向居民用户供电
- 配电自动化：实现故障定位、隔离和恢复供电
用电系统包括各类电力负荷：
- 工业负荷：电动机、电弧炉、电解槽等，功率大，负荷特性复杂
- 商业负荷：照明、空调、电梯等，与营业时间相关
- 居民负荷：家用电器，受生活习惯和气候影响
- 电动汽车：新型移动负荷，具有充电灵活性和双向能量流动能力
2.2 电力系统元件模型
2.2.1 同步发电机模型
同步发电机是电力系统中最主要的发电设备，其数学模型包括电气部分和机械部分。
电气部分采用派克变换后的d-q坐标系方程：
电压方程：
v_d = -R_s i_d + dψ_d/dt - ωψ_q
v_q = -R_s i_q + dψ_q/dt + ωψ_d
v_f = R_f i_f + dψ_f/dt
磁链方程：
ψ_d = -L_d i_d + L_{md} i_f
ψ_q = -L_q i_q
ψ_f = -L_{md} i_d + L_f i_f
其中：
- vd,vqv_d, v_qvd,vq：定子d轴、q轴电压
- id,iqi_d, i_qid,iq：定子d轴、q轴电流
- ψd,ψqψ_d, ψ_qψd,ψq：定子d轴、q轴磁链
- RsR_sRs：定子电阻
- Ld,LqL_d, L_qLd,Lq：d轴、q轴同步电感
- LmdL_{md}Lmd：d轴互感
- ωωω：转子电角速度
机械部分采用摇摆方程描述：
2H/ω_s d²δ/dt² = P_m - P_e - D(ω - ω_s)/ω_s
其中：
- HHH：惯性常数（秒）
- δδδ：转子功角（弧度）
- PmP_mPm：机械功率（标幺值）
- PeP_ePe：电磁功率（标幺值）
- DDD：阻尼系数
- ωsω_sωs：同步角速度
2.2.2 变压器模型
变压器是电力系统中实现电压变换的关键设备。在稳态分析中，通常采用Π型等效电路。
双绕组变压器等效电路参数：
Z_T = R_T + jX_T = (P_k/1000S_N) + j(V_k%/100)(V_N²/S_N)
Y_m = G_m - jB_m = (P_0/1000V_N²) - j(I_0%/100)(S_N/V_N²)
其中：
- RTR_TRT：等效电阻（短路损耗）
- XTX_TXT：等效电抗（短路电压百分比）
- GmG_mGm：励磁电导（空载损耗）
- BmB_mBm：励磁电纳（空载电流百分比）
- SNS_NSN：额定容量
- VNV_NVN：额定电压
三绕组变压器采用星形等效电路，需要考虑各绕组之间的阻抗：
Z_{12} = Z_1 + Z_2
Z_{23} = Z_2 + Z_3
Z_{31} = Z_3 + Z_1
解得各绕组等效阻抗：
Z_1 = (Z_{12} + Z_{31} - Z_{23})/2
Z_2 = (Z_{12} + Z_{23} - Z_{31})/2
Z_3 = (Z_{23} + Z_{31} - Z_{12})/2
2.2.3 输电线路模型
输电线路的参数包括电阻、电抗、电导和电纳，分别对应线路的四种物理效应。
线路参数计算：
单位长度电阻：
r = ρ/S (Ω/km)
单位长度电抗：
x = 0.1445lg(D_m/r) + 0.0157 (Ω/km)
单位长度电纳：
b = 7.58/lg(D_m/r) × 10⁻⁶ (S/km)
其中：
- ρρρ：导线电阻率
- SSS：导线截面积
- DmD_mDm：几何均距
- rrr：导线半径
线路等效电路：
短线路（<80km）：串联阻抗模型
Z = (r + jx)l
中长线路（80-250km）：Π型等效电路
Z' = Z sinh(γl)/(γl)
Y'/2 = Y/2 tanh(γl/2)/(γl/2)
长线路（>250km）：分布参数模型
γ = √(zy) = α + jβ
Z_c = √(z/y)
其中传播常数γγγ包括衰减常数ααα和相位常数βββ，特性阻抗ZcZ_cZc表征线路的波阻抗特性。
2.2.4 负荷模型
电力负荷是电力系统的服务对象，其特性对系统运行有重要影响。负荷模型分为静态模型和动态模型。
静态负荷模型：
多项式模型（ZIP模型）：
P = P_0[a_p(V/V_0)² + b_p(V/V_0) + c_p]
Q = Q_0[a_q(V/V_0)² + b_q(V/V_0) + c_q]
指数模型：
P = P_0(V/V_0)^{n_p}
Q = Q_0(V/V_0)^{n_q}
其中：
- P0,Q0P_0, Q_0P0,Q0：额定电压下的有功、无功功率
- V0V_0V0：额定电压
- ap,bp,cpa_p, b_p, c_pap,bp,cp：恒阻抗、恒电流、恒功率负荷比例
- np,nqn_p, n_qnp,nq：电压指数（通常1-3）
动态负荷模型：
感应电动机模型是动态负荷的主要组成部分：
机械运动方程：
T_J dω_r/dt = T_e - T_m
电磁转矩：
T_e = (V²r_r/s)/[(r_s + r_r/s)² + (x_s + x_r)²]
其中：
- TJT_JTJ：惯性时间常数
- ωrω_rωr：转子转速
- Te,TmT_e, T_mTe,Tm：电磁转矩和机械转矩
- sss：转差率
- rs,xsr_s, x_srs,xs：定子电阻、电抗
- rr,xrr_r, x_rrr,xr：转子电阻、电抗（归算到定子侧）
2.3 电力系统稳态分析
2.3.1 潮流计算基础
潮流计算是电力系统分析的基础，用于确定系统在稳态运行条件下的电压、功率分布。
节点功率方程：
对于n节点系统，节点i的注入功率为：
P_i = V_i Σ_{j=1}^n V_j(G_{ij}cosθ_{ij} + B_{ij}sinθ_{ij})
Q_i = V_i Σ_{j=1}^n V_j(G_{ij}sinθ_{ij} - B_{ij}cosθ_{ij})
其中：
- Pi,QiP_i, Q_iPi,Qi：节点i的注入有功、无功功率
- Vi,θiV_i, θ_iVi,θi：节点i的电压幅值和相角
- θij=θi−θjθ_{ij} = θ_i - θ_jθij=θi−θj：节点i、j之间的相角差
- Gij,BijG_{ij}, B_{ij}Gij,Bij：节点导纳矩阵的实部和虚部
节点分类：
| 节点类型 | 已知量 | 未知量 | 典型节点 |
|---|---|---|---|
| PQ节点 | P, Q | V, θ | 负荷节点、联络节点 |
| PV节点 | P, V | Q, θ | 发电机节点、调压节点 |
| 平衡节点 | V, θ | P, Q | 系统参考节点 |
2.3.2 牛顿-拉夫逊法
牛顿-拉夫逊法是求解非线性潮流方程的最常用方法，具有二阶收敛速度。
修正方程：
[ΔP] [H N][Δθ]
[ΔQ] = [J L][ΔV/V]
其中雅可比矩阵元素：
H_{ij} = ∂P_i/∂θ_j = -V_iV_j(G_{ij}sinθ_{ij} - B_{ij}cosθ_{ij})
N_{ij} = ∂P_i/∂V_j V_j = V_iV_j(G_{ij}cosθ_{ij} + B_{ij}sinθ_{ij})
J_{ij} = ∂Q_i/∂θ_j = -V_iV_j(G_{ij}cosθ_{ij} + B_{ij}sinθ_{ij})
L_{ij} = ∂Q_i/∂V_j V_j = -V_iV_j(G_{ij}sinθ_{ij} - B_{ij}cosθ_{ij})
计算步骤：
- 形成节点导纳矩阵Y
- 设定各节点电压初值（通常V=1.0, θ=0）
- 计算功率不平衡量ΔP, ΔQ
- 形成雅可比矩阵J
- 求解修正方程得Δθ, ΔV
- 更新电压：θ^{k+1} = θ^k + Δθ, V^{k+1} = V^k + ΔV
- 检查收敛条件：max|ΔP|, max|ΔQ| < ε
- 若不收敛，返回步骤3
2.3.3 快速解耦法
快速解耦法基于电力系统的物理特性，将有功-电压相角、无功-电压幅值解耦，大幅提高计算效率。
解耦原理：
- 输电网络中，电抗远大于电阻（X>>R）
- 有功功率主要受电压相角影响
- 无功功率主要受电压幅值影响
修正方程：
ΔP/V = B' Δθ
ΔQ/V = B'' ΔV
其中B’和B’'是简化的电纳矩阵：
- B’：忽略线路电阻和并联支路，仅考虑线路电抗
- B’'：忽略线路电阻和相角差的影响
特点：
- 雅可比矩阵为常数矩阵，只需一次因子分解
- 每次迭代只需前代回代求解
- 收敛速度略低于牛顿法，但单次迭代计算量小
- 适合在线计算和大型系统分析
3. 最优潮流与经济调度
3.1 最优潮流问题
最优潮流（Optimal Power Flow, OPF）是在满足各种运行约束的前提下，优化某一目标函数的潮流问题。
目标函数：
最常用的目标是最小化发电成本：
min f = Σ_{i∈G} C_i(P_{Gi}) = Σ_{i∈G} (a_i + b_iP_{Gi} + c_iP_{Gi}²)
其他常见目标：
- 最小化网损：min Σ_{ij} G_{ij}(V_i² + V_j² - 2V_iV_jcosθ_{ij})
- 最小化排放：min Σ_{i∈G} E_i(P_{Gi})
- 最大化社会效益：max Σ_{j∈D} U_j(P_{Dj}) - Σ_{i∈G} C_i(P_{Gi})
约束条件：
等式约束（潮流方程）：
P_{Gi} - P_{Di} = V_i Σ_{j=1}^n V_j(G_{ij}cosθ_{ij} + B_{ij}sinθ_{ij})
Q_{Gi} - Q_{Di} = V_i Σ_{j=1}^n V_j(G_{ij}sinθ_{ij} - B_{ij}cosθ_{ij})
不等式约束：
P_{Gi}^{min} ≤ P_{Gi} ≤ P_{Gi}^{max}
Q_{Gi}^{min} ≤ Q_{Gi} ≤ Q_{Gi}^{max}
V_i^{min} ≤ V_i ≤ V_i^{max}
|S_{ij}| ≤ S_{ij}^{max}
3.2 经典经济调度
经济调度（Economic Dispatch, ED）是在不考虑网络约束的情况下，优化发电机组出力的分配。
等微增率准则：
对于忽略网损和线路约束的情况，最优条件为各机组的边际成本相等：
dC_1/dP_{G1} = dC_2/dP_{G2} = ... = dC_n/dP_{Gn} = λ
其中λ为系统边际成本（影子价格）。
考虑网损的协调方程：
dC_i/dP_{Gi} × 1/(1 - ∂P_L/∂P_{Gi}) = λ
其中∂P_L/∂P_{Gi}为网损微增率，1/(1 - ∂P_L/∂P_{Gi})为网损修正系数。
λ迭代法：
- 设定λ的初值
- 根据等微增率准则计算各机组出力
- 检查出力是否越限，修正越限机组
- 计算总发电功率，检查功率平衡
- 若不满足，调整λ并返回步骤2
3.3 机组组合问题
机组组合（Unit Commitment, UC）是确定未来一段时间内各发电机组的开停机计划，是一个混合整数规划问题。
目标函数：
min Σ_{t=1}^T Σ_{i=1}^N [C_i(P_{Gi,t})u_{i,t} + S_{i,t}(1-u_{i,t-1})u_{i,t}]
其中：
- ui,tu_{i,t}ui,t：机组i在时段t的状态（0-停机，1-运行）
- Si,tS_{i,t}Si,t：机组i在时段t的启动成本
约束条件：
系统功率平衡：
Σ_{i=1}^N P_{Gi,t}u_{i,t} = D_t + P_{L,t}
旋转备用：
Σ_{i=1}^N min(P_{Gi}^{max} - P_{Gi,t}, R_i^{max})u_{i,t} ≥ R_t^{req}
最小开停机时间：
(u_{i,t} - u_{i,t-1})T_i^{on} ≤ Σ_{k=t}^{min(t+T_i^{on}-1,T)} u_{i,k}
(u_{i,t-1} - u_{i,t})T_i^{off} ≤ Σ_{k=t}^{min(t+T_i^{off}-1,T)} (1 - u_{i,k})
求解方法：
- 优先级列表法：按平均成本排序，简单快速
- 动态规划法：考虑时间耦合，适合小规模问题
- 拉格朗日松弛法：分解为单机组子问题，适合大规模系统
- 混合整数规划：精确求解，适合考虑复杂约束
4. 可再生能源与分布式能源
4.1 风力发电建模
风力发电是将风能转换为电能的技术，具有清洁、可再生的特点，但也存在间歇性和波动性。
风力机功率特性：
风力机的输出功率与风速的关系：
P_w =
{
0, v < v_{ci} 或 v > v_{co}
0.5ρAC_p(λ,β)v³, v_{ci} ≤ v ≤ v_r
P_r, v_r < v ≤ v_{co}
}
其中：
- vci,vcov_{ci}, v_{co}vci,vco：切入风速、切出风速（通常3-4m/s和20-25m/s）
- vrv_rvr：额定风速（通常12-15m/s）
- PrP_rPr：额定功率
- ρρρ：空气密度
- AAA：风轮扫风面积
- CpC_pCp：风能利用系数（理论上限为Betz极限0.593）
风速模型：
Weibull分布描述风速的统计特性：
f(v) = (k/c)(v/c)^{k-1}exp[-(v/c)^k]
其中：
- kkk：形状参数（通常1.5-3）
- ccc：尺度参数，与平均风速相关
风电场聚合模型：
对于大型风电场，需要考虑尾流效应和电气集电系统：
P_{wf} = Σ_{i=1}^{N_w} P_{w,i} × (1 - δ_{wake,i}) × (1 - δ_{elec,i})
其中：
- δwakeδ_{wake}δwake：尾流损失（通常5-15%）
- δelecδ_{elec}δelec：电气损失（通常2-5%）
4.2 光伏发电建模
光伏发电利用半导体材料的光生伏特效应将太阳能直接转换为电能。
光伏电池模型：
单二极管模型：
I = I_{ph} - I_0[exp((V + IR_s)/(nV_t)) - 1] - (V + IR_s)/R_{sh}
其中：
- IphI_{ph}Iph：光生电流，与辐照度成正比
- I0I_0I0：反向饱和电流
- RsR_sRs：串联电阻
- RshR_{sh}Rsh：并联电阻
- nnn：理想因子
- Vt=kT/qV_t = kT/qVt=kT/q：热电压
功率输出：
光伏阵列的输出功率：
P_{pv} = η_{mppt} × η_{inv} × N_s × N_p × P_{module}
其中：
- ηmpptη_{mppt}ηmppt：最大功率点跟踪效率（通常95-99%）
- ηinvη_{inv}ηinv：逆变器效率（通常95-98%）
- Ns,NpN_s, N_pNs,Np：串联、并联组件数
太阳辐射模型：
晴空辐射模型（Clear Sky Model）：
G_{cs} = G_{sc}(1 + 0.033cos(360n/365)) × (cosφcosδcosω + sinφsinδ)
其中：
- GscG_{sc}Gsc：太阳常数（1367 W/m²）
- nnn：一年中的第几天
- φφφ：当地纬度
- δδδ：太阳赤纬角
- ωωω：时角
4.3 储能系统建模
储能系统是平衡可再生能源波动、提高电网灵活性的关键技术。
电池储能模型：
状态 of Charge（SOC）动态：
SOC(t+1) = SOC(t) + η_c P_c(t)Δt/E_{cap} - P_d(t)Δt/(η_d E_{cap})
约束条件：
SOC_{min} ≤ SOC(t) ≤ SOC_{max}
0 ≤ P_c(t) ≤ P_{c}^{max}
0 ≤ P_d(t) ≤ P_{d}^{max}
P_c(t) × P_d(t) = 0 （不能同时充放电）
其中：
- ηc,ηdη_c, η_dηc,ηd：充放电效率
- EcapE_{cap}Ecap：额定容量
- Pc,PdP_c, P_dPc,Pd：充放电功率
储能应用场景：
- 能量时移（Energy Arbitrage）：在电价低时充电，电价高时放电
- 调频服务（Frequency Regulation）：快速响应频率偏差
- 可再生能源平滑：平抑风电、光伏的功率波动
- 备用容量：提供旋转备用或替代备用
- 延缓电网升级：在负荷增长区域替代输配电扩容
4.4 电动汽车与V2G
电动汽车（EV）作为新型移动负荷和分布式储能，为电网运行带来新的机遇和挑战。
充电负荷建模：
电动汽车充电需求：
P_{ev}(t) = Σ_{i=1}^{N_{ev}} P_{ch,i}(t) × x_i(t)
其中：
- xi(t)x_i(t)xi(t)：车辆i在t时刻的充电状态（0/1）
- Pch,iP_{ch,i}Pch,i：充电功率
V2G（Vehicle-to-Grid）：
电动汽车向电网反馈电能：
P_{v2g}(t) = Σ_{i=1}^{N_{v2g}} P_{dis,i}(t) × y_i(t)
V2G的优势：
- 提供调频、调峰等辅助服务
- 提高可再生能源消纳能力
- 降低车主用电成本
- 延缓电网投资
智能充电策略：
- 分时电价引导：根据电价信号调整充电时间
- 需求响应：在电网高峰时段降低充电功率
- 车队优化：协调大规模电动汽车的充放电
- 与可再生能源协调：在风电、光伏大发时段充电
5. 智能电网技术
5.1 高级量测体系
高级量测体系（Advanced Metering Infrastructure, AMI）是智能电网的数据基础，包括智能电表、通信网络和数据分析系统。
智能电表功能：
- 双向通信：远程抄表、费率下发、控制指令
- 高精度测量：电压、电流、功率、功率因数
- 事件记录：停电、电压暂降、窃电检测
- 负荷控制：远程通断、需求响应
数据分析应用：
负荷预测：
P_{load}(t) = f(T(t), H(t), DOW(t), HOD(t), Hist(t))
非侵入式负荷监测（NILM）：
P_{total}(t) = Σ_{i=1}^N P_i(t) × s_i(t)
其中si(t)s_i(t)si(t)为设备i在t时刻的状态。
5.2 需求响应
需求响应（Demand Response, DR）是通过价格信号或激励措施引导用户调整用电行为，实现供需平衡。
价格型需求响应：
分时电价（Time-of-Use, TOU）：
π(t) = π_{base} × [1 + α_{peak}I_{peak}(t) + α_{valley}I_{valley}(t)]
实时电价（Real-Time Pricing, RTP）：
π(t) = λ(t) = dC_{total}/dD(t)
激励型需求响应：
直接负荷控制（Direct Load Control, DLC）：
P_{load}(t) = P_{base}(t) × (1 - u_{ctrl}(t) × α_{red})
可中断负荷（Interruptible Load, IL）：
max U = Σ_t [π(t)P(t) - C_{curt}(t)]
s.t. P(t) ≥ P_{min}(t)
5.3 微电网与分布式能源
微电网是由分布式电源、储能、负荷和控制系统组成的局部电力系统，可以并网或孤岛运行。
微电网结构：
交流微电网：
- 各DG通过逆变器接入交流母线
- 传统电力设备兼容性好
- 需要同步控制
直流微电网：
- 光伏、储能、电动汽车等直流设备直接接入
- 减少变换环节，提高效率
- 控制相对简单
交直流混合微电网：
- 结合交流和直流的优点
- 根据负荷类型选择供电方式
- 灵活性高
微电网控制策略：
主从控制（Master-Slave）：
f = f_0 - k_p(P - P_0)
V = V_0 - k_q(Q - Q_0)
对等控制（Peer-to-Peer）：
f_i = f_0 - m_i(P_i - P_i^{ref})
V_i = V_0 - n_i(Q_i - Q_i^{ref})
下垂控制（Droop Control）：
f = f_{max} - (f_{max} - f_{min})/(P_{max} - P_{min}) × (P - P_{min})
5.4 主动配电网
主动配电网（Active Distribution Network, ADN）是具有主动控制和优化能力的现代配电网，能够大规模接纳分布式能源。
主动配电网特征：
- 双向潮流：从传统的单向供电变为双向能量流动
- 主动控制：通过无功补偿、电压调节等手段主动管理网络
- 信息集成：实时监测和通信系统支撑智能决策
- 市场参与：分布式能源参与电力市场交易
电压控制策略：
协调电压控制：
min Σ_{i∈N} (V_i - V_i^{ref})²
s.t. V_i^{min} ≤ V_i ≤ V_i^{max}
Q_{DG,i}^{min} ≤ Q_{DG,i} ≤ Q_{DG,i}^{max}
Q_{CB,j} ∈ {0, Q_{step}, 2Q_{step}, ...}
其中：
- QDGQ_{DG}QDG：分布式电源的无功出力
- QCBQ_{CB}QCB：电容器组的无功补偿
6. 电力市场与经济运行
6.1 电力市场结构
电力市场是通过市场机制实现电力资源优化配置的制度安排，包括电能市场、辅助服务市场和容量市场。
市场类型：
-
现货市场（Spot Market）：
- 日前市场（Day-Ahead）：提前一天确定发电计划和电价
- 实时市场（Real-Time）：平衡实际供需偏差
-
辅助服务市场（Ancillary Service Market）：
- 调频服务：一次调频、二次调频、三次调频
- 备用服务：旋转备用、非旋转备用、替代备用
- 无功支持：电压调节、无功补偿
-
容量市场（Capacity Market）：
- 确保长期供电可靠性
- 回收固定成本
- 激励新建发电容量
市场出清模型：
社会福利最大化：
max Σ_{j∈D} B_j(P_{Dj}) - Σ_{i∈G} C_i(P_{Gi})
s.t. Σ_{i∈G} P_{Gi} = Σ_{j∈D} P_{Dj}
P_{Gi}^{min} ≤ P_{Gi} ≤ P_{Gi}^{max}
P_{Dj}^{min} ≤ P_{Dj} ≤ P_{Dj}^{max}
其中：
- BjB_jBj：负荷j的效用函数（购电意愿）
- CiC_iCi：机组i的成本函数
6.2 节点电价
节点电价（Locational Marginal Price, LMP）反映了在特定节点增加单位负荷的边际成本，包括电能成本、阻塞成本和网损成本。
LMP分解：
LMP_i = λ_{sys} + λ_{congestion,i} + λ_{loss,i}
其中：
- λsysλ_{sys}λsys：系统边际成本（参考节点LMP）
- λcongestion,iλ_{congestion,i}λcongestion,i：节点i的阻塞分量
- λloss,iλ_{loss,i}λloss,i：节点i的网损分量
LMP计算：
通过对偶变量获得：
LMP_i = ∂L/∂P_{Di} = μ_{balance} + Σ_{k∈K} μ_k × SF_{i,k}
其中：
- μbalanceμ_{balance}μbalance：功率平衡约束的影子价格
- μkμ_kμk：线路k约束的影子价格
- SFi,kSF_{i,k}SFi,k：节点i对线路k的功率转移分布因子
6.3 可再生能源参与市场
随着可再生能源比例提高，其参与电力市场的方式也在不断创新。
可再生能源的出力特性：
- 零边际成本：风、光资源免费
- 间歇性和波动性：受天气条件影响
- 预测不确定性：日前预测误差通常5-15%
市场参与模式：
-
必须运行（Must-Take）：
- 可再生能源优先消纳
- 不参与竞价，按预测出力上网
- 通过补贴或溢价机制补偿
-
参与竞价（Price-Taker）：
- 按零报价参与市场
- 按市场出清价结算
- 承担预测偏差惩罚
-
自调度（Self-Scheduling）：
- 与储能、需求响应联合优化
- 参与能量和辅助服务市场
- 最大化整体收益
可再生能源的不确定性管理：
场景分析法：
min Σ_{s=1}^{N_s} π_s [Σ_{i∈G} C_i(P_{Gi,s}) + C_{curt}(P_{curt,s}) + C_{short}(P_{short,s})]
其中：
- πsπ_sπs：场景s的概率
- PcurtP_{curt}Pcurt：弃风/弃光功率
- PshortP_{short}Pshort：供电不足功率
7. 电力系统动态与稳定性
7.1 暂态稳定性
暂态稳定性是指系统在遭受大扰动（如短路故障、线路跳闸）后，各同步电机保持同步运行的能力。
摇摆方程：
M_i d²δ_i/dt² = P_{mi} - P_{ei} - D_i dδ_i/dt
其中：
- Mi=2Hi/ωsM_i = 2H_i/ω_sMi=2Hi/ωs：惯性常数
- Pei=Σj=1nEiEj(Gijcosδij+Bijsinδij)P_{ei} = Σ_{j=1}^n E_iE_j(G_{ij}cosδ_{ij} + B_{ij}sinδ_{ij})Pei=Σj=1nEiEj(Gijcosδij+Bijsinδij)：电磁功率
等面积法则：
对于单机无穷大系统，暂态稳定判据：
A_{acc} = ∫_{δ_0}^{δ_c} (P_m - P_e^{fault})dδ （加速面积）
A_{dec} = ∫_{δ_c}^{δ_{max}} (P_e^{post} - P_m)dδ （减速面积）
稳定条件：A_{acc} ≤ A_{dec}
暂态稳定评估方法：
-
时域仿真法：
- 数值积分摇摆方程
- 考虑详细模型和复杂故障
- 计算量大，但精度高
-
直接法：
- 构造李雅普诺夫函数
- 快速评估稳定裕度
- 模型简化，保守性较强
-
混合法：
- 结合时域仿真和直接法
- 先用直接法筛选严重故障
- 再对关键故障进行时域仿真
7.2 小干扰稳定性
小干扰稳定性是指系统在遭受小扰动（如负荷波动）后，恢复到原运行状态的能力。
线性化模型：
在平衡点线性化：
Δẋ = AΔx + BΔu
Δy = CΔx + DΔu
其中状态变量包括各发电机的功角、转速、励磁电压等。
特征值分析：
系统矩阵A的特征值决定稳定性：
|A - λI| = 0
- 实部为负：稳定
- 实部为正：不稳定
- 实部为零：临界稳定
机电振荡模式：
-
局部模式（Local Mode）：
- 频率0.7-2.0 Hz
- 一台发电机对系统其余部分的振荡
-
区域间模式（Inter-Area Mode）：
- 频率0.1-0.7 Hz
- 一个区域的发电机群对另一个区域的振荡
-
控制模式（Control Mode）：
- 与励磁系统、调速系统相关
- 频率0.2-3.0 Hz
阻尼控制：
电力系统稳定器（PSS）：
V_{PSS} = K_{PSS} × [ω(s)/(1 + sT_W)] × [(1 + sT_1)(1 + sT_3)]/[(1 + sT_2)(1 + sT_4)]
通过提供附加阻尼转矩，抑制低频振荡。
7.3 电压稳定性
电压稳定性是指系统在扰动后维持所有节点电压在可接受范围内的能力。
电压崩溃机理：
负荷持续增长导致：
- 无功需求增加
- 输电网络无功损耗增加
- 发电机无功出力达到极限
- 电压急剧下降，形成崩溃
PV曲线分析：
对于简单系统：
V = √[E²/2 - QX ± √(E⁴/4 - X²P² - XE²Q)]
最大功率传输点（鼻点）：
P_{max} = E²/(2X)
电压稳定指标：
-
灵敏度指标：
dV/dQ, dQ/dV -
奇异值指标：
σ_{min}(J_{QV}) -
裕度指标：
MW Margin = P_{nose} - P_{operating}
电压控制措施：
- 发电机无功调节
- 电容器组投切
- 有载调压变压器（OLTC）调节
- 静止无功补偿器（SVC）
- 静止同步补偿器（STATCOM）
- 低压减载（UVLS）
8. 多能源系统综合优化
8.1 电-热-气综合能源系统
综合能源系统（Integrated Energy System, IES）将电力、热力、燃气等多种能源形式协同优化，提高能源利用效率。
能源枢纽模型：
能源枢纽（Energy Hub）是IES的基本单元：
[P_e] [η_{CHP,e} 0 η_{HP} ][P_{CHP}]
[P_h] = [η_{CHP,h} η_{GB} η_{HP} ][P_{GB} ]
[P_g] [-1 -1 0 ][P_{HP} ]
其中：
- CHP：热电联产机组
- GB：燃气锅炉
- HP：热泵
电-热联合优化：
目标函数：
min Σ_{t=1}^T [C_{elec}(t) + C_{gas}(t) + C_{OM}(t)]
约束条件：
P_{load}(t) = P_{grid}(t) + P_{CHP,e}(t) + P_{HP,e}(t)
Q_{load}(t) = Q_{CHP,h}(t) + Q_{GB}(t) + Q_{HP,h}(t) - Q_{storage}(t)
Q_{storage}(t+1) = Q_{storage}(t) + η_{ch}Q_{ch}(t) - Q_{dis}(t)/η_{dis}
热电联产（CHP）：
CHP同时产生电能和热能，能源利用效率高：
η_{total} = (P_e + Q_h)/P_{fuel}
典型效率可达80-90%，远高于分产（发电效率35-50%，锅炉效率80-90%）。
8.2 氢能系统
氢能作为一种清洁、高效的二次能源，在综合能源系统中发挥重要作用。
电解水制氢：
电解槽功率消耗：
P_{elec} = m_{H2} × HHV/η_{elec}
其中：
- mH2m_{H2}mH2：氢气质量流量
- HHVHHVHHV：氢气高热值（141.9 MJ/kg）
- ηelecη_{elec}ηelec：电解效率（通常60-80%）
氢燃料电池：
燃料电池发电：
P_{FC} = m_{H2} × LHV × η_{FC}
其中：
- LHVLHVLHV：氢气低热值（120.1 MJ/kg）
- ηFCη_{FC}ηFC：燃料电池效率（通常40-60%）
氢能储能：
氢能作为长周期储能：
E_{storage} = m_{H2} × HHV
优势：
- 储能容量大，适合季节性储能
- 能量密度高
- 可跨能源部门使用（电力、交通、工业）
8.3 碳捕集与利用
碳捕集、利用与封存（CCUS）是实现碳中和的重要技术手段。
碳捕集电厂：
碳捕集能耗：
P_{CCS} = α × m_{CO2}
其中ααα为捕集能耗系数（通常0.2-0.4 MWh/tCO2）。
净发电功率：
P_{net} = P_{gross} - P_{CCS}
碳捕集与可再生能源协同：
利用可再生能源为碳捕集提供能量：
P_{CCS}(t) = P_{RE}(t) × η_{CCS}
实现负排放：
E_{net} = E_{biomass} - E_{CCS} < 0
9. 不确定性分析与鲁棒优化
9.1 可再生能源不确定性建模
可再生能源出力的不确定性是电力系统规划和运行的主要挑战。
概率预测：
点预测：
P̂_{RE}(t+k|t) = E[P_{RE}(t+k)|信息t]
区间预测：
P[P_{RE}^{min}(t+k) ≤ P_{RE}(t+k) ≤ P_{RE}^{max}(t+k)] = 1-α
概率密度预测：
f_{P_{RE}}(p; t+k|t)
场景生成方法：
蒙特卡洛模拟：
P_{RE,s} = P̂_{RE} + ε_s, ε_s ~ N(0, σ²)
拉丁超立方采样（LHS）：
分层采样，提高采样效率
场景削减：
min Σ_{s=1}^{N_{red}} π_s d(s, s̄)
9.2 鲁棒优化方法
鲁棒优化考虑最坏情况，保证解在不确定性范围内的可行性。
两阶段鲁棒优化：
min_x c^T x + max_{u∈U} min_y d^T y
s.t. Ax + By ≥ h - Cu, ∀u∈U
其中：
- xxx：第一阶段决策（ here-and-now）
- yyy：第二阶段决策（wait-and-see）
- uuu：不确定参数
- UUU：不确定性集合
不确定性集合：
盒式不确定性集合：
U = {u | u_i^{min} ≤ u_i ≤ u_i^{max}, ∀i}
多面体不确定性集合：
U = {u | Au ≤ b}
椭球不确定性集合：
U = {u | (u-ū)^T Σ^{-1} (u-ū) ≤ Γ²}
9.3 随机规划方法
随机规划通过场景描述不确定性，优化期望目标。
两阶段随机规划：
min_x c^T x + E_ξ[Q(x, ξ)]
s.t. Ax ≥ b
其中：
- Q(x,ξ)=minydTy∣Wy≥h(ξ)−T(ξ)xQ(x, ξ) = min_y {d^T y | Wy ≥ h(ξ) - T(ξ)x}Q(x,ξ)=minydTy∣Wy≥h(ξ)−T(ξ)x
- ξξξ：随机变量
场景法：
min c^T x + Σ_{s=1}^{N_s} π_s d^T y_s
s.t. Ax ≥ b
Wy_s ≥ h_s - T_s x, ∀s
机会约束规划：
min f(x)
s.t. P[Ax ≥ b(ξ)] ≥ 1-ε
其中εεε为允许违反约束的概率。
10. 仿真案例与Python实现
10.1 案例1：IEEE 14节点系统潮流计算
IEEE 14节点系统是电力系统分析的标准测试系统，包含5台发电机、14条母线和20条支路。
系统参数：
节点数据：
| 节点 | 类型 | 电压幅值 | 电压相角 | 有功负荷 | 无功负荷 |
|---|---|---|---|---|---|
| 1 | 平衡 | 1.060 | 0.000 | 0.000 | 0.000 |
| 2 | PV | 1.045 | - | 21.700 | 12.700 |
| 3 | PV | 1.010 | - | 94.200 | 19.000 |
| … | … | … | … | … | … |
支路数据：
| 从节点 | 到节点 | 电阻 | 电抗 | 电纳 |
|---|---|---|---|---|
| 1 | 2 | 0.01938 | 0.05917 | 0.02640 |
| 1 | 5 | 0.05403 | 0.22304 | 0.02460 |
| … | … | … | … | … |
Python实现：
import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse import csr_matrix
from scipy.sparse.linalg import spsolve
class PowerFlow:
def __init__(self, n_bus):
self.n_bus = n_bus
self.Y = np.zeros((n_bus, n_bus), dtype=complex)
self.V = np.ones(n_bus, dtype=complex)
self.P = np.zeros(n_bus)
self.Q = np.zeros(n_bus)
self.bus_type = np.ones(n_bus) # 1:PQ, 2:PV, 3:Slack
def add_line(self, i, j, r, x, b):
"""添加输电线路"""
z = complex(r, x)
y = 1/z
self.Y[i,j] -= y
self.Y[j,i] -= y
self.Y[i,i] += y + complex(0, b/2)
self.Y[j,j] += y + complex(0, b/2)
def newton_raphson(self, tol=1e-6, max_iter=20):
"""牛顿-拉夫逊法求解潮流"""
for iteration in range(max_iter):
# 计算功率不平衡
dP, dQ = self.power_mismatch()
# 检查收敛
if max(np.max(np.abs(dP)), np.max(np.abs(dQ))) < tol:
print(f"收敛于第{iteration+1}次迭代")
return True
# 形成雅可比矩阵
J = self.jacobian()
# 求解修正方程
dx = np.linalg.solve(J, np.concatenate([dP, dQ]))
# 更新电压
n_pq = np.sum(self.bus_type == 1)
dtheta = dx[:self.n_bus-1]
dV = dx[self.n_bus-1:]
self.V[1:] *= np.exp(1j * dtheta)
self.V[self.bus_type == 1] *= (1 + dV)
print("达到最大迭代次数，未收敛")
return False
def power_mismatch(self):
"""计算功率不平衡量"""
S_calc = self.V * np.conj(self.Y @ self.V)
P_calc = np.real(S_calc)
Q_calc = np.imag(S_calc)
dP = self.P - P_calc
dQ = self.Q - Q_calc
# 排除平衡节点和PV节点的无功
dP = dP[1:] # 排除平衡节点
dQ = dQ[self.bus_type == 1] # 仅PQ节点
return dP, dQ
def jacobian(self):
"""形成雅可比矩阵"""
n = self.n_bus
n_pq = np.sum(self.bus_type == 1)
J = np.zeros((n-1+n_pq, n-1+n_pq))
# 计算雅可比矩阵元素
for i in range(1, n):
for j in range(1, n):
if i == j:
J[i-1, j-1] = -np.imag(self.V[i] * np.conj(self.Y[i,:] @ self.V))
else:
J[i-1, j-1] = np.abs(self.V[i]) * np.abs(self.V[j]) * \
(self.Y[i,j].real * np.sin(np.angle(self.V[i]) - np.angle(self.V[j])) - \
self.Y[i,j].imag * np.cos(np.angle(self.V[i]) - np.angle(self.V[j])))
return J
# 创建IEEE 14节点系统
pf = PowerFlow(14)
# 添加支路数据（简化示例）
pf.add_line(0, 1, 0.01938, 0.05917, 0.0528)
pf.add_line(0, 4, 0.05403, 0.22304, 0.0492)
# ... 添加其他支路
# 设置节点数据
pf.P = np.array([0, 21.7, 94.2, 47.8, 7.6, 11.2, 0, 0, 29.5, 9, 3.5, 6.1, 13.5, 14.9]) / 100
pf.Q = np.array([0, 12.7, 19.0, -3.9, 1.6, 7.5, 0, 0, 16.6, 5.8, 1.8, 1.6, 5.8, 5.0]) / 100
# 设置节点类型
pf.bus_type[0] = 3 # 平衡节点
pf.bus_type[1] = 2 # PV节点
pf.bus_type[2] = 2 # PV节点
# 求解潮流
pf.newton_raphson()
# 输出结果
print("节点电压:")
for i in range(14):
print(f"节点{i+1}: |V|={np.abs(pf.V[i]):.4f}, θ={np.angle(pf.V[i])*180/np.pi:.2f}°")
10.2 案例2：经济调度优化
考虑3台发电机组的经济调度问题，优化各机组出力以最小化总发电成本。
机组成本函数：
| 机组 | a ($/MW²h) | b ($/MWh) | c ($/h) | Pmin (MW) | Pmax (MW) |
|---|---|---|---|---|---|
| G1 | 0.001562 | 7.92 | 561 | 150 | 600 |
| G2 | 0.001940 | 7.85 | 310 | 100 | 400 |
| G3 | 0.004820 | 7.97 | 78 | 50 | 200 |
Python实现：
from scipy.optimize import minimize
import numpy as np
# 机组参数
a = np.array([0.001562, 0.001940, 0.004820])
b = np.array([7.92, 7.85, 7.97])
c = np.array([561, 310, 78])
P_min = np.array([150, 100, 50])
P_max = np.array([600, 400, 200])
# 负荷需求
P_load = 850 # MW
def cost_function(P):
"""发电成本函数"""
return np.sum(a * P**2 + b * P + c)
def power_balance(P):
"""功率平衡约束"""
return np.sum(P) - P_load
# 等微增率法求解
lambda_min = max(b)
lambda_max = max(2 * a * P_max + b)
def solve_ed(lambda_val):
"""给定λ，求解各机组出力"""
P = (lambda_val - b) / (2 * a)
# 处理越限
P = np.maximum(P, P_min)
P = np.minimum(P, P_max)
return P
# λ迭代
for _ in range(100):
lambda_mid = (lambda_min + lambda_max) / 2
P = solve_ed(lambda_mid)
if np.sum(P) > P_load:
lambda_max = lambda_mid
else:
lambda_min = lambda_mid
if abs(np.sum(P) - P_load) < 0.01:
break
print("经济调度结果:")
print(f"系统边际成本λ = {lambda_mid:.4f} $/MWh")
for i in range(3):
print(f"机组{i+1}: P = {P[i]:.2f} MW, 边际成本 = {2*a[i]*P[i] + b[i]:.4f} $/MWh")
print(f"总发电成本 = {cost_function(P):.2f} $/h")
10.3 案例3：风电功率预测
利用时间序列方法预测风电功率，为电力系统调度提供依据。
Python实现：
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
# 生成模拟风电数据
np.random.seed(42)
n_hours = 168 # 一周数据
t = np.arange(n_hours)
# 风速（Weibull分布 + 日周期 + 趋势）
k, c = 2.0, 8.0 # Weibull参数
base_wind = c * np.random.weibull(k, n_hours)
daily_cycle = 2 * np.sin(2 * np.pi * t / 24)
trend = 0.01 * t
wind_speed = base_wind + daily_cycle + trend
wind_speed = np.maximum(wind_speed, 0)
# 风电功率（简化模型）
v_ci, v_r, v_co = 3, 12, 25 # 切入、额定、切出风速
P_r = 100 # 额定功率
P_wind = np.zeros(n_hours)
for i in range(n_hours):
v = wind_speed[i]
if v < v_ci or v > v_co:
P_wind[i] = 0
elif v < v_r:
P_wind[i] = P_r * ((v - v_ci)/(v_r - v_ci))**3
else:
P_wind[i] = P_r
# 简单预测模型（持续性模型 + 趋势）
def persistence_forecast(P_hist, horizon=24):
"""持续性预测"""
P_pred = np.zeros(horizon)
for h in range(horizon):
if h < len(P_hist):
P_pred[h] = P_hist[-1] # 用最新值预测
else:
P_pred[h] = P_pred[h-1]
return P_pred
# 预测未来24小时
P_forecast = persistence_forecast(P_wind, 24)
# 可视化
plt.figure(figsize=(14, 8))
plt.subplot(2, 1, 1)
plt.plot(t, wind_speed, 'b-', label='Wind Speed')
plt.axhline(y=v_ci, color='g', linestyle='--', label='Cut-in')
plt.axhline(y=v_r, color='r', linestyle='--', label='Rated')
plt.axhline(y=v_co, color='k', linestyle='--', label='Cut-out')
plt.xlabel('Hour')
plt.ylabel('Wind Speed (m/s)')
plt.title('Wind Speed Time Series')
plt.legend()
plt.grid(True, alpha=0.3)
plt.subplot(2, 1, 2)
plt.plot(t, P_wind, 'b-', label='Actual')
plt.plot(np.arange(n_hours, n_hours+24), P_forecast, 'r--', label='Forecast')
plt.xlabel('Hour')
plt.ylabel('Wind Power (MW)')
plt.title('Wind Power and Forecast')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('wind_power_forecast.png', dpi=150)
plt.show()
10.4 案例4：微电网能量管理
考虑包含光伏、储能、负荷的微电网，优化储能充放电策略。
Python实现：
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import linprog
# 24小时数据
hours = 24
# 光伏出力（归一化）
pv_profile = np.array([0, 0, 0, 0, 0, 0.05, 0.15, 0.35, 0.6, 0.75,
0.85, 0.9, 0.88, 0.8, 0.65, 0.45, 0.25, 0.1,
0.02, 0, 0, 0, 0, 0])
P_pv_max = 100 # kW
P_pv = P_pv_max * pv_profile
# 负荷曲线
load_profile = np.array([0.4, 0.35, 0.33, 0.32, 0.35, 0.45, 0.6, 0.75,
0.8, 0.78, 0.75, 0.72, 0.7, 0.72, 0.75, 0.8,
0.85, 0.9, 0.88, 0.75, 0.65, 0.55, 0.48, 0.42])
P_load_max = 150 # kW
P_load = P_load_max * load_profile
# 电价（分时电价）
price = np.array([0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.5, 0.5, 0.8, 0.8,
0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.5, 0.5, 0.8, 0.8,
0.8, 0.5, 0.5, 0.3]) # $/kWh
# 储能参数
E_max = 200 # kWh
P_ch_max = 50 # kW
P_dis_max = 50 # kW
eta_ch = 0.95 # 充电效率
eta_dis = 0.95 # 放电效率
SOC_min = 0.1
SOC_max = 0.9
SOC_0 = 0.5 # 初始SOC
# 简化优化：贪心算法
SOC = np.zeros(hours + 1)
SOC[0] = SOC_0
P_ch = np.zeros(hours)
P_dis = np.zeros(hours)
P_grid = np.zeros(hours)
for t in range(hours):
net_power = P_load[t] - P_pv[t] # 净负荷
if net_power > 0: # 需要供电
# 优先使用储能放电
P_dis[t] = min(net_power / eta_dis, P_dis_max,
(SOC[t] - SOC_min) * E_max)
remaining = net_power - P_dis[t] * eta_dis
P_grid[t] = remaining # 剩余从电网购买
else: # 有多余光伏
surplus = -net_power
# 优先给储能充电
P_ch[t] = min(surplus * eta_ch, P_ch_max,
(SOC_max - SOC[t]) * E_max / eta_ch)
# 多余功率弃光（或上网）
# 更新SOC
SOC[t+1] = SOC[t] + (P_ch[t] * eta_ch - P_dis[t] / eta_dis) / E_max
# 计算成本
cost_grid = np.sum(P_grid * price)
revenue_pv_sell = np.sum(np.maximum(P_pv - P_load + P_dis/eta_dis - P_ch*eta_ch, 0) * price * 0.5)
print("微电网能量管理结果:")
print(f"从电网购电成本: ${cost_grid:.2f}")
print(f"光伏上网收益: ${revenue_pv_sell:.2f}")
print(f"净成本: ${cost_grid - revenue_pv_sell:.2f}")
# 可视化
fig, axes = plt.subplots(4, 1, figsize=(14, 12))
ax1 = axes[0]
ax1.plot(range(hours), P_pv, 'y-', label='PV', linewidth=2)
ax1.plot(range(hours), P_load, 'b-', label='Load', linewidth=2)
ax1.set_ylabel('Power (kW)')
ax1.set_title('Microgrid Power Profiles')
ax1.legend()
ax1.grid(True, alpha=0.3)
ax2 = axes[1]
ax2.bar(range(hours), P_ch, alpha=0.7, label='Charge', color='green')
ax2.bar(range(hours), -P_dis, alpha=0.7, label='Discharge', color='red')
ax2.set_ylabel('Power (kW)')
ax2.set_title('Battery Operation')
ax2.legend()
ax2.grid(True, alpha=0.3)
ax2.axhline(y=0, color='k', linewidth=0.5)
ax3 = axes[2]
ax3.plot(range(hours), SOC[:-1], 'g-', linewidth=2)
ax3.axhline(y=SOC_max, color='r', linestyle='--', label='SOC_max')
ax3.axhline(y=SOC_min, color='r', linestyle='--', label='SOC_min')
ax3.set_ylabel('SOC')
ax3.set_title('Battery State of Charge')
ax3.legend()
ax3.grid(True, alpha=0.3)
ax3.set_ylim(0, 1)
ax4 = axes[3]
ax4.bar(range(hours), P_grid, alpha=0.7, color='steelblue')
ax4.plot(range(hours), price * 100, 'r-', label='Price (x100)', linewidth=2)
ax4.set_xlabel('Hour')
ax4.set_ylabel('Power (kW) / Price')
ax4.set_title('Grid Import and Electricity Price')
ax4.legend()
ax4.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('microgrid_management.png', dpi=150)
plt.show()
10.5 案例5：电力系统暂态稳定性分析
分析简单电力系统的暂态稳定性，计算临界切除时间。
Python实现：
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
# 单机无穷大系统参数
H = 5.0 # 惯性常数 (s)
D = 0.5 # 阻尼系数
f_n = 50 # 额定频率 (Hz)
omega_s = 2 * np.pi * f_n # 同步角速度
# 系统等值电抗
X_d = 1.0 # 发电机暂态电抗
X_T = 0.1 # 变压器电抗
X_L = 0.5 # 线路电抗
X_total = X_d + X_T + X_L
# 运行参数
E = 1.1 # 发电机内电势
V_inf = 1.0 # 无穷大母线电压
P_m = 0.8 # 机械功率（标幺值）
# 计算功角特性
P_max_pre = E * V_inf / X_total # 故障前最大功率
P_max_fault = 0 # 三相短路时，P_max = 0
P_max_post = E * V_inf / (X_total * 1.5) # 故障后（单回线运行）
# 初始功角
delta_0 = np.arcsin(P_m / P_max_pre)
print(f"初始功角: {delta_0 * 180 / np.pi:.2f}°")
# 摇摆方程
def swing_equation(y, t, P_max, P_m, H, D, omega_s):
delta, omega = y
ddelta_dt = (omega - omega_s)
domega_dt = (omega_s / (2 * H)) * (P_m - P_max * np.sin(delta) - D * (omega - omega_s) / omega_s)
return [ddelta_dt, domega_dt]
# 仿真场景：三相短路，切除时间不同
t_clear_list = [0.1, 0.2, 0.3] # 切除时间 (s)
t_end = 5.0
t = np.linspace(0, t_end, 1000)
plt.figure(figsize=(14, 10))
for idx, t_clear in enumerate(t_clear_list):
# 分段仿真
# 阶段1: 故障前稳态 (已处于稳态)
# 阶段2: 故障期间
t_fault = np.linspace(0, t_clear, int(1000 * t_clear / t_end))
y0 = [delta_0, omega_s]
sol_fault = odeint(swing_equation, y0, t_fault,
args=(P_max_fault, P_m, H, D, omega_s))
# 阶段3: 故障后
t_post = np.linspace(t_clear, t_end, int(1000 * (t_end - t_clear) / t_end))
y0_post = sol_fault[-1]
sol_post = odeint(swing_equation, y0_post, t_post,
args=(P_max_post, P_m, H, D, omega_s))
# 合并结果
delta_full = np.concatenate([sol_fault[:, 0], sol_post[:, 0]])
omega_full = np.concatenate([sol_fault[:, 1], sol_post[:, 1]])
t_full = np.concatenate([t_fault, t_post])
# 绘制功角曲线
plt.subplot(2, 1, 1)
plt.plot(t_full, delta_full * 180 / np.pi,
label=f't_clear = {t_clear}s')
# 绘制转速曲线
plt.subplot(2, 1, 2)
plt.plot(t_full, (omega_full - omega_s) * 180 / np.pi,
label=f't_clear = {t_clear}s')
# 添加临界功角线
delta_critical = np.pi - np.arcsin(P_m / P_max_post)
plt.subplot(2, 1, 1)
plt.axhline(y=delta_critical * 180 / np.pi, color='r',
linestyle='--', label=f'Critical angle = {delta_critical * 180 / np.pi:.1f}°')
plt.xlabel('Time (s)')
plt.ylabel('Rotor Angle (°)')
plt.title('Transient Stability Analysis - Rotor Angle')
plt.legend()
plt.grid(True, alpha=0.3)
plt.subplot(2, 1, 2)
plt.xlabel('Time (s)')
plt.ylabel('Speed Deviation (°/s)')
plt.title('Rotor Speed Deviation')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('transient_stability.png', dpi=150)
plt.show()
10.6 案例6：电动汽车充电优化
优化大规模电动汽车的充电调度，降低电网峰荷和充电成本。
Python实现：
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
# 仿真参数
n_hours = 24
n_ev = 1000 # 电动汽车数量
# 到达时间和离开时间分布
np.random.seed(42)
arrival_hour = np.random.normal(18, 2, n_ev) % 24 # 主要傍晚到达
departure_hour = (arrival_hour + np.random.normal(8, 1, n_ev)) % 24 # 次日离开
# 电池参数
battery_capacity = 60 # kWh
initial_soc = np.random.uniform(0.2, 0.5, n_ev) # 到达时SOC
target_soc = 0.9 # 目标SOC
charging_power_max = 7 # kW（家用慢充）
charging_efficiency = 0.95
# 电价（分时电价）
price = np.array([0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.5, 0.5, 0.8, 0.8,
0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.5, 0.5, 0.8, 0.8,
0.8, 0.5, 0.5, 0.3])
# 基础负荷
base_load = np.array([40, 38, 36, 35, 36, 42, 55, 70, 80, 78,
75, 72, 70, 72, 75, 80, 85, 90, 88, 75,
65, 55, 48, 42]) # MW
# 简化优化：贪心算法
P_ch_opt = np.zeros((n_ev, n_hours))
def is_parked(i, t):
"""判断车辆i在时刻t是否停放"""
arr = arrival_hour[i]
dep = departure_hour[i]
if arr < dep:
return arr <= t < dep
else: # 跨午夜
return t >= arr or t < dep
for i in range(n_ev):
energy_needed = (target_soc - initial_soc[i]) * battery_capacity
parking_hours = [t for t in range(n_hours) if is_parked(i, t)]
if len(parking_hours) == 0:
continue
# 按电价排序充电时段
parking_prices = [(t, price[t]) for t in parking_hours]
parking_prices.sort(key=lambda x: x[1])
remaining_energy = energy_needed
for t, _ in parking_prices:
if remaining_energy <= 0:
break
max_charge = min(charging_power_max, remaining_energy / charging_efficiency)
P_ch_opt[i, t] = max_charge
remaining_energy -= max_charge * charging_efficiency
# 计算结果
total_charging = np.sum(P_ch_opt, axis=0) / 1000 # MW
print("电动汽车充电优化结果:")
print(f"总充电需求: {np.sum((target_soc - initial_soc) * battery_capacity):.2f} MWh")
print(f"总充电成本: ${np.sum(P_ch_opt * price / 1000):.2f}")
print(f"无EV时的峰荷: {np.max(base_load):.2f} MW")
print(f"有EV时的峰荷: {np.max(base_load + total_charging):.2f} MW")
# 可视化
fig, axes = plt.subplots(3, 1, figsize=(14, 12))
ax1 = axes[0]
ax1.plot(range(n_hours), base_load, 'b-', label='Base Load', linewidth=2)
ax1.plot(range(n_hours), base_load + total_charging, 'r-',
label='Base + EV Charging', linewidth=2)
ax1.set_ylabel('Power (MW)')
ax1.set_title('Load Profile with EV Charging')
ax1.legend()
ax1.grid(True, alpha=0.3)
ax2 = axes[1]
ax2.bar(range(n_hours), total_charging, alpha=0.7, color='green')
ax2.plot(range(n_hours), price * 20, 'r-', label='Price (x20)', linewidth=2)
ax2.set_ylabel('Power (MW) / Price')
ax2.set_title('EV Charging Power and Electricity Price')
ax2.legend()
ax2.grid(True, alpha=0.3)
ax3 = axes[2]
# 显示部分车辆的充电计划
for i in range(min(10, n_ev)):
ax3.plot(range(n_hours), P_ch_opt[i, :], alpha=0.5)
ax3.set_xlabel('Hour')
ax3.set_ylabel('Charging Power (kW)')
ax3.set_title('Individual EV Charging Schedules (Sample)')
ax3.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('ev_charging_optimization.png', dpi=150)
plt.show()
AtomGit 是由开放原子开源基金会联合 CSDN 等生态伙伴共同推出的新一代开源与人工智能协作平台。平台坚持“开放、中立、公益”的理念，把代码托管、模型共享、数据集托管、智能体开发体验和算力服务整合在一起，为开发者提供从开发、训练到部署的一站式体验。
更多推荐
3
0- 0
已为社区贡献41条内容
所有评论(0)