准备修改成通用程序。

2.收敛条件强制让不等式满足要求

.gitignore

@@ -1 +1,2 @@
-venv
+venv
+*.xls

auto_differentiation.py

@@ -0,0 +1,319 @@
+# 利用自动微分计算
+import datetime
+from typing import List
+
+import sympy
+import data
+import exp
+import math
+import main
+import numpy as np
+
+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
+
+
+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
+l_array = data.l_array
+t_data = data.t
+conductor_n = data.conductor_n
+lambda_i_array = data.lambda_i_array
+# TODO: 暂时没考虑荷载变化
+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)
+fx_d_delta_l_i_sigma_i = exp.get_lambdify_d_delta_l_i_sigma_i(d_delta_l_i_sigma_i)
+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"),
+)
+symbol_sigma_i1 = exp.fun_sigma_i1(delta_Li)
+delta_Li_i = sympy.symbols("delta_Li_i")
+d_sigma_i1_d_l_i = sympy.diff(symbol_sigma_i1, delta_Li_i)
+
+
+fx_d_sigma_i1_d_l_i = exp.get_lambdify_d_sigma_i1_d_l_i(d_sigma_i1_d_l_i, delta_Li)
+
+
+# 一共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):
+
+    val_list = []
+    for i in range(span_count):
+        val = sympy.Float(
+            fx_d_delta_l_i_sigma_i(
+                val_delta_l_li[i],
+                l_array[i],
+                lambda_i_array[i],
+                alpha_data,
+                elastic_data,
+                t_e_data,
+                t_data,
+                lambda_m_data,
+                t_m_data,
+                sigma_m_data,
+                val_sigma_i[i],
+                math.atan(h_array[i] / l_array[i]),
+            )
+        )
+        # manual_val = exp.manual_diff_delta_li_sigma_i(
+        #     h_array[i],
+        #     l_array[i],
+        #     lambda_m_data,
+        #     lambda_i_array[i],
+        #     val_sigma_i[i],
+        #     sigma_m_data,
+        #     alpha_data,
+        #     t_data,
+        #     t_e_data,
+        #     t_m_data,
+        #     elastic_data,
+        # )
+        # if math.fabs(val - manual_val) > 1e-5:
+        #     raise Exception("d_delta_l_i_sigma_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):
+
+    row = []
+    for i in range(span_count - 1):
+        col = []
+        for j in range(span_count):
+            if i < j:
+                col.append(0)
+            else:
+                _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
+                _val = sympy.Float(
+                    fx_d_sigma_i1_d_l_i(
+                        string_g / conductor_n,
+                        area_data,
+                        lambda_i_array[i],
+                        lambda_i_array[i + 1],
+                        val_sigma_i[i],
+                        h_array[i],
+                        h_array[i + 1],
+                        l_array[i],
+                        l_array[i + 1],
+                        string_length,
+                        val_sigma_i[i + 1],
+                        math.atan(h_array[i] / l_array[i]),
+                        math.atan(h_array[i + 1] / l_array[i + 1]),
+                        *_val_delta_l_li
+                    )
+                )
+
+                # manual_val = exp.manual_diff_sigma_i1_d_l_i(
+                #     h_array[i],
+                #     l_array[i],
+                #     h_array[i + 1],
+                #     l_array[i + 1],
+                #     string_g / conductor_n,
+                #     area_data,
+                #     lambda_i_array[i],
+                #     lambda_i_array[i + 1],
+                #     val_sigma_i[i],
+                #     string_length,
+                #     math.fsum(val_delta_l_li[0 : i + 1]),
+                # )
+                # if math.fabs(manual_val - _val) > 1e-5:
+                #     raise Exception("d_sigma_i1_delta_L_i 自动和手动微分不匹配")
+                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():
+    starttime = datetime.datetime.now()
+    # 初始化
+    val_delta_li = [data.string_length / (span_count + 1) for _ in range(span_count)]
+    val_sigma_i = [sigma_m_data for _ in range(span_count)]
+    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 = []
+        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],
+                    )
+                )
+        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)
+        AA = np.array(A.tolist(), dtype=np.float64)
+        b = np.array(b_list, dtype=np.float64)
+        x = np.linalg.solve(-AA, b)
+        x_list = x
+        # 强制要求等式满足
+        abs_min: List[float] = [
+            math.fabs(_x) for _x in [*x_list, *fx_delta_Li, *fx_sigma_i1]
+        ]
+        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}次迭代收敛,最大偏差{bias}".format(loop=loop, bias=abs_min[-1]))
+        # print(val_delta_li)
+        # print(val_sigma_i)
+        verify(val_delta_li, val_sigma_i)
+        endtime = datetime.datetime.now()
+        # print('执行时间{t}'.format(t=(endtime-starttime).microseconds/1000))
+        pass
+
+
+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.")

