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