unbanlanced_tension/auto_differentiation.py

# 利用自动微分计算
import sympy
import data
import exp
import math
import main

sympy.init_printing()
# h_i 悬点高差
# l_i 悬点档距
# _alpha 导线膨胀系数  1/°C
# _elastic 弹性系数 N/mm2
# _t_e 架线时考虑初伸长的降温，取正值。单位°C
# lambda_i 计算不平衡张力时导线比载 N/(m.mm)
# sigma_i 计算不平衡张力时最低点水平应力 单位N/mm2
# t_i 计算不平衡张力时导线温度 单位°C
# _lambda_m 导线架线时时导线比载 N/(m.mm)
# _sigma_m 导线架线时时最低点水平应力 单位N/mm2
# _t_m 导线架线时时导线温度 单位°C


delta_Li__1 = sympy.symbols(
    "delta_Li:{span_count}".format(span_count=data.span_count - 1)
)
delta_Li = (
    *delta_Li__1,
    sympy.symbols("delta_Li_i"),
)
# sigma_i = sympy.symbols("sigma_i:{span_count}".format(span_count=data.span_count))


loop_end = data.loop_end  # 最大循环次数
# 架线时的状态
# 取外过无风
string_length = data.string_length  # 串长 单位m
string_g = data.string_g  # 串重 单位N
t_m_data = data.t_m  # 导线架设时的气温。单位°C
t_e_data = data.t_e  # 架线时考虑初伸长的降温，取正值。单位°C
alpha_data = data.alpha  # 导线膨胀系数  1/°C
elastic_data = data.elastic  # 弹性系数 N/mm2
area_data = data.area  # 导线面积 mm2
lambda_m_data = data.lambda_m  # 导线比载 N/(m.mm)
sigma_m_data = data.sigma_m  # 架线时，初伸长未释放前的最低点水平应力。单位N/mm2
span_count = data.span_count  # 几个档距
# n个档距,n-1个直线塔
h_array = data.h_array
hi_matrix = sympy.Matrix(h_array)
# sympy.pprint(hi_matrix)
l_array = data.l_array
l_matrix = sympy.Matrix(l_array)
t_data = data.t
conductor_n = data.conductor_n
# ti_matrix = sympy.Matrix(t_array)
lambda_i_array = data.lambda_i_array
# TODO: 暂时没考虑荷载变化
lambda_m_matrix = sympy.Matrix(lambda_i_array)
lambda_i_matrix = sympy.Matrix(lambda_i_array)


symbol_delta_l_i = exp.delta_li()
sigma_i = sympy.symbols("sigma_i")
d_delta_l_i_sigma_i = sympy.diff(symbol_delta_l_i, sigma_i)
symbol_sigma_i1 = exp.fun_sigma_i1(delta_Li)
d_sigma_i1_sigma_i = sympy.diff(symbol_sigma_i1, sigma_i)
# sigma_i1 = sympy.symbols("sigma_i1")
# d_sigma_i1_sigma_i1 = sympy.diff(symbol_sigma_i1, sigma_i)
delta_Li_i = sympy.symbols("delta_Li_i")
d_sigma_i1_d_l_i = sympy.diff(symbol_sigma_i1, delta_Li_i)


# 一共2n个变量，n个delta_Li，n个sigma_i
# 分 [
#     A B
#     C D
#     E1 E2
#   ]
# 6块


# B为dΔli/dσi
def evaluate_d_delta_l_i_sigma_i(val_delta_l_li, val_sigma_i):
    (
        delta_l_i,
        l_i,
        lambda_i,
        alpha,
        E,
        t_e,
        t_i,
        lambda_m,
        t_m,
        sigma_m,
        _sigma_i,
        beta_i,
    ) = sympy.symbols(
        """ 
        delta_l_i,
        l_i,
        lambda_i,
        alpha,
        E,
        t_e,
        t_i,
        lambda_m,
        t_m,
        sigma_m,
        sigma_i,
        beta_i
        """
    )
    val_list = []
    for i in range(span_count):
        val = d_delta_l_i_sigma_i.subs(
            [
                (delta_l_i, val_delta_l_li[i]),
                (l_i, l_array[i]),
                (lambda_i, lambda_i_array[i]),
                (alpha, alpha_data),
                (E, elastic_data),
                (t_e, t_e_data),
                (t_i, t_data),
                (lambda_m, lambda_m_data),
                (t_m, t_m_data),
                (sigma_m, sigma_m_data),
                (_sigma_i, val_sigma_i[i]),
                (beta_i, math.atan(h_array[i] / l_array[i])),
            ]
        )
        val_list.append(val)
    return val_list