data.py

@@ -0,0 +1,117 @@
+# # loop_end = 10000000  # 最大循环次数
+# loop_end = 1000000  # 最大循环次数
+# # 架线时的状态
+# # 取外过无风
+# string_length = 9.2  # 串长 单位m
+# string_g = 60 * 9.8  # 串重 单位N
+# t_m = 15  # 导线架设时的气温。单位°C
+# t_e = 20  # 架线时考虑初伸长的降温，取正值。单位°C
+# alpha = 0.0000155  # 导线膨胀系数  1/°C
+# elastic = 95900  # 弹性系数 N/mm2
+# area = 154.48  # 导线面积 mm2
+# lambda_m = 14.8129 / area  # 导线比载 N/(m.mm)
+# lambda_i_array = [
+#     lambda_m * 0.9,
+#     lambda_m * 1.2,
+#     lambda_m,
+#     lambda_m,
+#     lambda_m,
+#     lambda_m * 0.9,
+#     lambda_m * 1.5,
+#     lambda_m,
+#     lambda_m,
+#     lambda_m,
+#     lambda_m,
+#     lambda_m,
+#     lambda_m,
+#     lambda_m * 0.9,
+#     lambda_m * 1.3,
+# ]
+# # 取400m代表档距下
+# sigma_m = 28517 / area  # 架线时，初伸长未释放前的最低点水平应力。单位N/mm2
+# span_count = 14  # 几个档距
+# # n个档距,n-1个直线塔
+# h_array = [
+#     0,
+#     0,
+#     0,
+#     0,
+#     0,
+#     0,
+#     0,
+#     0,
+#     0,
+#     0,
+#     0,
+#     0,
+#     0,
+#     0,
+#     0,
+# ]
+# l_array = [
+#     400,
+#     300,
+#     300,
+#     500,
+#     300,
+#     400,
+#     300,
+#     300,
+#     500,
+#     300,
+#     300,
+#     500,
+#     300,
+#     400,
+#     300,
+# ]
+# t = 15
+# epsilon = 1e-4  # 收敛判据
+# conductor_n = 6  # 导线分裂数
+
+
+loop_end = 5000000  # 最大循环次数
+# 架线时的状态
+# 取外过无风
+string_length = 0.5  # 串长 单位m
+string_g = 30 * 9.8  # 串重 单位N
+t_m = -25  # 导线架设时的气温。单位°C
+t_e = 0  # 架线时考虑初伸长的降温，取正值。单位°C
+alpha = 0.0000155  # 导线膨胀系数  1/°C
+elastic = 95900  # 弹性系数 N/mm2
+area = 154.48  # 导线面积 mm2
+lambda_m = 7.3256 / area  # 导线比载 N/(m.mm)
+lambda_i_array = [
+    14.7012 / area,
+    14.7012 / area,
+    14.7012 / area,
+    8.8007 / area,
+    14.7012 / area,
+    14.7012 / area,
+    14.7012 / area,
+]
+# 取400m代表档距下
+sigma_m = 17449 / area  # 架线时，初伸长未释放前的最低点水平应力。单位N/mm2
+span_count = 7  # 几个档距
+# n个档距,n-1个直线塔
+h_array = [
+    0,
+    0,
+    0,
+    0,
+    0,
+    0,
+    0,
+]
+l_array = [
+    703,
+    720,
+    587,
+    620,
+    539,
+    450,
+    611,
+]
+t = 40
+epsilon = 1e-4  # 收敛判据
+conductor_n = 1  # 导线分裂数

exp.py

