egm/main.py

import math
import os.path
import sys
import tomli
from loguru import logger
from core import *
import timeit


# 打印参数
def parameter_display(para_dis: Parameter):
    logger.info(f"额定电压 kV {para_dis.rated_voltage}")
    logger.info(f"导线弧垂 m {para_dis.h_c_sag}")
    logger.info(f"地线弧垂 m {para_dis.h_g_sag}")
    logger.info(f"全塔高 m {para_dis.h_arm[0]}")
    logger.info(f"串绝缘距离 m {para_dis.insulator_c_len}")
    logger.info(f"导线串长 m {para_dis.string_c_len}")
    logger.info(f"地线串长 m {para_dis.string_g_len}")
    logger.info(f"挂点垂直坐标 m {para_dis.h_arm}")
    logger.info(f"挂点水平坐标 m {para_dis.gc_x}")
    logger.info(f"地面倾角 ° {[an * 180 / math.pi for an in para_dis.ground_angels]}")
    logger.info(f"海拔高度 m {para_dis.altitude}")
    if para_dis.ng > 0:
        logger.info("不采用雷暴日计算地闪密度和雷电流密度")
        logger.info(f"地闪密度 次/(每平方公里·每年) {para_dis.ng}")
        logger.info(f"概率密度曲线系数a {para_dis.Ip_a}")
        logger.info(f"概率密度曲线系数b {para_dis.Ip_b}")
        pass
    else:
        logger.info(f"雷暴日 d {para_dis.td}")


def read_parameter(toml_file_path):
    with open(toml_file_path, "rb") as toml_fs:
        toml_dict = tomli.load(toml_fs)
        toml_parameter = toml_dict["parameter"]
        para.h_g_sag = toml_parameter["h_g_sag"]  # 地线弧垂
        para.h_c_sag = toml_parameter["h_c_sag"]  # 导线弧垂
        # para.h_whole = toml_parameter["h_whole"]  # 杆塔全高
        para.td = toml_parameter["td"]  # 雷暴日
        para.insulator_c_len = toml_parameter["insulator_c_len"]  # 串子绝缘长度
        para.string_c_len = toml_parameter["string_c_len"]
        para.string_g_len = toml_parameter["string_g_len"]
        para.gc_x = toml_parameter["gc_x"]  # 导、地线水平坐标
        para.ground_angels = [
            angel / 180 * math.pi for angel in toml_parameter["ground_angels"]
        ]  # 地面倾角，向下为正
        para.h_arm = toml_parameter["h_arm"]
        para.altitude = toml_parameter["altitude"]
        para.rated_voltage = toml_parameter["rated_voltage"]
        toml_advance = toml_dict["advance"]
        para.ng = toml_advance["ng"]  # 地闪密度
        para.Ip_a = toml_advance["Ip_a"]  # 概率密度曲线系数a
        para.Ip_b = toml_advance["Ip_b"]  # 概率密度曲线系数b
        toml_optional = toml_dict["optional"]
        para.voltage_n = toml_optional["voltage_n"]  # 工作电压分成多少份来计算
        para.max_i = toml_optional["max_i"]


