自动微分收敛了。

auto_differentiation.py

@@ -0,0 +1,352 @@
+# 利用自动微分计算
+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.elastic  # 弹性系数 N/mm2
+area = 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),
+                (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),
+                        (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,
+                    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,
+                        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)
+
+
+solve()
+print("Finished.")

data.py

@@ -1,3 +1,4 @@
+#loop_end = 10000000  # 最大循环次数
 loop_end = 1000000  # 最大循环次数
 # 架线时的状态
 # 取外过无风
@@ -9,11 +10,13 @@ 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.3,lambda_m,lambda_m,lambda_m]
 # 取400m代表档距下
 sigma_m = 28517 / area  # 架线时，初伸长未释放前的最低点水平应力。单位N/mm2
-span_count = 4  # 几个档距
+span_count = 3 # 几个档距
 # n个档距,n-1个直线塔
-h_array = [0, 0, 0, 0]
-l_array = [200, 300, 1400, 500]
-t_array = [15, 15, 15, 15]
-epslon = 1e-3  # 收敛判据
+h_array = [0, 0, 0, 0, 0]
+l_array = [400, 300, 300, 500, 300]
+t = 15
+epsilon = 1e-4  # 收敛判据
+conductor_n = 6  # 导线分裂数

exp.py

@@ -0,0 +1,118 @@
+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,"""
+    )
+    # beta_i = sympy.atan(h_i / l_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)
+        )
+    )
+    # d_delta_li_sigma_i = sympy.diff(_delta_li, sigma_i)
+    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

main.py

@@ -19,13 +19,13 @@ import data
 
 def delta_li(
     h_i: float,
-    l_i: float,
+    l_i,
     lambda_i: float,
     _alpha: float,
     _elastic: float,
     _t_e: float,
     t_i: float,
-    sigma_i: float,
+    sigma_i,
     _lambda_m: float,
     _t_m: float,
     _sigma_m: float,
@@ -71,7 +71,7 @@ def fun_sigma_i1(
     beta_i1 = math.atan(h_i1 / l_i1)
     _sigma_i1 = (
         (
-            g_i / 2 / area #g_i传入时已考虑导线分裂数
+            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
@@ -99,21 +99,21 @@ def cal_loop():
     # n个档距,n-1个直线塔
     h_array = data.h_array
     l_array = data.l_array
-    t_array = data.t_array
-    lambda_array = [lambda_m for j in range(span_count)]
+    t_i = data.t
+    lambda_i_array = data.lambda_i_array
     loop_count = 1
-    sigma_0 = sigma_m * 0.5
+    sigma_0 = sigma_m * 0.2
     while True:
         b_i = 0
         # 一次增加0.1N
-        sigma_0 = sigma_0 + 0.01 / data.area
+        sigma_0 = sigma_0 + 0.1 / data.area
         sigma_array = [sigma_0 for j in range(span_count)]
         delta_l_i_array = []
-        for i in range(span_count - 1):
+        for i in range(span_count):
             h_i = h_array[i]
             l_i = l_array[i]
-            lambda_i = lambda_array[i]
-            t_i = t_array[i]
+            lambda_i = lambda_i_array[i]
+            t_i = t_i
             sigma_i = sigma_array[i]
             _delta_l_i = delta_li(
                 h_i,
@@ -134,33 +134,33 @@ def cal_loop():
             g_i = string_g / data.conductor_n
             h_i1 = h_array[i + 1]
             l_i1 = l_array[i + 1]
-            lambda_i1 = lambda_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
-        loop_count += 1
+            if i<span_count-1:
+                lambda_i1 = lambda_i_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=delta_l_i_array[i]))
+                print("第{i}串偏移值为{bias}".format(i=i, bias=math.fsum(delta_l_i_array[0:i])))
             verify(
                 area,
                 h_array,
@@ -169,8 +169,8 @@ def cal_loop():
                 string_g / data.conductor_n,
                 sigma_array,
                 delta_l_i_array,
-                lambda_array,
-                t_array,
+                lambda_i_array,
+                t_i,
                 alpha,
                 elastic,
                 t_e,
@@ -179,6 +179,7 @@ def cal_loop():
                 sigma_m,
             )
             break
+        loop_count += 1
         if loop_count >= loop_end:
             print("!!!未找到解。")
             print(sigma_array)
@@ -196,7 +197,7 @@ def verify(
     sigma_array: [float],
     delta_l_i_array: [float],
     lambda_array: [float],
-    t_array,
+    t_i,
     alpha,
     elastic,
     t_e,
@@ -208,17 +209,12 @@ def verify(
     b_i = 0
     for i in range(len(delta_l_i_array)):
         sigma_i = sigma_array[i]
-        sigma_i1 = sigma_array[i + 1]
-        left_equ = sigma_array[i + 1]
         _delta_l_i = delta_l_i_array[i]
-        t_i = t_array[i]
+        t_i = t_i
         # 此处用新版大手册p329页(5-58)校验偏移值。
         lambda_i = lambda_array[i]
-        lambda_i1 = lambda_array[i + 1]
         h_i = h_array[i]
-        h_i1 = h_array[i + 1]
         l_i = l_array[i]
-        l_i1 = l_array[i + 1]
         cal_delta_l_i = delta_li(
             h_i,
             l_i,
@@ -234,26 +230,49 @@ def verify(
         )
         if math.fabs(cal_delta_l_i - _delta_l_i) > 1e-4:
             print("!!!偏移等式不满足。")
-        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]
-        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) > 1e-4:
-            print(math.fabs(right_equ - left_equ))
-            print("!!!应力等式不满足")
-            return
+        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) > 1e-4:
+                print(math.fabs(right_equ - left_equ))
+                print("!!!应力等式不满足")
+                return
     print("等式满足。")
 
+if __name__=='__main__':
+    cal_loop()
+import sympy
 
-cal_loop()
-
-
-print("Finished.")
+# sympy.init_printing(use_unicode=False)
+# a, b, c, d = sympy.symbols("a b c d")
+# m, w, n, t = sympy.symbols("m w n t")
+# b1, b2, b3 = sympy.symbols("b1 b2 b3")
+# x = sympy.symbols("x")
+# fx = (a * sympy.exp(-x / b1) + b * sympy.exp(-t / b1) + c) * sympy.cos(
+#     w * x + d
+# ) + m * sympy.exp(-t / b3) * sympy.sin(w * x + n)
+# int_f = sympy.integrate(fx, x)
+# sympy.pprint(int_f, use_unicode=True)
+# sympy.print_mathml(int_f)
+# import sympy.matrices
+# u_list=sympy.symbols("u1:101")
+# sympy.pprint(u_list)
+# m=sympy.Matrix(u_list)
+# sympy.pprint(m)
+# print("Finished.")