320 lines
10 KiB
Python
320 lines
10 KiB
Python
# 利用自动微分计算
|
||
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.")
|