def egm():
    if len(sys.argv) < 2:
        toml_file_path = r"article.toml"
    else:
        toml_file_path = sys.argv[1]
    if not os.path.exists(toml_file_path):
        logger.info(f"无法找到数据文件{toml_file_path}，程序退出。")
        sys.exit(0)
    logger.info(f"读取文件{toml_file_path}")
    read_parameter(toml_file_path)
    #########################################################
    # 以上是需要设置的参数
    parameter_display(para)
    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
    h_c_sag = para.h_c_sag
    gc_x = para.gc_x
    h_arm = para.h_arm
    gc_y = [
        h_whole - string_g_len - h_g_sag * 2 / 3,  # 地线对地平均高
    ]
    if len(h_arm) > 1:
        for hoo in h_arm[1:]:
            gc_y.append(hoo - string_c_len - h_c_sag * 2 / 3)
    if len(gc_y) > 2:  # 双回路
        phase_n = 3  # 边相导线数量
    else:
        phase_n = 1
    # 地闪密度 利用Q╱GDW 11452-2015 架空输电线路防雷导则的公式 Ng=0.023*Td^(1.3) 20天雷暴日地闪密度为1.13
    td = para.td
    ng = func_ng(td)
    avr_n_sf = 0  # 考虑电压的影响计算的跳闸率
    ground_angels = para.ground_angels
    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):
            logger.info("导线可能掉地面下了，程序退出。")
            return 0
        for phase_conductor_foo in range(phase_n):
            exposed_curve_shielded = False
            rs_x = gc_x[phase_conductor_foo]
            rs_y = gc_y[phase_conductor_foo]
            rc_x = gc_x[phase_conductor_foo + 1]
            rc_y = gc_y[phase_conductor_foo + 1]
            if phase_n == 1:
                rg_type = "g"
            if phase_n > 1:  # 多回路
                if phase_conductor_foo < 2:
                    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 = (
                math.atan(
                    (rc_x - rs_x)
                    / (
                        (h_arm[0] - string_g_len - h_arm[phase_conductor_foo + 1])
                        + string_c_len
                    )
                )
                * 180
                / math.pi
            )  # 保护角
            logger.info(f"保护角{shield_angle:.3f}°")
            logger.debug(f"最低相防护标识{rg_type}")
            for u_bar in range(voltage_n):  # 计算不同工作电压下的跳闸率
                u_ph = (
                    math.sqrt(2)
                    * 750
                    * math.cos(2 * math.pi / voltage_n * u_bar)
                    / 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)
                _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)
                ):  # 雷电流
                    # logger.info(f"尝试计算电流为{i_bar:.2f}")
                    rs = rs_fun(i_bar)
                    rc = rc_fun(i_bar, u_ph)
                    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(
                        rs, rc, rs_x, rs_y, rc_x, rc_y
                    )
                    i_max = i_bar
                    if not rg_rc_circle_intersection:  # if circle_intersection is []
                        logger.debug("保护弧和暴露弧无交点，检查设置参数。")
                        continue
                    circle_rc_line_or_rg_intersection = None
                    if rg_type == "g":
                        circle_rc_line_or_rg_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(
                            rg, rc, rg_x, rg_y, rc_x, rc_y
                        )
                    if not circle_rc_line_or_rg_intersection:
                        # 暴露弧和捕捉弧无交点
                        if rg_type == "g":
                            if rg_line_func(rc_x) > rc_y:
                                i_min = i_bar  # 用于后面判断最小和最大电流是否相等，相等意味着暴露弧一直被屏蔽
                                logger.info(f"捕捉面在暴露弧之上，设置最小电流为{i_min:.2f}")
                            else:
                                logger.info("暴露弧和地面捕捉弧无交点，检查设置参数。")
                            continue
                        else:
                            logger.info("上面的导地线无法保护下面的导地线，检查设置参数。")
                            continue
                    min_distance_intersection = (
                        np.sum(
                            (
                                np.array(rg_rc_circle_intersection)
                                - np.array(circle_rc_line_or_rg_intersection)
                            )
                            ** 2
                        )
                        ** 0.5
                    )  # 计算两圆交点和地面直线交点的最小距离
                    if min_distance_intersection < 0.1:
                        break  # 已经找到了最大电流
                    # 判断是否以完全被保护
                    if (
                        rg_rc_circle_intersection[1]
                        < circle_rc_line_or_rg_intersection[1]
                    ):
                        circle_rs_line_or_rg_intersection = None
                        if rg_type == "g":
                            circle_rs_line_or_rg_intersection = (
                                solve_circle_line_intersection(
                                    rs, rs_x, rs_y, rg_line_func
                                )  # 保护弧和捕雷弧的交点
                            )
                        if rg_type == "c":
                            circle_rs_line_or_rg_intersection = (
                                solve_circle_intersection(
                                    rs, rg, rs_x, rs_y, rg_x, rg_y
                                )
                            )
                        # 判断与保护弧的交点是否在暴露弧外面
                        distance = (
                            np.sum(
                                (
                                    np.array(circle_rs_line_or_rg_intersection)
                                    - np.array([rc_x, rc_y])
                                )
                                ** 2
                            )
                            ** 0.5
                        )
                        if distance > rc:
                            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")
                # 判断是否导线已经被完全保护
                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
                curt_fineness = 0.1  # 电流积分细度
                if i_min > i_max or abs(i_min - i_max) < curt_fineness:
                    logger.info("最大电流小于等于最小电流，没有暴露弧。")
                    continue
                # 开始积分
                curt_segment_n = int((i_max - i_min) / curt_fineness)  # 分成多少份
                i_curt_samples, d_curt = np.linspace(
                    i_min, i_max, curt_segment_n + 1, retstep=True
                )
                bd_area_vec = np.vectorize(bd_area)
                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)
                )
                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)
                #     cal_bd_second = bd_area(i_curt + d_curt, u_ph, dgc, h_gav, h_cav)
                #     cal_thunder_density_first = thunder_density(i_curt)
                #     cal_thunder_density_second = thunder_density(i_curt + d_curt)
                #     calculus += (
                #         (
                #             cal_bd_first * cal_thunder_density_first
                #             + cal_bd_second * cal_thunder_density_second
                #         )
                #         / 2
                #         * d_curt
                #     )
                #     if abs(calculus-0.05812740052770032)<1e-5:
                #         abc=123
                #         pass
                rated_voltage = para.rated_voltage
                n_sf = (
                    2
                    * ng
                    / 10
                    * calculus
                    * arc_possibility(rated_voltage, insulator_c_len)
                )
                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"不同相跳闸率是{np.array2string(np.mean(n_sf_phases,axis=1),precision=16)}次/(km·a)"
        )


def speed():
    a = 0
    for bar in range(100000000):
        a += bar


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.")