@@ -0,0 +1,273 @@
+import sympy
+import math
+
+# h 悬点高差
+# 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
+
+
+def delta_li():
+    (
+        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,"""
+    )
+    _delta_li = delta_l_i - (
+        l_i
+        / ((sympy.cos(beta_i) ** 2) * (1 + (lambda_i * l_i / sigma_i) ** 2 / 8))
+        * (
+            (l_i * sympy.cos(beta_i)) ** 2
+            / 24
+            * ((lambda_m / sigma_m) ** 2 - (lambda_i / sigma_i) ** 2)
+            + ((sigma_i - sigma_m) / E / sympy.cos(beta_i))
+            + alpha * (t_i + t_e - t_m)
+        )
+    )
+    return _delta_li
+
+
+# area 导线截面 单位mm2
+# sigma_i 第i档内水平应力 单位N/mm2
+# b_i 悬垂串沿线路方向水平偏移距离，沿大号方向为正，反之为负。 单位m
+# stringlen_i 第i基直线塔串长 单位m
+# G_i 第i基直线塔串重 单位N
+# h_i 悬垂串处千中垂位置时，,第 i 基对第 i-1 杆塔上导线悬挂点间的高差大号比小号杆塔悬挂点高者h本身为正值，反之为负值。
+# lambda_i 第i档导线比载 N/(m.mm)
+# h_i1 悬垂串处千中垂位置时，第 i+1 基对第 i 杆塔上导线悬挂点间的高差大号比小号杆塔悬挂点高者h本身为正值，反之为负值。
+# lambda_i1
+def fun_sigma_i1(delta_Li):
+    (
+        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
+        """
+    )
+
+    def b_i():
+        _t = sympy.Float(0)
+        for f in delta_Li:
+            _t += f
+        return _t
+
+    _sigma_i1 = sigma_i1 - (
+        (
+            G_i / 2 / A  # G_i传入时已考虑导线分裂数
+            + lambda_i * l_i / 2 / sympy.cos(beta_i)
+            + lambda_i1 * l_i1 / 2 / sympy.cos(beta_i1)
+            + sigma_i * h_i / l_i
+        )
+        + sigma_i / b_i() * sympy.sqrt(stringlen_i ** 2 - b_i() ** 2)
+    ) / (sympy.sqrt(stringlen_i ** 2 - b_i() ** 2) / b_i() + h_i1 / l_i1)
+    return _sigma_i1
+
+
+def manual_diff_delta_li_sigma_i(
+    h_i, l_i, lambda_m, lambda_i, sigma_i, sigma_m, alpha, t_i, t_e, t_m, E
+):
+    beta_i = math.atan(h_i / l_i)
+    A = (
+        (l_i * math.cos(beta_i)) ** 2
+        / 24
+        * ((lambda_m / sigma_m) ** 2 - (lambda_i / sigma_i) ** 2)
+        + (sigma_i - sigma_m) / E / math.cos(beta_i)
+        + alpha * (t_i + t_e - t_m)
+    )
+    B = 1 + (lambda_i * l_i / sigma_i) ** 2 / 8
+    dA_dsigma_i = ((l_i * math.cos(beta_i)) ** 2) / 24 * 2 * lambda_i ** 2 * (
+        sigma_i ** (-3)
+    ) + 1 / E / math.cos(beta_i)
+    dB_dsigma_i = -2 * (lambda_i * l_i) ** 2 / 8 * (sigma_i ** (-3))
+    _t = -l_i / (math.cos(beta_i) ** 2) * (dA_dsigma_i * B - A * dB_dsigma_i) / (B ** 2)
+    return _t
+
+
+def manual_diff_sigma_i1_d_l_i(
+    h_i, l_i, h_i1, l_i1, Gi, A, lambda_i, lambda_i1, sigma_i, stringlen, b_i
+):
+    beta_i = math.atan(h_i / l_i)
+    beta_i1 = math.atan(h_i1 / l_i1)
+    A = (
+        Gi / 2 / A
+        + lambda_i * l_i / 2 / math.cos(beta_i)
+        + lambda_i1 * l_i1 / 2 / math.cos(beta_i1)
+        + sigma_i * h_i / l_i
+        + sigma_i / b_i * ((stringlen ** 2 - b_i ** 2) ** 0.5)
+    )
+    B = ((stringlen ** 2 - b_i ** 2) ** 0.5) / b_i + h_i1 / l_i1
+    dA_d_delta_L1 = (
+        sigma_i
+        * (
+            -((stringlen ** 2 - b_i ** 2) ** -0.5) * (b_i ** 2)
+            - (stringlen ** 2 - b_i ** 2) ** 0.5
+        )
+        / (b_i ** 2)
+    )
+    dB_d_delta_L1 = (
+        1
+        / (b_i ** 2)
+        * (
+            -((stringlen ** 2 - b_i ** 2) ** -0.5) * (b_i ** 2)
+            - (stringlen ** 2 - b_i ** 2) ** 0.5
+        )
+    )
+    _t = -(dA_d_delta_L1 * B - A * dB_d_delta_L1) / (B ** 2)
+    return _t
+
+
+def get_lambdify_d_delta_l_i_sigma_i(_d_delta_l_i_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
+        """
+    )
+    return sympy.lambdify(
+        [
+            delta_l_i,
+            l_i,
+            lambda_i,
+            alpha,
+            E,
+            t_e,
+            t_i,
+            lambda_m,
+            t_m,
+            sigma_m,
+            _sigma_i,
+            beta_i,
+        ],
+        _d_delta_l_i_sigma_i,
+    )
+
+
+def get_lambdify_d_sigma_i1_d_l_i(_d_sigma_i1_d_l_i, delta_Li):
+    (
+        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,
+        """
+    )
+
+    return sympy.lambdify(
+        [
+            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,
+            *delta_Li,
+        ],
+        _d_sigma_i1_d_l_i,
+    )

