利用numpy增强数值稳定性

auto_differentiation.py

@@ -1,9 +1,12 @@
 # 利用自动微分计算
+from typing import List
+
 import sympy
 import data
 import exp
 import math
 import main
+import numpy as np
 
 sympy.init_printing()
 # h_i 悬点高差
@@ -26,7 +29,6 @@ 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  # 最大循环次数
@@ -44,19 +46,13 @@ sigma_m_data = data.sigma_m  # 架线时，初伸长未释放前的最低点水
 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)
@@ -113,11 +109,7 @@ def get_lambdify_d_delta_l_i_sigma_i(_d_delta_l_i_sigma_i):
 
 fx_d_delta_l_i_sigma_i = get_lambdify_d_delta_l_i_sigma_i(d_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)
 
@@ -315,7 +307,7 @@ def evaluate_d_sigma_i1_d_delta_sigma_i(val_delta_li):
 
 def solve():
     # 初始化
-    val_delta_li = [1 for _ in range(span_count)]
+    val_delta_li = [0.1 for _ in range(span_count)]
     # val_delta_li = [0.15931589580330385, -0.19294219226439696, 0.035236744603670586, -0.03124795962518869, 0.034422847061933395, 0.1639551737321582, -0.18876605469917845, 0.03862307099222713, -0.01858088432390902]
 
     # val_sigma_i = [175.44277576372207, 176.07105062365383, 175.94897602628689, 175.95574563415994, 175.83116207025824, 175.84778616876434, 176.5132398968081, 176.4405408454816, 176.51871673188586]
@@ -346,7 +338,6 @@ def solve():
         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]
@@ -381,35 +372,24 @@ def solve():
                         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)
+        # b = sympy.Matrix(b_list)
-        # sympy.pprint(b)
+        # x = sympy.linsolve((-A, b))
-        x = sympy.linsolve((-A, b))
+        # x_list = [float(_x) for _x in list(x)[0]]
-        x_list = [float(_x) for _x in list(x)[0]]
+        # numpy
-        abs_min = [math.fabs(_x) for _x in x_list]
+        # print(A.tolist())
+        AA = np.array(A.tolist(), np.float64)
+        # AA=np.array([[1, 0, 0, -0.0104519104486432, 0, 0], [0, 1, 0, 0, -0.00782784930594165, 0], [0, 0, 1, 0, 0, -0.00640057303461462], [-3.78731866382623, 0, 0, -1.00000000000000, 1, 0], [-3.47643075688507, -3.47643075688507, 0, 0, -1.00000000000000, 1], [1, 1, 1, 0, 0, 0]])
+        # print(b.transpose().tolist()[0])
+        b = np.array(b_list, dtype=np.float64)
+        # BB=[-0.159573673689794, 0.140383878899664, 0.00728144380483030, -0.378687120179166, -0.694957566239367, 0.300000000000000]
+        x = np.linalg.solve(-AA, b)
+        x_list = x
+        abs_min: List[float] = [math.fabs(_x) for _x in x_list]
         abs_min.sort()
         if abs_min[-1] < 1e-5:
             break

data.py

@@ -12,15 +12,20 @@ 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 * 1.2,
+    lambda_m,
+    lambda_m,
+    lambda_m,
+    lambda_m * 0.9,
+    lambda_m * 1.1,
+    lambda_m,
+    lambda_m,
+    lambda_m,
     lambda_m,
     lambda_m,
     lambda_m,
     lambda_m * 0.9,
     lambda_m * 1.3,
-    lambda_m,
-    lambda_m,
-    lambda_m,
 ]
 # 取400m代表档距下
 sigma_m = 28517 / area  # 架线时，初伸长未释放前的最低点水平应力。单位N/mm2
@@ -37,8 +42,29 @@ h_array = [
     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,
 ]
-l_array = [400, 300, 300, 500, 300, 400, 300, 300, 500, 300]
 t = 15
 epsilon = 1e-4  # 收敛判据
 conductor_n = 6  # 导线分裂数

main.py

@@ -69,6 +69,7 @@ def fun_sigma_i1(
 ) -> float:
     beta_i = math.atan(h_i / l_i)
     beta_i1 = math.atan(h_i1 / l_i1)
+    try:
         _sigma_i1 = (
             (
                 g_i / 2 / area  # g_i传入时已考虑导线分裂数
@@ -78,6 +79,8 @@ def fun_sigma_i1(
             )
             + 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)
+    except:
+        a=12
     return _sigma_i1