egm/main.py

import sys

from loguru import logger
from core import *
import timeit


def egm():
    h_g_avr_sag = 11.67 * 2 / 3
    h_c_avr_sag = 14.43 * 2 / 3
    h_whole = 130  # 杆塔全高
    voltage_n = 3  # 工作电压分成多少份来计算
    td = 20  # 雷暴日
    insulator_c_len = 6.6  # 串子绝缘长度
    string_c_len = 9.2
    string_g_len = 0.5
    gc_x = [17.9, 17, 15, 17.0]

    # 以后考虑地形角度,地面线
    def ground_surface(x):
        return 0

    gc_y = [
        h_whole - string_g_len - h_g_avr_sag,  # 地线对地平均高
        # h_whole - string_c_len - h_c_avr_sag - 2.7,  # 导线对地平均高
        h_whole - string_c_len - h_c_avr_sag - 20,  # 导线对地平均高
        # h_whole - string_c_len - h_c_avr_sag - 35.7,  # 导线对地平均高
    ]
    if len(gc_y) > 2:  # 双回路
        phase_n = 3  # 边相导线数量
    else:
        phase_n = 1
    #########################################################
    rg_type = None
    # 以上是需要设置的参数
    avr_n_sf = 0  # 考虑电压的影响计算的跳闸率
    rg_x = None
    rg_y = None
    cad = Draw()
    # 跳闸率 利用Q╱GDW 11452-2015 架空输电线路防雷导则的公式 Ng=0.023*Td^(1.3) 20天雷暴日地闪密度为1.13
    ng = func_ng(td)
    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) / ((rs_y - rc_y) + 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
            i_min = min_i(insulator_c_len, u_ph / 1.732)
            _min_i = i_min  # 尝试的最小电流
            _max_i = 200  # 尝试的最大电流
            # 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)
                #######
                # 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, rg, rc_x, rc_y
                    )
                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 > 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, rg, rs_x, rs_y)
                        )
                    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("暴露弧已经完全被屏蔽")
                        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, 2)
            cad.draw(i_max, u_ph, rs_x, rs_y, rc_x, rc_y, rg_x, rg_y, rg_type, 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)
            cal_bd_np = (
                bd_area_vec(
                    i_curt_samples,
                    u_ph,
                    rc_x,
                    rc_y,
                    rs_x,
                    rs_y,
                    rg_x,
                    rg_y,
                    ground_surface,
                    rg_type,
                )
                * thunder_density(i_curt_samples)
            )
            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
            n_sf = 2 * ng / 10 * calculus * arc_possibility(750, 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:.6}")
    logger.info(f"跳闸率是{avr_n_sf:.6f}")
    logger.info(f"不同相跳闸率是{np.array2string(np.mean(n_sf_phases,axis=1),precision=6)}")


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