From 7a5bb05f58b22b906ea97931841d451570f66ac7 Mon Sep 17 00:00:00 2001 From: dmy Date: Wed, 6 Nov 2024 23:32:51 +0800 Subject: [PATCH] =?UTF-8?q?1.=E6=AF=8F=E4=B8=80=E6=AC=A1=E8=AE=A1=E7=AE=97?= =?UTF-8?q?=E9=83=BD=E9=87=8D=E6=96=B0=E7=94=BB=E7=94=BB=E3=80=82?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .gitignore | 13 ++++ animation.py | 11 +-- main.py | 198 ++++++++++++++++++++++++++++----------------------- plot.py | 100 ++++++++++++++++++++++++++ 4 files changed, 227 insertions(+), 95 deletions(-) create mode 100644 .gitignore create mode 100644 plot.py diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..f805294 --- /dev/null +++ b/.gitignore @@ -0,0 +1,13 @@ +*.dxf +build +__pycache__ +CSharp +.idea +dist +*.spec +*.dwg +历史 +.venv +*.toml +launch.json +settings.json \ No newline at end of file diff --git a/animation.py b/animation.py index af63352..5f56906 100644 --- a/animation.py +++ b/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 diff --git a/main.py b/main.py index a03f3b4..157586a 100644 --- a/main.py +++ b/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.") diff --git a/plot.py b/plot.py new file mode 100644 index 0000000..bce5869 --- /dev/null +++ b/plot.py @@ -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()