1.每一次计算都重新画画。
This commit is contained in:
dmy
parent
13e25832ed
commit
7a5bb05f58
|
|
@ -0,0 +1,13 @@
|
|||
*.dxf
|
||||
build
|
||||
__pycache__
|
||||
CSharp
|
||||
.idea
|
||||
dist
|
||||
*.spec
|
||||
*.dwg
|
||||
历史
|
||||
.venv
|
||||
*.toml
|
||||
launch.json
|
||||
settings.json
|
||||
11
animation.py
11
animation.py
|
|
@ -25,18 +25,20 @@ class Animation:
|
|||
if not cls._disable:
|
||||
# print("desc")
|
||||
return func(cls, *args, **kwargs)
|
||||
return not_run(cls,*args, **kwargs)
|
||||
return not_run(cls, *args, **kwargs)
|
||||
|
||||
return wrapTheFunction
|
||||
|
||||
def disable(self, _disable):
|
||||
self._disable = _disable
|
||||
|
||||
@switch_decorator
|
||||
def init_fig(self):
|
||||
ax = self._ax
|
||||
ax.set_aspect(1)
|
||||
ax.set_xlim([-500, 500])
|
||||
ax.set_ylim([-500, 500])
|
||||
|
||||
@switch_decorator
|
||||
def show(self):
|
||||
self._fig.show()
|
||||
|
|
@ -47,10 +49,12 @@ class Animation:
|
|||
x = np.linspace(0, 300)
|
||||
y = line_func(x)
|
||||
ax.plot(x, y)
|
||||
|
||||
@switch_decorator
|
||||
def add_rs(self, rs, rs_x, rs_y):
|
||||
ax = self._ax
|
||||
ax.add_artist(plt.Circle((rs_x, rs_y), rs, fill=False))
|
||||
|
||||
@switch_decorator
|
||||
def add_rc(self, rc, rc_x, rc_y):
|
||||
ax = self._ax
|
||||
|
|
@ -71,15 +75,14 @@ class Animation:
|
|||
ax.plot([rc_x, intersection_x1], [rc_y, intersection_y1], color="red")
|
||||
ax.plot([rc_x, intersection_x2], [rc_y, intersection_y2], color="red")
|
||||
pass
|
||||
|
||||
@switch_decorator
|
||||
def clear(self):
|
||||
ax = self._ax
|
||||
# fig = self._fig
|
||||
ax.cla()
|
||||
|
||||
@switch_decorator
|
||||
def pause(self):
|
||||
# fig=self._fig
|
||||
# print('tick')
|
||||
ax = self._ax
|
||||
self._ticks += 1
|
||||
ticks = self._ticks
|
||||
|
|
|
|||
198
main.py
198
main.py
|
|
@ -1,10 +1,12 @@
|
|||
import math
|
||||
import os.path
|
||||
import sys
|
||||
import time
|
||||
import tomli
|
||||
from loguru import logger
|
||||
from core import *
|
||||
import timeit
|
||||
from animation import Animation
|
||||
|
||||
|
||||
# 打印参数
|
||||
|
|
@ -59,7 +61,7 @@ def read_parameter(toml_file_path):
|
|||
|
||||
def egm():
|
||||
if len(sys.argv) < 2:
|
||||
toml_file_path = r"article.toml"
|
||||
toml_file_path = r"内自500kV-ZCK上相.toml"
|
||||
else:
|
||||
toml_file_path = sys.argv[1]
|
||||
if not os.path.exists(toml_file_path):
|
||||
|
|
@ -70,7 +72,7 @@ def egm():
|
|||
#########################################################
|
||||
# 以上是需要设置的参数
|
||||
parameter_display(para)
|
||||
h_whole = para.h_arm[0] # 塔全高
|
||||
h_whole = para.h_arm[0] # 挂点高
|
||||
string_g_len = para.string_g_len
|
||||
string_c_len = para.string_c_len
|
||||
h_g_sag = para.h_g_sag
|
||||
|
|
@ -82,7 +84,7 @@ def egm():
|
|||
]
|
||||
if len(h_arm) > 1:
|
||||
for hoo in h_arm[1:]:
|
||||
gc_y.append(hoo - string_c_len - h_c_sag * 2 / 3)
|
||||
gc_y.append(hoo - string_c_len - h_c_sag * 2 / 3) # 导线平均高
|
||||
if len(gc_y) > 2: # 双回路
|
||||
phase_n = 3 # 边相导线数量
|
||||
else:
|
||||
|
|
@ -92,12 +94,15 @@ def egm():
|
|||
ng = func_ng(td)
|
||||
avr_n_sf = 0 # 考虑电压的影响计算的跳闸率
|
||||
ground_angels = para.ground_angels
|
||||
# 初始化动画
|
||||
animate = Animation()
|
||||
animate.disable(False)
|
||||
# animate.show()
|
||||
for ground_angel in ground_angels:
|
||||
logger.info(f"地面倾角{ground_angel/math.pi*180:.3f}°")
|
||||
rg_type = None
|
||||
rg_x = None
|
||||
rg_y = None
|
||||
cad = Draw()
|
||||
voltage_n = para.voltage_n
|
||||
n_sf_phases = np.zeros((phase_n, voltage_n)) # 存储每一相的跳闸率
|
||||
if np.any(np.array(gc_y) < 0):
|
||||
|
|
@ -113,13 +118,13 @@ def egm():
|
|||
rg_type = "g"
|
||||
if phase_n > 1: # 多回路
|
||||
if phase_conductor_foo < 2:
|
||||
rg_type = "c" # 捕捉弧有下面一相导线的击距代替
|
||||
rg_type = "c" # 捕捉弧由下面一相导线的击距代替
|
||||
rg_x = gc_x[phase_conductor_foo + 2]
|
||||
rg_y = gc_y[phase_conductor_foo + 2]
|
||||
else:
|
||||
rg_type = "g"
|
||||
# TODO 保护角公式可能有问题，后面改
|
||||
shield_angle = (
|
||||
shield_angle_at_avg_height = (
|
||||
math.atan(
|
||||
(rc_x - rs_x)
|
||||
/ (
|
||||
|
|
@ -129,78 +134,113 @@ def egm():
|
|||
)
|
||||
* 180
|
||||
/ math.pi
|
||||
) # 保护角
|
||||
logger.info(f"保护角{shield_angle:.3f}°")
|
||||
) # 挂点处保护角
|
||||
logger.info(f"挂点处保护角{shield_angle_at_avg_height:.3f}°")
|
||||
logger.debug(f"最低相防护标识{rg_type}")
|
||||
rated_voltage = para.rated_voltage
|
||||
for u_bar in range(voltage_n): # 计算不同工作电压下的跳闸率
|
||||
u_ph = (
|
||||
math.sqrt(2)
|
||||
* 750
|
||||
* math.cos(2 * math.pi / voltage_n * u_bar)
|
||||
/ 1.732
|
||||
) # 运行相电压
|
||||
# TODO 需要区分交、直流
|
||||
# u_ph = (
|
||||
# math.sqrt(2)
|
||||
# * rated_voltage
|
||||
# * math.cos(2 * math.pi / voltage_n * u_bar)
|
||||
# / 1.732
|
||||
# ) # 运行相电压
|
||||
u_ph = rated_voltage / 1.732
|
||||
logger.info(f"计算第{phase_conductor_foo + 1}相，电压为{u_ph:.2f}kV")
|
||||
# 迭代法计算最大电流
|
||||
i_max = 0
|
||||
insulator_c_len = para.insulator_c_len
|
||||
i_min = min_i(insulator_c_len, u_ph / 1.732)
|
||||
# i_min = min_i(insulator_c_len, u_ph / 1.732)
|
||||
# TODO 需要考虑交、直流
|
||||
i_min = min_i(insulator_c_len, u_ph)
|
||||
_min_i = i_min # 尝试的最小电流
|
||||
_max_i = para.max_i # 尝试的最大电流
|
||||
# cad.draw(i_min, u_ph, rs_x, rs_y, rc_x, rc_y, rg_x, rg_y, rg_type, 2)
|
||||
for i_bar in np.linspace(
|
||||
_min_i, _max_i, int((_max_i - _min_i) / 0.1)
|
||||
_min_i, _max_i, int((_max_i - _min_i) / 1)
|
||||
): # 雷电流
|
||||
# logger.info(f"尝试计算电流为{i_bar:.2f}")
|
||||
logger.info(f"尝试计算电流为{i_bar:.2f}")
|
||||
rs = rs_fun(i_bar)
|
||||
animate.add_rs(rs, rs_x, rs_y)
|
||||
rc = rc_fun(i_bar, u_ph)
|
||||
animate.add_rc(rc, rc_x, rc_y)
|
||||
rg = rg_fun(i_bar, rc_y, u_ph, typ=rg_type)
|
||||
rg_line_func = None
|
||||
if rg_type == "g":
|
||||
rg_line_func = rg_line_function_factory(rg, ground_angel)
|
||||
#######
|
||||
# cccCount += 1
|
||||
# if cccCount % 30 == 0:
|
||||
# import core
|
||||
#
|
||||
# core.gMSP.add_circle((0, h_gav), rs)
|
||||
# core.gMSP.add_circle(
|
||||
# (dgc, h_cav), rc_fun(i_bar, -u_ph), dxfattribs={"color": 4}
|
||||
# )
|
||||
# core.gMSP.add_circle((dgc, h_cav), rc)
|
||||
#######
|
||||
rg_rc_circle_intersection = solve_circle_intersection(
|
||||
animate.add_rg_line(rg_line_func)
|
||||
rs_rc_circle_intersection = solve_circle_intersection(
|
||||
rs, rc, rs_x, rs_y, rc_x, rc_y
|
||||
)
|
||||
i_max = i_bar
|
||||
if not rg_rc_circle_intersection: # if circle_intersection is []
|
||||
if not rs_rc_circle_intersection: # if circle_intersection is []
|
||||
logger.debug("保护弧和暴露弧无交点，检查设置参数。")
|
||||
continue
|
||||
circle_rc_line_or_rg_intersection = None
|
||||
circle_rc_or_rg_line_intersection = None
|
||||
if rg_type == "g":
|
||||
circle_rc_line_or_rg_intersection = (
|
||||
circle_rc_or_rg_line_intersection = (
|
||||
solve_circle_line_intersection(rc, rc_x, rc_y, rg_line_func)
|
||||
)
|
||||
elif rg_type == "c":
|
||||
circle_rc_line_or_rg_intersection = solve_circle_intersection(
|
||||
circle_rc_or_rg_line_intersection = solve_circle_intersection(
|
||||
rg, rc, rg_x, rg_y, rc_x, rc_y
|
||||
)
|
||||
if not circle_rc_line_or_rg_intersection:
|
||||
if not circle_rc_or_rg_line_intersection:
|
||||
# 暴露弧和捕捉弧无交点
|
||||
if rg_type == "g":
|
||||
if rg_line_func(rc_x) > rc_y:
|
||||
i_min = i_bar # 用于后面判断最小和最大电流是否相等，相等意味着暴露弧一直被屏蔽
|
||||
logger.info(f"捕捉面在暴露弧之上，设置最小电流为{i_min:.2f}")
|
||||
logger.info(
|
||||
f"捕捉面在暴露弧之上，设置最小电流为{i_min:.2f}"
|
||||
)
|
||||
else:
|
||||
logger.info("暴露弧和地面捕捉弧无交点，检查设置参数。")
|
||||
continue
|
||||
else:
|
||||
logger.info("上面的导地线无法保护下面的导地线，检查设置参数。")
|
||||
logger.info(
|
||||
"上面的导地线无法保护下面的导地线，检查设置参数。"
|
||||
)
|
||||
continue
|
||||
animate.add_expose_area(
|
||||
rc_x,
|
||||
rc_y,
|
||||
*rs_rc_circle_intersection,
|
||||
*circle_rc_or_rg_line_intersection,
|
||||
)
|
||||
cad = Draw()
|
||||
cad.draw(
|
||||
i_min,
|
||||
u_ph,
|
||||
rs_x,
|
||||
rs_y,
|
||||
rc_x,
|
||||
rc_y,
|
||||
rg_x,
|
||||
rg_y,
|
||||
rg_type,
|
||||
ground_angel,
|
||||
2,
|
||||
) # 最小电流时
|
||||
cad.draw(
|
||||
i_max,
|
||||
u_ph,
|
||||
rs_x,
|
||||
rs_y,
|
||||
rc_x,
|
||||
rc_y,
|
||||
rg_x,
|
||||
rg_y,
|
||||
rg_type,
|
||||
ground_angel,
|
||||
6,
|
||||
) # 最大电流时
|
||||
cad.save_as(f"egm{phase_conductor_foo + 1}.dxf")
|
||||
min_distance_intersection = (
|
||||
np.sum(
|
||||
(
|
||||
np.array(rg_rc_circle_intersection)
|
||||
- np.array(circle_rc_line_or_rg_intersection)
|
||||
np.array(rs_rc_circle_intersection)
|
||||
- np.array(circle_rc_or_rg_line_intersection)
|
||||
)
|
||||
** 2
|
||||
)
|
||||
|
|
@ -210,8 +250,8 @@ def egm():
|
|||
break # 已经找到了最大电流
|
||||
# 判断是否以完全被保护
|
||||
if (
|
||||
rg_rc_circle_intersection[1]
|
||||
< circle_rc_line_or_rg_intersection[1]
|
||||
rs_rc_circle_intersection[1]
|
||||
< circle_rc_or_rg_line_intersection[1]
|
||||
):
|
||||
circle_rs_line_or_rg_intersection = None
|
||||
if rg_type == "g":
|
||||
|
|
@ -241,43 +281,16 @@ def egm():
|
|||
logger.info(f"电流为{i_bar}kV时，暴露弧已经完全被屏蔽")
|
||||
exposed_curve_shielded = True
|
||||
break
|
||||
# if phase_conductor_foo == 2:
|
||||
cad.draw(
|
||||
i_min,
|
||||
u_ph,
|
||||
rs_x,
|
||||
rs_y,
|
||||
rc_x,
|
||||
rc_y,
|
||||
rg_x,
|
||||
rg_y,
|
||||
rg_type,
|
||||
ground_angel,
|
||||
2,
|
||||
)
|
||||
cad.draw(
|
||||
i_max,
|
||||
u_ph,
|
||||
rs_x,
|
||||
rs_y,
|
||||
rc_x,
|
||||
rc_y,
|
||||
rg_x,
|
||||
rg_y,
|
||||
rg_type,
|
||||
ground_angel,
|
||||
6,
|
||||
)
|
||||
cad.save_as(f"egm{phase_conductor_foo + 1}.dxf")
|
||||
animate.pause()
|
||||
# 判断是否导线已经被完全保护
|
||||
if abs(i_max - _max_i) < 1e-5:
|
||||
logger.info("无法找到最大电流，可能是杆塔较高。")
|
||||
logger.info(f"最大电流设置为自然界最大电流{i_max}kA")
|
||||
logger.info(f"最大电流为{i_max:.2f}")
|
||||
logger.info(f"最小电流为{i_min:.2f}")
|
||||
if exposed_curve_shielded:
|
||||
logger.info("暴露弧已经完全被屏蔽,不会跳闸。")
|
||||
continue
|
||||
# if exposed_curve_shielded:
|
||||
# logger.info("暴露弧已经完全被屏蔽,不会跳闸。")
|
||||
# continue
|
||||
curt_fineness = 0.1 # 电流积分细度
|
||||
if i_min > i_max or abs(i_min - i_max) < curt_fineness:
|
||||
logger.info("最大电流小于等于最小电流，没有暴露弧。")
|
||||
|
|
@ -291,21 +304,22 @@ def egm():
|
|||
td = para.td
|
||||
ip_a = para.Ip_a
|
||||
ip_b = para.Ip_b
|
||||
cal_bd_np = (
|
||||
bd_area_vec(
|
||||
i_curt_samples,
|
||||
u_ph,
|
||||
rc_x,
|
||||
rc_y,
|
||||
rs_x,
|
||||
rs_y,
|
||||
rg_x,
|
||||
rg_y,
|
||||
ground_angel,
|
||||
rg_type,
|
||||
)
|
||||
* thunder_density(i_curt_samples, td, ip_a, ip_b)
|
||||
bd_area_vec_result = bd_area_vec(
|
||||
i_curt_samples,
|
||||
u_ph,
|
||||
rc_x,
|
||||
rc_y,
|
||||
rs_x,
|
||||
rs_y,
|
||||
rg_x,
|
||||
rg_y,
|
||||
ground_angel,
|
||||
rg_type,
|
||||
)
|
||||
thunder_density_result = thunder_density(
|
||||
i_curt_samples, td, ip_a, ip_b
|
||||
) # 雷电流幅值密度函数
|
||||
cal_bd_np = bd_area_vec_result * thunder_density_result
|
||||
calculus = np.sum(cal_bd_np[:-1] + cal_bd_np[1:]) / 2 * d_curt
|
||||
# for i_curt in i_curt_samples[:-1]:
|
||||
# cal_bd_first = bd_area(i_curt, u_ph, dgc, h_gav, h_cav)
|
||||
|
|
@ -333,10 +347,10 @@ def egm():
|
|||
)
|
||||
avr_n_sf += n_sf / voltage_n
|
||||
n_sf_phases[phase_conductor_foo][u_bar] = n_sf
|
||||
logger.info(f"工作电压为{u_ph:.2f}kV时,跳闸率是{n_sf:.16f}次/(km·a)")
|
||||
logger.info(f"线路跳闸率是{avr_n_sf:.16f}次/(km·a)")
|
||||
logger.info(f"工作电压为{u_ph:.2f}kV时,跳闸率是{n_sf:.16f}次/(100km·a)")
|
||||
logger.info(f"线路跳闸率是{avr_n_sf:.16f}次/(100km·a)")
|
||||
logger.info(
|
||||
f"不同相跳闸率是{np.array2string(np.mean(n_sf_phases,axis=1),precision=16)}次/(km·a)"
|
||||
f"不同相跳闸率是{np.array2string(np.mean(n_sf_phases,axis=1),precision=16)}次/(100km·a)"
|
||||
)
|
||||
|
||||
|
||||
|
|
@ -349,6 +363,8 @@ def speed():
|
|||
if __name__ == "__main__":
|
||||
logger.remove()
|
||||
logger.add(sys.stderr, level="DEBUG")
|
||||
run_time = timeit.timeit("egm()", globals=globals(), number=1)
|
||||
print(f"运行时间:{run_time:.2f}s")
|
||||
print("Finished.")
|
||||
egm()
|
||||
# run_time = timeit.timeit("egm()", globals=globals(), number=1)
|
||||
# logger.info(f"运行时间:{run_time:.2f}s")
|
||||
# input('enter any key to exit.')
|
||||
logger.info("Finished.")
|
||||
|
|
|
|||
|
|
@ -0,0 +1,100 @@
|
|||
import matplotlib
|
||||
from plot_data import *
|
||||
import matplotlib.pyplot as plt
|
||||
import matplotlib.ticker as mticker
|
||||
|
||||
matplotlib.use("Qt5Agg")
|
||||
# 解决中文乱码
|
||||
plt.rcParams["font.sans-serif"] = ["simsun"]
|
||||
plt.rcParams["font.family"] = "sans-serif"
|
||||
# plt.rcParams["font.weight"] = "bold"
|
||||
# 解决负号无法显示的问题
|
||||
plt.rcParams["axes.unicode_minus"] = False
|
||||
plt.rcParams["savefig.dpi"] = 1200 # 图片像素
|
||||
# plt.savefig("port.png", dpi=600, bbox_inches="tight")
|
||||
fontsize = 12
|
||||
################################################
|
||||
witdh_of_bar=0.3
|
||||
color=plt.cm.BuPu(np.linspace(152/255, 251/255,152/255))
|
||||
percent1 = data_150m塔高_不同地线保护角[:, 1] / data_150m塔高_不同地线保护角[:, 0]
|
||||
# percent1 = data_66m串长_不同塔高[:, 1] / data_66m串长_不同塔高[:, 0]
|
||||
# percent2 = data_68m串长_不同塔高[:, 1] / data_68m串长_不同塔高[:, 0]
|
||||
fig, ax = plt.subplots()
|
||||
x = np.arange(len(category_names_150m塔高_不同地线保护角)) # the label locations
|
||||
p1 = ax.bar(category_names_150m塔高_不同地线保护角, percent1, witdh_of_bar, label="绕击/反击跳闸率比值",color=color,hatch='-')
|
||||
# p1 = ax.bar(x - 0.3 / 2, percent1, 0.3, label="6.6m绝缘距离")
|
||||
# p2 = ax.bar(x + 0.3 / 2, percent2, 0.3, label="6.8m绝缘距离")
|
||||
ax.xaxis.set_major_locator(mticker.FixedLocator(x))
|
||||
ax.set_xticklabels(category_names_150m塔高_不同地线保护角)
|
||||
ax.set_ylabel("比值", fontsize=fontsize)
|
||||
ax.set_xlabel("地线保护角(°)", fontsize=fontsize)
|
||||
# ax.set_xlabel("接地电阻(Ω)", fontsize=fontsize)
|
||||
plt.xticks(fontsize=fontsize)
|
||||
plt.yticks(fontsize=fontsize)
|
||||
ax.bar_label(p1, padding=0, fontsize=fontsize)
|
||||
# ax.bar_label(p2, padding=0, fontsize=fontsize)
|
||||
ax.legend(fontsize=fontsize)
|
||||
|
||||
fig.tight_layout()
|
||||
plt.show()
|
||||
|
||||
|
||||
# results = {
|
||||
# "100m": 100 * data[0, :] / np.sum(data[0, :]),
|
||||
# "110m": data[1, :] / np.sum(data[1, :]),
|
||||
# "120m": data[2, :] / np.sum(data[2, :]),
|
||||
# "130m": data[3, :] / np.sum(data[3, :]),
|
||||
# "140m": data[4, :] / np.sum(data[4, :]),
|
||||
# "150m": data[5, :] / np.sum(data[5, :]),
|
||||
# }
|
||||
|
||||
|
||||
# def survey(results, category_names):
|
||||
# """
|
||||
# Parameters
|
||||
# ----------
|
||||
# results : dict
|
||||
# A mapping from question labels to a list of answers per category.
|
||||
# It is assumed all lists contain the same number of entries and that
|
||||
# it matches the length of *category_names*.
|
||||
# category_names : list of str
|
||||
# The category labels.
|
||||
# """
|
||||
# labels = list(results.keys())
|
||||
# data = np.array(list(results.values()))
|
||||
# data_cum = data.cumsum(axis=1)
|
||||
# category_colors = plt.get_cmap("RdYlGn")(np.linspace(0.15, 0.85, data.shape[1]))
|
||||
#
|
||||
# fig, ax = plt.subplots(figsize=(9.2, 5))
|
||||
# ax.invert_yaxis()
|
||||
# ax.xaxis.set_visible(False)
|
||||
# ax.set_xlim(0, np.sum(data, axis=1).max())
|
||||
#
|
||||
# for i, (colname, color) in enumerate(zip(category_names, category_colors)):
|
||||
# widths = data[:, i]
|
||||
# starts = data_cum[:, i] - widths
|
||||
# rects = ax.barh(
|
||||
# labels, widths, left=starts, height=0.5, label=colname, color=color
|
||||
# )
|
||||
#
|
||||
# r, g, b, _ = color
|
||||
# text_color = "white" if r * g * b < 0.5 else "darkgrey"
|
||||
# ax.bar_label(rects, label_type="center", color=text_color)
|
||||
# ax.legend(
|
||||
# ncol=len(category_names),
|
||||
# bbox_to_anchor=(0, 1),
|
||||
# loc="lower left",
|
||||
# fontsize="small",
|
||||
# )
|
||||
#
|
||||
# return fig, ax
|
||||
|
||||
# percent=data/np.sum(data,axis=1)[:,None]*100
|
||||
# percent = data[:, 1] / data[:, 0]
|
||||
# plt.bar(category_names, percent, 0.3, label="黑")
|
||||
# # plt.bar(category_names, percent[:,0], 0.2, label="r")
|
||||
#
|
||||
# # plt.bar(category_names, [0.014094 / 100, 0.025094 / 100], 0.2, label="h")
|
||||
# plt.legend()
|
||||
# # survey(results, category_names)
|
||||
# plt.show()
|
||||
Loading…
Reference in New Issue