import os
os.environ['KMP_DUPLICATE_LIB_OK'] = 'TRUE' # 添加到代码开头
import cvxpy as cp
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import warnings
from datetime import datetime, timedelta
import os
# ==============================================================================
# 网络配置文件解析函数,按标签段读取自定义 .e 网络配置文件,把支路参数、光伏/储能/台区节点号以及平衡节点号解析成结构化字典,供主程序直接调用。
# ==============================================================================
def load_network_config(config_filepath):
config_data = {
'branch': [],
'pv_nodes': [],
'ess_nodes': [],
'district_nodes': [],
'slack_bus_node': None # 初始化 slack_bus_node
}
current_section = None
with open(config_filepath, 'r', encoding='utf-8') as f:
for line in f:
line = line.strip()
if not line or line.startswith('//'):
continue
if line.startswith('<') and not line.startswith('</'):
if line == '<distributed_network>':
current_section = 'branch'
elif line == '<pv_point>':
current_section = 'pv'
elif line == '<ess_point>':
current_section = 'ess'
elif line == '<region_point>':
current_section = 'district'
elif line == '<slack_bus>':
current_section = 'slack' # 新增识别
else:
current_section = None
continue
if line.startswith('</'):
current_section = None
continue
if line.startswith('#'):
parts = line.lstrip('#').strip().split()
if not parts: continue
if current_section == 'branch' and len(parts) >= 5:
config_data['branch'].append([int(parts[1]), int(parts[2]), float(parts[3]), float(parts[4])])
elif current_section == 'pv':
config_data['pv_nodes'].append(int(parts[0]))
elif current_section == 'ess':
config_data['ess_nodes'].append(int(parts[0]))
elif current_section == 'district':
config_data['district_nodes'].append(int(parts[0]))
elif current_section == 'slack': # 新增解析
config_data['slack_bus_node'] = int(parts[0])
return config_data
# ==============================================================================
# 一、数据输入(一定修改为正确的文件读取路径)
# ==============================================================================
# --- 使用相对路径,请确保数据文件与.py文件在同一目录下,或提供正确绝对路径 ---
import argparse
parser = argparse.ArgumentParser(description='电压优化算法')
parser.add_argument('--data', required=True, help='数据文件路径(.e)')
parser.add_argument('--config', required=True, help='配置文件路径(.e)')
parser.add_argument('--output', required=True, help='输出目录路径')
args = parser.parse_args()
data_file_path = args.data
config_file_path = args.config
hdr_pv, rows_pv = None, []
hdr_load, rows_load = None, []
try:
with open(data_file_path, 'r', encoding='utf-8', errors='ignore') as f:
current = None
for line in f:
line = line.rstrip('\n')
if line.strip() == '<point_output_pv>':
current = 'pv'
elif line.strip() == '<load_point>':
current = 'load' #通过current状态变量标记当前读取的数据类型(PV或Load)到新的标签段(如<other_section>)时重置状态
elif line.startswith('<') and line.endswith('>'):
current = None
else:
if current == 'pv':
if line.startswith('@'):
hdr_pv = [h for h in line.lstrip('@').strip().split() if h]
elif line.startswith('#'):
vals = [v for v in line.lstrip('#').strip().split() if v]
if len(vals) == len(hdr_pv): rows_pv.append(vals)
elif current == 'load':
if line.startswith('@'):
hdr_load = [h for h in line.lstrip('@').strip().split() if h]
elif line.startswith('#'):
vals = [v for v in line.lstrip('#').strip().split() if v]
if len(vals) == len(hdr_load): rows_load.append(vals)
except FileNotFoundError:
print(f"错误:数据文件 '{data_file_path}' 未找到。请检查路径。")
exit()
df_pv = pd.DataFrame(rows_pv, columns=hdr_pv)
df_pv['Time'] = pd.to_datetime(df_pv['Time'])
df_pv['Pv'] = pd.to_numeric(df_pv['Pv'], errors='coerce')
df_pv['Power_pv_pre'] = pd.to_numeric(df_pv['Power_pv_pre'], errors='coerce')
df_load = pd.DataFrame(rows_load, columns=hdr_load)
df_load['Time'] = pd.to_datetime(df_load['Time']) # 确保负荷时间也被解析
df_load['Load_Bus_Pre'] = pd.to_numeric(df_load['Load_Bus_Pre'], errors='coerce')
#将时间列转换为Pandas的datetime类型,errors='coerce'自动将无效数值转为NaN(如空字符串或非数字字符)
# ==============================================================================
# 二、主程序
# ==============================================================================
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
warnings.filterwarnings("ignore", category=UserWarning)
#设置Matplotlib中文字体显示(避免中文乱码),关闭负数显示警告,忽略UserWarning类警告
#SimHei是Windows系统黑体字库,unicode_minus解决负号显示为方块的问题,警告过滤可避免后续优化计算时的冗余提示
# ------------------- 1. 参数设置 -------------------
T = 96
print(f"正在从 '{config_file_path}' 加载网络配置...")
try:
network_config = load_network_config(config_file_path)
except FileNotFoundError:
print(f"错误:配置文件 '{config_file_path}' 未找到。")
exit()
# --- 1.2 动态定义支路、节点、设备 ---
raw_branch_data = network_config['branch']
branch_list = [[row[0], row[1], row[2] + 1j * row[3]] for row in raw_branch_data]
branch = np.array(branch_list)#将支路电阻(R)和电抗(X)合并为复数阻抗
# 从拓扑中获取所有节点,确定节点数nb和支路数nl,np.real(branch[:, 0:2])提取所有首末节点编号,np.unique去重后得到全网节点列表
all_node_ids = sorted(np.unique(np.real(branch[:, 0:2]).astype(int)))
nb = len(all_node_ids)
nl = len(raw_branch_data)
# 读取并验证平衡节点,必须有且仅有一个平衡节点(参考节点)
slack_bus_id = network_config.get('slack_bus_node')
if slack_bus_id is None:
raise ValueError("错误:配置文件中未通过 <slack_bus> 指定平衡节点!")
# 识别负荷节点(所有非平衡节点)
load_node_ids = sorted([nid for nid in all_node_ids if nid != slack_bus_id])#非平衡节点即负荷节点
pv_nodes = network_config['pv_nodes']#光伏数量
ess_nodes = network_config['ess_nodes']#储能数量
district_nodes = network_config['district_nodes']#台区数量
npv = len(pv_nodes)
ness = len(ess_nodes)
ndistrict = len(district_nodes)
print("\n网络配置加载成功:")
print(f" - 系统节点: {all_node_ids} (共 {nb} 个)")
print(f" - 平衡节点: {slack_bus_id}")
print(f" - 负荷节点: {load_node_ids} (共 {len(load_node_ids)} 个)")
print(f" - 支路数 (nl): {nl}")
print(f" - 主网光伏数量 (npv): {npv}, 接入节点: {pv_nodes}")
print(f" - 主网储能数量 (ness): {ness}, 接入节点: {ess_nodes}")
print(f" - 台区数量 (ndistrict): {ndistrict}, 接入节点: {district_nodes}\n")
# --- 1.3 精确分配负荷和光伏数据 ---
df_load['Bus_i'] = pd.to_numeric(df_load['Bus_i'])
load_data_dict = {
bus_id: group['Load_Bus_Pre'].iloc[0:T].values / 1000#groupby('Bus_i')按节点分组,iloc[0:T]取前T个时间点(T=96对应24小时15分钟间隔),/1000单位转换
for bus_id, group in df_load.groupby('Bus_i')
}
for bus_id, data in load_data_dict.items():
if len(data) < T:
raise ValueError(f"数据文件中负荷组 {bus_id} 的数据点不足 {T} 个。")
# 创建 pload 矩阵并按顺序分配
pload = np.zeros((nb, T))
node_to_idx = {node_id: i for i, node_id in enumerate(all_node_ids)}
print("\n正在为负荷节点精确分配负荷数据...")
load_groups_in_file = sorted(load_data_dict.keys())
if len(load_groups_in_file) < len(load_node_ids):
raise ValueError(
f"数据文件中有 {len(load_groups_in_file)} 组负荷数据,不足以分配给 {len(load_node_ids)} 个负荷节点。")
#pload矩阵维度:(节点数, 时间步长),映射逻辑:将数据文件中的负荷ID映射到网络拓扑节点
for i, target_node_id in enumerate(load_node_ids):
source_load_id = load_groups_in_file[i]
target_idx = node_to_idx[target_node_id]
pload[target_idx, :] = load_data_dict[source_load_id]
print(f" - 负荷组 {source_load_id} (来自文件) -> 节点 {target_node_id}")
# 动态为每个光伏站分配功率曲线
print(f"\n正在为 {npv} 个光伏站分配功率曲线...")
p_pv = np.zeros((npv, T))
for i in range(npv):
pv_station_id = (i % 2) + 1 # 示例文件中只有1和2,这里做循环使用
station_df = df_pv[df_pv['Pv'] == pv_station_id]
if len(station_df) < T:
raise ValueError(f"数据文件 '{data_file_path}' 中光伏站 {pv_station_id} 的数据点不足 {T} 个。")
power_curve = station_df['Power_pv_pre'].iloc[0:T].values / 1000
p_pv[i, :] = np.clip(power_curve, a_min=0, a_max=None)#np.clip确保无负功率(夜间可能出现的仪器误差)
print("光伏功率曲线分配完成。")
# --- 1.4 基于有功功率大小分档设定无功负荷 ---
print("\n正在根据功率因数计算无功负荷...")
load_thresholds = {'heavy': 0.200, 'light': 0.050}
power_factors = {'heavy': 0.85, 'medium': 0.92, 'light': 0.98}#根据有功负荷大小分三档设置功率因数(PF)
tan_phi_map = {level: np.tan(np.arccos(pf)) for level, pf in power_factors.items()}#通过tan(arccos(PF))计算无功/有功比
conditions = [
pload > load_thresholds['heavy'],
(pload > load_thresholds['light']) & (pload <= load_thresholds['heavy']),
pload <= load_thresholds['light']
]#定义无功计算方式
choices = [pload * tan_phi_map['heavy'], pload * tan_phi_map['medium'], pload * tan_phi_map['light']]
qload = np.select(conditions, choices, default=0.0)
print("无功负荷计算完成。")
# --- 1.5 在优化前进行数据时间对齐 ---
print("\n--- 正在进行数据时间对齐 ---")
# 1. 计算光伏数据的偏移量
start_time_pv = df_pv['Time'].iloc[0]
shift_pv = (start_time_pv.hour * 60 + start_time_pv.minute) // 15#将数据起始时间转换为15分钟间隔的步长数(如08:30→8×4 + 2=34个步长)
print(f"光伏数据起始时间: {start_time_pv.strftime('%H:%M')},对应时间步偏移量: {shift_pv}")
# 2. 计算负荷数据的偏移量
start_time_load = df_load['Time'].iloc[0]
shift_load = (start_time_load.hour * 60 + start_time_load.minute) // 15
print(f"负荷数据起始时间: {start_time_load.strftime('%H:%M')},对应时间步偏移量: {shift_load}")
# 3. 使用np.roll进行循环移位,将数据对齐到00:00开始,确保所有数据从00:00开始对齐
p_pv = np.roll(p_pv, shift=shift_pv, axis=1)
pload = np.roll(pload, shift=shift_load, axis=1)
qload = np.roll(qload, shift=shift_load, axis=1)
print("数据已对齐至00:00-23:45时间序列。优化将基于对齐后的数据进行。")
print("---------------------------------\n")
# --- 1.6 设置其他参数 ---
branch[:, 2] = branch[:, 2] / (10 ** 2)#阻抗标幺化
R = np.real(branch[:, 2])#支路电阻
X = np.imag(branch[:, 2])#支路电抗
c_ess = [1] * ness
p_essmax = [0.3] * ness
s_essmax = [0.3 * 1.1] * ness
eta_ch = [0.9] * ness
eta_dis = [0.9] * ness
p_district_max = [1.0] * ndistrict
q_district_max = [0.5] * ndistrict
# ------------------- 2. 拓扑结构 -------------------
upstream = np.zeros((nb, nl))#创建 nb×nl 全零矩阵,行对应节点,列对应支路;后续标记“支路 i 的首端是哪个节点”。
dnstream = np.zeros((nb, nl))
for i in range(nl):#逐条支路进行处理
from_node_id = int(np.real(branch[i, 0]))#取出第 i 条支路的首端节点编号。
to_node_id = int(np.real(branch[i, 1]))#取出第 i 条支路的末端节点编号。
from_node_idx = node_to_idx[from_node_id]#把首端节点编号映射到 0‥nb-1 的行索引。
to_node_idx = node_to_idx[to_node_id]#把末端节点编号映射到 0‥nb-1 的行索引。
upstream[from_node_idx, i] = 1#在 upstream 矩阵的 (首端节点行, 支路列) 置 1,表示“节点是该支路的起点”。
dnstream[to_node_idx, i] = 1#在 dnstream 矩阵的 (末端节点行, 支路列) 置 1,表示“节点是该支路的终点”。
Vmax = 1.06 * 1.06 * np.ones((nb, T))#给所有节点、所有时段设置电压幅值平方的上限 (1.06 p.u.)²。
Vmin = 0.94 * 0.94 * np.ones((nb, T))#给所有节点、所有时段设置电压幅值平方的下限 (0.94 p.u.)²。
Pgmax = np.zeros((nb, T))#建立四个 nb×T 零矩阵,用于后续设定各节点有功/无功出力上下界。
Qgmax = np.zeros((nb, T))
Pgmin = np.zeros((nb, T))
Qgmin = np.zeros((nb, T))
slack_bus_idx = node_to_idx[slack_bus_id]#找到平衡节点在矩阵中的行索引。
Pgmax[slack_bus_idx, :] = 10#仅对平衡节点设定大范围的出力上下界(±10 MW,±10 MVAr),其余节点保持 0,后续再细调。
Pgmin[slack_bus_idx, :] = -10 # (MW)
Qgmax[slack_bus_idx, :] = 10
Qgmin[slack_bus_idx, :] = -10 # (MVar)
# ------------------- 3. 定义决策变量 -------------------
V = cp.Variable((nb, T))#全网各节点各时段电压幅值平方(标幺)。
I = cp.Variable((nl, T))#各支路各时段电流幅值平方(标幺)。
P = cp.Variable((nl, T))#各支路各时段有功功率。
Q = cp.Variable((nl, T))#各支路各时段无功功率。
Pg = cp.Variable((nb, T))#各节点各时段注入有功功率(含光伏、储能、平衡机等)。
Qg = cp.Variable((nb, T))#各节点各时段注入无功功率。
q_pv = cp.Variable((npv, T))#每个光伏站各时段的无功出力。
p_ch = cp.Variable((ness, T), nonneg=True)#每套储能各时段的充电有功功率(≥0)。
p_dis = cp.Variable((ness, T), nonneg=True)
u_ch = cp.Variable((ness, T), boolean=True)
u_dis = cp.Variable((ness, T), boolean=True)
e_ess = cp.Variable((ness, T))#每套储能各时段的剩余电量(能量)。
q_ess = cp.Variable((ness, T))#每套储能各时段的无功出力。
p_pv_utilized = cp.Variable((nb, T), nonneg=True)#各节点实际被消纳的光伏有功功率(≥0)。
P_district_exchange = cp.Variable((ndistrict, T))#各台区与主网交换的有功功率。
Q_district_exchange = cp.Variable((ndistrict, T))#各台区与主网交换的无功功率。
# ------------------- 4. 设置约束 -------------------
constraints = []
pv_nodes_idx = [node_to_idx[node] for node in pv_nodes]
ess_nodes_idx = [node_to_idx[node] for node in ess_nodes]
district_nodes_idx = [node_to_idx[node] for node in district_nodes]
# 4.1 主网光伏约束
#把光伏无功 q_pv 映射到对应节点,限制无功不超过有功 ±30%,并把实际被利用的光伏有功 p_pv_utilized 约束在 0 与可用光伏之间。
P_pv_matrix = np.zeros((nb, T))
for i, node_idx in enumerate(pv_nodes_idx): P_pv_matrix[node_idx, :] = p_pv[i, :]
Q_pv = cp.Variable((nb, T))
for i, node_idx in enumerate(pv_nodes_idx): constraints += [Q_pv[node_idx, :] == q_pv[i, :]]
non_pv_nodes_idx = [i for i in range(nb) if i not in pv_nodes_idx]
for node_idx in non_pv_nodes_idx: constraints += [Q_pv[node_idx, :] == 0]
constraints += [-0.3 * P_pv_matrix[pv_nodes_idx, :] <= q_pv, q_pv <= 0.3 * P_pv_matrix[pv_nodes_idx, :]]
constraints += [p_pv_utilized[pv_nodes_idx, :] <= P_pv_matrix[pv_nodes_idx, :]]
constraints += [p_pv_utilized[non_pv_nodes_idx, :] == 0, p_pv_utilized >= 0]
# 4.2 主网储能约束
#通过 0-1 变量 u_ch / u_dis 保证“充放电互锁”,功率、容量、SOC 均限制在设定范围,且保持日循环(首尾 SOC 相等)。
P_dis_matrix = cp.Constant(np.zeros((nb, T)))
P_ch_matrix = cp.Constant(np.zeros((nb, T)))
Q_ess_matrix = cp.Constant(np.zeros((nb, T)))
for i, node_idx in enumerate(ess_nodes_idx):
# 使用辅助矩阵来构建,避免CVXPY的+=性能问题
temp_p_dis = np.zeros((nb, T))
temp_p_dis[node_idx, :] = p_dis[i, :].value if p_dis.value is not None else 0
P_dis_matrix += cp.vstack(
[cp.Constant(np.zeros((node_idx, T))), p_dis[i:i + 1, :], cp.Constant(np.zeros((nb - node_idx - 1, T)))])
P_ch_matrix += cp.vstack(
[cp.Constant(np.zeros((node_idx, T))), p_ch[i:i + 1, :], cp.Constant(np.zeros((nb - node_idx - 1, T)))])
Q_ess_matrix += cp.vstack(
[cp.Constant(np.zeros((node_idx, T))), q_ess[i:i + 1, :], cp.Constant(np.zeros((nb - node_idx - 1, T)))])
constraints += [u_ch + u_dis <= 1]
for i in range(ness):
constraints += [0 <= p_dis[i, :], p_dis[i, :] <= u_dis[i, :] * p_essmax[i]]
constraints += [0 <= p_ch[i, :], p_ch[i, :] <= u_ch[i, :] * p_essmax[i]]
for t in range(T):
constraints += [cp.norm(cp.vstack([p_dis[i, t] - p_ch[i, t], q_ess[i, t]])) <= s_essmax[i]]
for i in range(ness):
constraints += [e_ess[i, 0] == c_ess[i] / 2]
for t in range(T - 1):
constraints += [
e_ess[i, t + 1] == e_ess[i, t] + eta_ch[i] * p_ch[i, t] * 0.25 - (p_dis[i, t] / eta_dis[i]) * 0.25]
constraints += [e_ess[i, 0] == e_ess[i, T - 1]]
constraints += [c_ess[i] * 0.2 <= e_ess[i, :], e_ess[i, :] <= c_ess[i] * 0.9]
# 4.3 台区约束
#将台区与主网的交换功率 P_district_exchange / Q_district_exchange 限制在给定上下界,并映射到对应节点。
P_district_grid = cp.Constant(np.zeros((nb, T)))
Q_district_grid = cp.Constant(np.zeros((nb, T)))
for d, node_idx in enumerate(district_nodes_idx):
constraints += [-p_district_max[d] <= P_district_exchange[d, :], P_district_exchange[d, :] <= p_district_max[d]]
constraints += [-q_district_max[d] <= Q_district_exchange[d, :], Q_district_exchange[d, :] <= q_district_max[d]]
P_district_grid += cp.vstack([cp.Constant(np.zeros((node_idx, T))), P_district_exchange[d:d + 1, :],
cp.Constant(np.zeros((nb - node_idx - 1, T)))])
Q_district_grid += cp.vstack([cp.Constant(np.zeros((node_idx, T))), Q_district_exchange[d:d + 1, :],
cp.Constant(np.zeros((nb - node_idx - 1, T)))])
# 4.4 主网潮流约束 (修正了潮流方程中的储能和台区项)
#用线性化 DistFlow 描述支路功率平衡、节点功率平衡、电压降落;对每条支路每时段施加二阶锥松弛,使模型可解且保持凸性。
Pin = -upstream @ P + upstream @ cp.multiply(I, R[:, None] * np.ones((1, T))) + dnstream @ P
Qin = -upstream @ Q + upstream @ cp.multiply(I, X[:, None] * np.ones((1, T))) + dnstream @ Q
# 节点功率平衡方程
constraints += [Pin + pload - Pg - p_pv_utilized - P_dis_matrix + P_ch_matrix - P_district_grid == 0]
constraints += [Qin + qload - Qg - Q_pv - Q_ess_matrix - Q_district_grid == 0]
# 潮流方程
for i in range(nl):
from_node_idx = node_to_idx[int(np.real(branch[i, 0]))]
to_node_idx = node_to_idx[int(np.real(branch[i, 1]))]
constraints += [V[to_node_idx, :] == V[from_node_idx, :] - 2 * R[i] * P[i, :] - 2 * X[i] * Q[i, :] + (
R[i] ** 2 + X[i] ** 2) * I[i, :]]
# 二阶锥约束
for i in range(nl):
from_node_idx = node_to_idx[int(np.real(branch[i, 0]))]
for t in range(T):
constraints += [cp.norm(cp.vstack([2 * P[i, t], 2 * Q[i, t], I[i, t] - V[from_node_idx, t]])) <= I[i, t] + V[
from_node_idx, t]]
# 4.5 通用约束
#统一给所有节点加上电压、功率、电流的硬边界,并把平衡节点电压固定为 1 p.u.。
constraints += [Vmin <= V, V <= Vmax]
constraints += [Pgmin <= Pg, Pg <= Pgmax]
constraints += [Qgmin <= Qg, Qg <= Qgmax]
constraints += [0 <= I, I <= 1000] # 放宽电流约束以防无解
constraints += [V[slack_bus_idx, :] == 1.0]
# ------------------- 5. 设置目标函数 -------------------
optimization_target = "voltage"
objective = cp.Minimize(cp.sum(cp.abs(V - 1)))#电压优化目标:最小化所有节点电压幅值平方与1 p.u.的绝对偏差总和
# ------------------- 6. 求解优化问题 -------------------
#把目标函数和所有约束塞进 cp.Problem,调用 Gurobi 求解器,设置 3 % 的 MIPGap,开始求解。
problem = cp.Problem(objective, constraints)
print("\n正在调用求解器求解...")
#problem.solve(solver=cp.GUROBI, verbose=True, MIPGap=0.03)
# 使用 SCIP 求解器,设置 Gap = 0.03
problem.solve(solver=cp.SCIP, verbose=True, scip_params={"limits/gap": 0.03})
# ------------------- 7. 结果分析 -------------------
print("\n--- 结果分析 ---")
if problem.status in [cp.OPTIMAL, cp.OPTIMAL_INACCURATE]:#判断求解器返回的状态是不是“已最优”或“近似最优”。
print("成功求解!")
if optimization_target == "voltage":
voltage_deviation = cp.sum(cp.abs(V - 1)).value#把目标函数值取出来,即所有节点电压与 1 p.u. 的绝对偏差总和。
print("所有节点的电压偏差总和:", voltage_deviation)
elif optimization_target == "pv_utilization":
total_pv_available = np.sum(P_pv_matrix)#计算主网光伏理论可发总量。
if total_pv_available > 1e-6:#确保分母不为 0。
pv_utilization_rate_main = np.sum(p_pv_utilized.value[pv_nodes_idx, :]) / total_pv_available#用“实际被利用光伏总量 / 理论可发总量”得出主网光伏消纳率。
print(f"主网光伏消纳率: {pv_utilization_rate_main:.3%}")
else:
print("主网光伏可用总量为0,无法计算消纳率。")
else:
print("求解出错:", problem.status)
# ------------------- 8. 生成 .e 文件 -------------------
if problem.status in [cp.OPTIMAL, cp.OPTIMAL_INACCURATE]:#仅在 OPTIMAL / OPTIMAL_INACCURATE 状态下才继续;
output_dir = args.output#定义输出文件夹的绝对路径。
os.makedirs(output_dir, exist_ok=True)#如果该路径不存在就自动创建,若已存在则静默跳过。
start_time = datetime.now().replace(hour=0, minute=0, second=0, microsecond=0)#把当前日期截成“今天 00:00:00”,作为 24 h 时间轴的起点。
time_steps = [start_time + timedelta(minutes=15 * t) for t in range(T)]#以 15 min 为步长,从 00:00 开始生成 T=96 个 datetime 对象,供后续逐时刻写文件使用。
# 台区 .e 文件
region_filename = os.path.join(output_dir, f"region_action_plan_{start_time.strftime('%Y%m%d%H%M')}.e")#把目录、前缀、时间戳拼成完整的台区计划文件路径。
with open(region_filename, 'w', encoding='utf-8') as f:#以 UTF-8 编码、写模式打开文件,准备写入。
f.write("<region_action_plan>\n")#写入段起始标签 <region_action_plan>。
f.write("//\t时间\t台区接入节点号\t台区有功功率(kW)\t台区无功功率(kVar)\n")#写入注释行,说明后面各列含义。
f.write("@\ttime\tregion_node_id\tregion_active_power\tregion_reactive_power\n")#写入列名行(以 @ 开头,后续解析用)。
for i, node_id in enumerate(district_nodes):#遍历每一个台区节点编号。
for t in range(T):#对该节点遍历 0–95 共 96 个时段。
time_str = time_steps[t].strftime('%Y-%m-%dT%H:%M:%S')#把当前时段转成 “YYYY-MM-DDTHH:MM:SS” 字符串格式。
p_exchange = P_district_exchange.value[i, t] * 1000#把求解得到的台区有功交换功率从 MW → kW。
q_exchange = Q_district_exchange.value[i, t] * 1000#把无功交换功率从 MVAr → kVAr。
f.write(f"#\t{time_str}\t{node_id}\t{p_exchange:.4f}\t{q_exchange:.4f}\n")#以 # 开头写一行数据:时间、节点号、有功(kW)、无功(kVAr),均保留 4 位小数。
f.write("</region_action_plan>\n")#写入段结束标签。
print(f".e 文件已生成: {region_filename}")#提示用户文件已成功写出,并显示完整路径。
# 光伏出力计划 .e 文件
pv_filename = os.path.join(output_dir, f"pv_action_plan_{start_time.strftime('%Y%m%d%H%M')}.e")
with open(pv_filename, 'w', encoding='utf-8') as f:
f.write("<pv_action_plan>\n")
f.write("//\t时间\t光伏接入节点号\t光伏输出有功(kW)\t光伏输出无功(kVar)\n")
f.write("@\ttime\tpv_node_id\tpv_active_power\tpv_reactive_power\n")
for i, node_id in enumerate(pv_nodes):
node_idx = pv_nodes_idx[i]
for t in range(T):
time_str = time_steps[t].strftime('%Y-%m-%dT%H:%M:%S')
p_pv_util = p_pv_utilized.value[node_idx, t] * 1000
q_pv_val = q_pv.value[i, t] * 1000
f.write(f"#\t{time_str}\t{node_id}\t{p_pv_util:.4f}\t{q_pv_val:.4f}\n")
f.write("</pv_action_plan>\n")
print(f".e 文件已生成: {pv_filename}")
# 储能计划动作 .e 文件
storage_filename = os.path.join(output_dir, f"storage_action_plan_{start_time.strftime('%Y%m%d%H%M')}.e")
with open(storage_filename, 'w', encoding='utf-8') as f:
f.write("<storage_action_plan>\n")
f.write("//\t时间\t储能接入节点号\t储能输出有功(kW)\t储能输出无功(kVar)\n")
f.write("@\ttime\tstorage_node_id\tstorage_active_power\tstorage_reactive_power\n")
for i, node_id in enumerate(ess_nodes):
for t in range(T):
time_str = time_steps[t].strftime('%Y-%m-%dT%H:%M:%S')
p_ess = (p_dis.value[i, t] - p_ch.value[i, t]) * 1000
q_ess_val = q_ess.value[i, t] * 1000
f.write(f"#\t{time_str}\t{node_id}\t{p_ess:.4f}\t{q_ess_val:.4f}\n")
f.write("</storage_action_plan>\n")
print(f".e 文件已生成: {storage_filename}")
else:
print("因求解失败,无法生成 .e 文件。")
# ------------------- 9. 可视化图形 -------------------
if problem.status in [cp.OPTIMAL, cp.OPTIMAL_INACCURATE]:
hours = np.linspace(0, 24, T, endpoint=False)
major_ticks = np.arange(0, 25, 4)
# 图表1&2:台区功率流动
# 移除台区数量条件限制,确保图表生成
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 10), sharex=True)
for i in range(ndistrict):
ax1.plot(hours, P_district_exchange.value[i, :], marker='o', linestyle='-',
label=f'台区 (节点 {district_nodes[i]})', markersize=4)
ax2.plot(hours, Q_district_exchange.value[i, :], marker='s', linestyle='--',
label=f'台区 (节点 {district_nodes[i]})', markersize=4)
ax1.set_title('台区24小时功率流动')
ax1.set_ylabel('有功功率 (MW)')
ax1.grid(True)
ax1.legend()
ax2.set_xlabel('时间 (小时)')
ax2.set_ylabel('无功功率 (MVar)')
ax2.grid(True)
ax2.legend()
plt.xticks(major_ticks)
plt.tight_layout()
plt.show()
这是其中的一个优化算法,在调用优化算法进行优化后,如何将优化结果(非图片格式)传送到GUI界面那里,然后在GUI完成数据分析再画图呢?
最新发布