# C为dσi1dΔli
# C只有n-1行
def evaluate_d_sigma_i1_d_delta_l_i(val_delta_l_li, val_sigma_i):
    (
        G_i,
        A,
        lambda_i,
        lambda_i1,
        _sigma_i,
        h_i,
        h_i1,
        l_i,
        l_i1,
        stringlen_i,
        _sigma_i1,
        beta_i,
        beta_i1,
    ) = sympy.symbols(
        """
        G_i,
        A,
        lambda_i,
        lambda_i1,
        sigma_i,
        h_i,
        h_i1,
        l_i,
        l_i1,
        stringlen_i,
        sigma_i1,
        beta_i,
        beta_i1,
        """
    )
    row = []
    for i in range(span_count - 1):
        col = []
        for j in range(span_count):
            if i < j:
                col.append(0)
            else:
                _val = d_sigma_i1_d_l_i.subs(
                    [
                        (G_i, string_g / conductor_n),
                        (A, area_data),
                        (lambda_i, lambda_i_array[i]),
                        (lambda_i1, lambda_i_array[i + 1]),
                        (_sigma_i, val_sigma_i[i]),
                        (h_i, h_array[i]),
                        (h_i1, h_array[i + 1]),
                        (l_i, l_array[i]),
                        (l_i1, l_array[i + 1]),
                        (stringlen_i, string_length),
                        (_sigma_i1, val_sigma_i[i + 1]),
                        (beta_i, math.atan(h_array[i] / l_array[i])),
                        (beta_i1, math.atan(h_array[i + 1] / l_array[i + 1])),
                    ]
                )
                _val_delta_l_li = list(val_delta_l_li)
                _val_delta_l_li[-1] = _val_delta_l_li[j]  # 把需要求导的Δlj放最后一个位置
                _val_delta_l_li[j] = 0
                # σi1的第i+1行至倒数第2行全部清0
                for k in range(i + 1, len(_val_delta_l_li) - 1):
                    _val_delta_l_li[k] = 0
                # if index == i:
                #     _val = _val.subs(li, val_delta_l_li[index])
                # if index > i:
                #     _val = _val.subs(li, 0)
                for index, li in enumerate(delta_Li):
                    _val = _val.subs(li, _val_delta_l_li[index])
                    pass
                col.append(_val)
        row.append(col)
    return sympy.Matrix(row)


# D为dΔσi1dσi
# D只有n-1行
def evaluate_d_sigma_i1_d_delta_sigma_i(val_delta_li):
    row = []
    for i in range(span_count - 1):
        col = []
        for j in range(span_count):
            if i == j:
                sum_delta_li = math.fsum(val_delta_li)
                _val = -(
                    (
                        h_array[i] / l_array[i]
                        + ((string_g / conductor_n) ** 2 - sum_delta_li ** 2) ** 0.5
                    )
                    / (
                        ((string_g / conductor_n) ** 2 - sum_delta_li ** 2) ** 0.5
                        + h_array[i + 1] / l_array[i + 1]
                    )
                )
                col.append(_val)
                continue
            if i == j - 1:
                col.append(1)
                continue
            col.append(0)
        row.append(col)
    return sympy.Matrix(row)