main.py

@@ -2,37 +2,277 @@
 # 新版输电线路大手册 P328
 
 import math
-
+import data
-# 架线时的状态
+import numpy as np
-t_m = 100  # 导线架设时的气温。单位°C
-sigma_m = 10  # 架线时，初伸长未释放前的最低点水平应力。单位N/mm2
-lambda_m = 10  # 导线比载 N/(m.mm)
-t_e = 25  # 架线时考虑初伸长的降温，取正值。单位°C
-
 
 # h_i 悬点高差
 # l_i 悬点档距
-# alpha 导线膨胀系数  1/°C
+# _alpha 导线膨胀系数  1/°C
-# E 弹性系数 N/mm2
+# _elastic 弹性系数 N/mm2
+# _t_e 架线时考虑初伸长的降温，取正值。单位°C
 # lambda_i 计算不平衡张力时导线比载 N/(m.mm)
 # sigma_i 计算不平衡张力时最低点水平应力 单位N/mm2
 # t_i 计算不平衡张力时导线温度 单位°C
-# lambda_m 导线架线时时导线比载 N/(m.mm)
+# _lambda_m 导线架线时时导线比载 N/(m.mm)
-# sigma_m 导线架线时时最低点水平应力 单位N/mm2
+# _sigma_m 导线架线时时最低点水平应力 单位N/mm2
-# t_m 导线架线时时导线温度 单位°C
+# _t_m 导线架线时时导线温度 单位°C
-def delta_li(h_i, l_i, lambda_i, alpha, E, t_i, sigma_i, lambda_m, t_m, sigma_m):
+
-    beta_i = math.atan(h_i / l_i)
+
-    t = (
+def delta_li(
+    h_i: float,
+    l_i,
+    lambda_i: float,
+    _alpha: float,
+    _elastic: float,
+    _t_e: float,
+    t_i: float,
+    sigma_i,
+    _lambda_m: float,
+    _t_m: float,
+    _sigma_m: float,
+) -> float:
+    beta_i = math.atan(h_i / l_i)  # 高差角
+    _delta_li = (
         l_i
         / ((math.cos(beta_i) ** 2) * (1 + (lambda_i * l_i / sigma_i) ** 2 / 8))
         * (
             (l_i * math.cos(beta_i)) ** 2
             / 24
-            * ((lambda_m / sigma_m) ** 2 - (lambda_i / sigma_i) ** 2)
+            * ((_lambda_m / _sigma_m) ** 2 - (lambda_i / sigma_i) ** 2)
-            + ((sigma_i - sigma_m) / E / math.cos(beta_i))
+            + ((sigma_i - _sigma_m) / _elastic / math.cos(beta_i))
-            + alpha * (t_i + t_e - t_m)
+            + _alpha * (t_i + _t_e - _t_m)
         )
     )
+    return _delta_li
 
 
-print("Finished.")
+# area 导线截面 单位mm2
+# sigma_i 第i档内水平应力 单位N/mm2
+# b_i 悬垂串沿线路方向水平偏移距离，沿大号方向为正，反之为负。 单位m
+# length_i 第i基直线塔串长 单位m
+# g_i 第i基直线塔串重 单位N
+# h_i 悬垂串处千中垂位置时，,第 i 基对第 i-1 杆塔上导线悬挂点间的高差大号比小号杆塔悬挂点高者h本身为正值，反之为负值。
+# lambda_i 第i档导线比载 N/(m.mm)
+# h_i1 悬垂串处千中垂位置时，第 i+1 基对第 i 杆塔上导线悬挂点间的高差大号比小号杆塔悬挂点高者h本身为正值，反之为负值。
+# lambda_i1
+def fun_sigma_i1(
+    area: float,
+    sigma_i: float,
+    b_i: float,
+    length_i: float,
+    g_i: float,
+    h_i: float,
+    l_i: float,
+    lambda_i: float,
+    h_i1: float,
+    l_i1: float,
+    lambda_i1: float,
+) -> float:
+    beta_i = math.atan(h_i / l_i)
+    beta_i1 = math.atan(h_i1 / l_i1)
+    _sigma_i1 = (
+        (
+            g_i / 2 / area  # g_i传入时已考虑导线分裂数
+            + lambda_i * l_i / 2 / math.cos(beta_i)
+            + lambda_i1 * l_i1 / 2 / math.cos(beta_i1)
+            + sigma_i * h_i / l_i
+        )
+        + sigma_i / b_i * math.sqrt(length_i ** 2 - b_i ** 2)
+    ) / (math.sqrt(length_i ** 2 - b_i ** 2) / b_i + h_i1 / l_i1)
+    return _sigma_i1
+
+
+# 求解循环。
+def cal_loop():
+    loop_end = data.loop_end  # 最大循环次数
+    # 架线时的状态
+    # 取外过无风
+    string_length = data.string_length  # 串长 单位m
+    string_g = data.string_g  # 串重 单位N
+    t_m = data.t_m  # 导线架设时的气温。单位°C
+    t_e = data.t_e  # 架线时考虑初伸长的降温，取正值。单位°C
+    alpha = data.alpha  # 导线膨胀系数  1/°C
+    elastic = data.elastic  # 弹性系数 N/mm2
+    area = data.area  # 导线面积 mm2
+    lambda_m = data.lambda_m  # 导线比载 N/(m.mm)
+    sigma_m = data.sigma_m  # 架线时，初伸长未释放前的最低点水平应力。单位N/mm2
+    span_count = data.span_count  # 几个档距
+    # n个档距,n-1个直线塔
+    h_array = data.h_array
+    l_array = data.l_array
+    t_i = data.t
+    lambda_i_array = data.lambda_i_array
+    loop_count = 1
+    sigma_0 = sigma_m * 0.2
+    while True:
+        b_i = 0
+        # 一次增加0.1N
+        sigma_0 = sigma_0 + 0.01 / data.area
+        sigma_array = [sigma_0 for _ in range(span_count)]
+        delta_l_i_array = []
+        for i in range(span_count):
+            h_i = h_array[i]
+            l_i = l_array[i]
+            lambda_i = lambda_i_array[i]
+            t_i = t_i
+            sigma_i = sigma_array[i]
+            _delta_l_i = delta_li(
+                h_i,
+                l_i,
+                lambda_i,
+                alpha,
+                elastic,
+                t_e,
+                t_i,
+                sigma_i,
+                lambda_m,
+                t_m,
+                sigma_m,
+            )
+            delta_l_i_array.append(_delta_l_i)
+            b_i += _delta_l_i
+            length_i = string_length
+            g_i = string_g / data.conductor_n
+            if i < span_count - 1:
+                lambda_i1 = lambda_i_array[i + 1]
+                h_i1 = h_array[i + 1]
+                l_i1 = l_array[i + 1]
+                try:
+                    sigma_i1 = fun_sigma_i1(
+                        area,
+                        sigma_i,
+                        b_i,
+                        length_i,
+                        g_i,
+                        h_i,
+                        l_i,
+                        lambda_i,
+                        h_i1,
+                        l_i1,
+                        lambda_i1,
+                    )
+                except ValueError:
+                    break
+                    pass
+                sigma_array[i + 1] = sigma_i1
+        if math.fabs(b_i) < data.epsilon:
+            print("迭代{loop_count}次找到解。".format(loop_count=loop_count))
+            print("悬垂串偏移累加bi为{b_i}".format(b_i=b_i))
+            for i in range(span_count):
+                print("第{i}档导线应力为{tension}".format(i=i, tension=sigma_array[i]))
+            for i in range(span_count - 1):
+                print(
+                    "第{i}串偏移值为{bias}".format(i=i, bias=math.fsum(delta_l_i_array[0:i]))
+                )
+            verify(
+                area,
+                h_array,
+                l_array,
+                string_length,
+                string_g / data.conductor_n,
+                sigma_array,
+                delta_l_i_array,
+                lambda_i_array,
+                t_i,
+                alpha,
+                elastic,
+                t_e,
+                lambda_m,
+                t_m,
+                sigma_m,
+            )
+            break
+        loop_count += 1
+        if loop_count >= loop_end:
+            print("!!!未找到解。")
+            print(sigma_array)
+            print(b_i)
+            break
+
+
+# 检验等式。
+def verify(
+    area: float,
+    h_array: [float],
+    l_array: [float],
+    string_length: [float],
+    string_g: [float],
+    sigma_array: [float],
+    delta_l_i_array: [float],
+    lambda_array: [float],
+    t_i: float,
+    alpha: float,
+    elastic: float,
+    t_e: float,
+    lambda_m: float,
+    t_m: float,
+    sigma_m: float,
+    epsilon: float = 1e-4,
+):
+    # 用新版大手册p329页(5-61)第一个公式校验
+    b_i = 0
+    if math.fabs((math.fsum(delta_l_i_array))) > 1e-5:
+        print("偏移累加不等于0")
+        return
+    for i in range(len(delta_l_i_array)):
+        sigma_i = sigma_array[i]
+        _delta_l_i = delta_l_i_array[i]
+        t_i = t_i
+        # 此处用新版大手册p329页(5-58)校验偏移值。
+        lambda_i = lambda_array[i]
+        h_i = h_array[i]
+        l_i = l_array[i]
+        cal_delta_l_i = delta_li(
+            h_i,
+            l_i,
+            lambda_i,
+            alpha,
+            elastic,
+            t_e,
+            t_i,
+            sigma_i,
+            lambda_m,
+            t_m,
+            sigma_m,
+        )
+        if math.fabs(cal_delta_l_i - _delta_l_i) > 1e-4:
+            print("!!!偏移等式不满足。")
+        if i < len(delta_l_i_array) - 1:
+            sigma_i1 = sigma_array[i + 1]
+            left_equ = sigma_array[i + 1]
+            lambda_i1 = lambda_array[i + 1]
+            h_i1 = h_array[i + 1]
+            l_i1 = l_array[i + 1]
+            beta_i = math.atan(h_i / l_i)
+            beta_i1 = math.atan(h_i1 / l_i1)
+            w_i = (
+                lambda_i * l_i / 2 / math.cos(beta_i)
+                + sigma_i * h_i / l_i
+                + (lambda_i1 * l_i1 / 2 / math.cos(beta_i1) - sigma_i1 * h_i1 / l_i1)
+            )
+            b_i += delta_l_i_array[i]
+            # 新版大手册p329 (5-61) 最上方公式
+            right_equ = sigma_i + b_i / math.sqrt(string_length ** 2 - b_i ** 2) * (
+                string_g / 2 / area + w_i  # string_g已在传入时考虑了导线分裂数
+            )
+            # TODO 等式允许误差是否可调？
+            if math.fabs(right_equ - left_equ) > epsilon:
+                print(math.fabs(right_equ - left_equ))
+                print("!!!应力等式不满足")
+                return
+    print("等式满足。")
+    print("sigma")
+    print(np.array(sigma_array) * area)
+    print("delta_li")
+    print(delta_l_i_array)
+    print("串偏移")
+    b_i = 0
+    b_i_list = []
+    for l_i in delta_l_i_array[:-1]:
+        b_i += l_i
+        b_i_list.append(b_i)
+    print(b_i_list)
+
+
+if __name__ == "__main__":
+    cal_loop()