From 2750267b945d32b0ed973cfed17912deb0c567a5 Mon Sep 17 00:00:00 2001
From: facat <facat@facat.com>
Date: Mon, 14 Dec 2020 10:31:29 +0800
Subject: [PATCH] =?UTF-8?q?1.=E5=A2=9E=E5=8A=A0=E4=B8=8E=E6=89=8B=E5=8A=A8?=
 =?UTF-8?q?=E5=BE=AE=E5=88=86=E5=85=AC=E5=BC=8F=E6=AF=94=E8=BE=83=E3=80=82?=
 =?UTF-8?q?=202.=E7=89=9B=E9=A1=BF=E6=B3=95=E5=AF=B9=E5=88=9D=E5=80=BC?=
 =?UTF-8?q?=E6=AF=94=E8=BE=83=E6=95=8F=E6=84=9F=E3=80=82?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 auto_differentiation.py | 42 ++++++++++++++++++++++++++-----
 data.py                 | 32 +++++++++++++++++++----
 exp.py                  | 56 ++++++++++++++++++++++++++++++++++++++++-
 main.py                 |  4 +++
 4 files changed, 122 insertions(+), 12 deletions(-)

diff --git a/auto_differentiation.py b/auto_differentiation.py
index 0174056..979a406 100644
--- a/auto_differentiation.py
+++ b/auto_differentiation.py
@@ -126,6 +126,21 @@ def evaluate_d_delta_l_i_sigma_i(val_delta_l_li, val_sigma_i):
                 (beta_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
 
@@ -201,6 +216,21 @@ def evaluate_d_sigma_i1_d_delta_l_i(val_delta_l_li, val_sigma_i):
                 for index, li in enumerate(delta_Li):
                     _val = _val.subs(li, _val_delta_l_li[index])
                     pass
+                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)
@@ -237,11 +267,11 @@ def evaluate_d_sigma_i1_d_delta_sigma_i(val_delta_li):
 
 def solve():
     # 初始化
-    val_delta_li = [0.1 for i in range(span_count)]
-    # val_delta_li = [0.15864687475316822, -0.1935189734784845, 0.03478489898855073]
+    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 = [sigma_m_data for _ in range(span_count)]
-    # val_sigma_i = [175.38451579479482, 176.01015153076614, 175.88355419459572]
+    val_sigma_i = [175.44277576372207, 176.07105062365383, 175.94897602628689, 175.95574563415994, 175.83116207025824, 175.84778616876434, 176.5132398968081, 176.4405408454816, 176.51871673188586]
+    # val_sigma_i = [sigma_m_data for _ in range(span_count)]
 
     loop = 0
     while True:
@@ -333,7 +363,7 @@ def solve():
         x_list = list(x)[0]
         abs_min = [math.fabs(_x) for _x in x_list]
         abs_min.sort()
-        if abs_min[-1] < 1e-5:
+        if abs_min[-1] < 1e-6:
             break
         print("最大偏差{max_dx}".format(max_dx=abs_min[-1]))
         # 更新变量
@@ -346,7 +376,7 @@ def solve():
         print(loop)
         print(val_delta_li)
         print(val_sigma_i)
-        verify(val_delta_li,val_sigma_i)
+        verify(val_delta_li, val_sigma_i)
 
 
 def verify(val_delta_li, val_sigma_i):
diff --git a/data.py b/data.py
index 45c8bb9..cd37172 100644
--- a/data.py
+++ b/data.py
@@ -1,4 +1,4 @@
-#loop_end = 10000000  # 最大循环次数
+# loop_end = 10000000  # 最大循环次数
 loop_end = 1000000  # 最大循环次数
 # 架线时的状态
 # 取外过无风
@@ -10,13 +10,35 @@ 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]
+lambda_i_array = [
+    lambda_m * 0.9,
+    lambda_m * 1.3,
+    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
-span_count = 3 # 几个档距
+span_count = 9  # 几个档距
 # n个档距,n-1个直线塔
-h_array = [0, 0, 0, 0, 0]
-l_array = [400, 300, 300, 500, 300]
+h_array = [
+    0,
+    0,
+    0,
+    0,
+    0,
+    0,
+    0,
+    0,
+    0,
+    0,
+]
+l_array = [400, 300, 300, 500, 300, 400, 300, 300, 500, 300]
 t = 15
 epsilon = 1e-4  # 收敛判据
 conductor_n = 6  # 导线分裂数
diff --git a/exp.py b/exp.py
index 83f7cb6..f967c4d 100644
--- a/exp.py
+++ b/exp.py
@@ -27,7 +27,7 @@ def delta_li():
         t_m,
         sigma_m,
         sigma_i,
-        beta_i
+        beta_i,
     ) = sympy.symbols(
         """ 
         delta_l_i,
@@ -106,6 +106,7 @@ def fun_sigma_i1(delta_Li):
         for f in delta_Li:
             _t += f
         return _t
+
     _sigma_i1 = sigma_i1 - (
         (
             G_i / 2 / A  # G_i传入时已考虑导线分裂数
@@ -116,3 +117,56 @@ def fun_sigma_i1(delta_Li):
         + 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
diff --git a/main.py b/main.py
index c3c3985..16aa0c4 100644
--- a/main.py
+++ b/main.py
@@ -257,6 +257,10 @@ def verify(
                 print("!!!应力等式不满足")
                 return
     print("等式满足。")
+    print('sigma')
+    print(sigma_array)
+    print('delta_li')
+    print(delta_l_i_array)
 
 
 if __name__ == "__main__":