def solve():
    # 初始化
    val_delta_li = [0.1 for i in range(span_count)]
    # val_delta_li = [0.15864687475316822, -0.1935189734784845, 0.03478489898855073]

    val_sigma_i = [sigma_m_data for _ in range(span_count)]
    # val_sigma_i = [175.38451579479482, 176.01015153076614, 175.88355419459572]

    loop = 0
    while True:
        loop += 1
        print("第{loop}次迭代".format(loop=loop))
        if loop >= 20:
            break
        # A为dΔli/dli
        M_A = sympy.eye(span_count)
        # B为dΔli/dσi
        M_B = sympy.diag(
            evaluate_d_delta_l_i_sigma_i(val_delta_li, val_sigma_i), unpack=True
        )
        # C为dΔσi1dli
        M_C = evaluate_d_sigma_i1_d_delta_l_i(val_delta_li, val_sigma_i)
        # D为dΔσi1dσi
        M_D = evaluate_d_sigma_i1_d_delta_sigma_i(val_delta_li)
        E1 = [1 for _ in range(span_count)]
        E2 = [0 for _ in range(span_count)]
        E = list(E1)
        E.extend(E2)
        M_E = sympy.Matrix([E])
        # 解方程
        A = sympy.Matrix([[M_A, M_B], [M_C, M_D], [M_E]])
        fx_delta_Li = []
        fx_sigma_i1 = []
        b_i = 0
        for i in range(span_count):
            fx_delta_Li.append(
                val_delta_li[i]
                - main.delta_li(
                    h_array[i],
                    l_array[i],
                    lambda_i_array[i],
                    alpha_data,
                    elastic_data,
                    t_e_data,
                    t_data,
                    val_sigma_i[i],
                    lambda_m_data,
                    t_m_data,
                    sigma_m_data,
                )
            )
            if i < span_count - 1:
                fx_sigma_i1.append(
                    val_sigma_i[i + 1]
                    - main.fun_sigma_i1(
                        area_data,
                        val_sigma_i[i],
                        math.fsum(val_delta_li[0 : i + 1]),
                        string_length,
                        string_g / conductor_n,
                        h_array[i],
                        l_array[i],
                        lambda_i_array[i],
                        h_array[i + 1],
                        l_array[i + 1],
                        lambda_i_array[i + 1],
                    )
                )

                # lambda_i1 = lambda_i_array[i + 1]
                # h_i1 = h_array[i + 1]
                # l_i1 = l_array[i + 1]
                # h_i = h_array[i]
                # l_i = l_array[i]
                # beta_i = math.atan(h_i / l_i)
                # beta_i1 = math.atan(h_i1 / l_i1)
                # w_i = (
                #         lambda_i_array[i] * l_array[i] / 2 / math.cos(beta_i)
                #         + val_sigma_i[i] * h_i / l_i
                #         + (lambda_i1 * l_i1 / 2 / math.cos(beta_i1) - val_sigma_i[i+1] * h_i1 / l_i1)
                # )
                # b_i += val_delta_li[i]
                # # 新版大手册p329 (5-61) 最上方公式
                # right_equ = val_sigma_i[i] + b_i / math.sqrt(string_length ** 2 - b_i ** 2) * (
                #         string_g/conductor_n / 2 / area + w_i  # string_g已在传入时考虑了导线分裂数
                # )

        fx_sum_Li = [math.fsum(val_delta_li)]
        b_list = []
        b_list.extend(fx_delta_Li)
        b_list.extend(fx_sigma_i1)
        b_list.extend(fx_sum_Li)
        b = sympy.Matrix(b_list)
        # sympy.pprint(b)
        x = sympy.linsolve((-A, b))
        x_list = list(x)[0]
        abs_min = [math.fabs(_x) for _x in x_list]
        abs_min.sort()
        if abs_min[-1] < 1e-5:
            break
        print("最大偏差{max_dx}".format(max_dx=abs_min[-1]))
        # 更新变量
        for i in range(span_count):
            val_delta_li[i] += x_list[i]
            val_sigma_i[i] += x_list[i + span_count]
    if loop >= loop_end:
        print("不收敛")
    else:
        print(loop)
        print(val_delta_li)
        print(val_sigma_i)
        verify(val_delta_li,val_sigma_i)


def verify(val_delta_li, val_sigma_i):
    main.verify(
        area_data,
        h_array,
        l_array,
        string_length,
        string_g / conductor_n,
        val_sigma_i,
        val_delta_li,
        lambda_i_array,
        t_data,
        alpha_data,
        elastic_data,
        t_e_data,
        lambda_m_data,
        t_m_data,
        sigma_m_data,
        1e-5,
    )


solve()

print("Finished.")