diff --git a/auto_differentiation.py b/auto_differentiation.py index 979a406..2d28603 100644 --- a/auto_differentiation.py +++ b/auto_differentiation.py @@ -60,25 +60,9 @@ 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): +def get_lambdify_d_delta_l_i_sigma_i(_d_delta_l_i_sigma_i): ( delta_l_i, l_i, @@ -108,46 +92,37 @@ def evaluate_d_delta_l_i_sigma_i(val_delta_l_li, val_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_data), - (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])), - ] - ) - 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 + return sympy.lambdify( + ( + delta_l_i, + l_i, + lambda_i, + alpha, + E, + t_e, + t_i, + lambda_m, + t_m, + sigma_m, + _sigma_i, + beta_i, + ), + _d_delta_l_i_sigma_i, + ) -# C为dσi1dΔli -# C只有n-1行 -def evaluate_d_sigma_i1_d_delta_l_i(val_delta_l_li, val_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) + + +def get_lambdify_d_sigma_i1_d_l_i(_d_sigma_i1_d_l_i): ( G_i, A, @@ -179,6 +154,84 @@ def evaluate_d_sigma_i1_d_delta_l_i(val_delta_l_li, val_sigma_i): beta_i1, """ ) + + return sympy.lambdify( + ( + 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, + *delta_Li, + ), + _d_sigma_i1_d_l_i, + ) + + +fx_d_sigma_i1_d_l_i = get_lambdify_d_sigma_i1_d_l_i(d_sigma_i1_d_l_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): + + 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 = [] @@ -186,36 +239,31 @@ def evaluate_d_sigma_i1_d_delta_l_i(val_delta_l_li, val_sigma_i): if i < j: col.append(0) else: - _val = d_sigma_i1_d_l_i.subs( - [ - (G_i, string_g / conductor_n), - (A, area_data), - (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 + _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], @@ -267,11 +315,11 @@ def evaluate_d_sigma_i1_d_delta_sigma_i(val_delta_li): def solve(): # 初始化 - val_delta_li = [0.1 for _ in range(span_count)] + val_delta_li = [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] - # val_sigma_i = [sigma_m_data for _ in range(span_count)] + # 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: @@ -360,10 +408,10 @@ def solve(): b = sympy.Matrix(b_list) # sympy.pprint(b) x = sympy.linsolve((-A, b)) - x_list = list(x)[0] + x_list = [float(_x) for _x in list(x)[0]] abs_min = [math.fabs(_x) for _x in x_list] abs_min.sort() - if abs_min[-1] < 1e-6: + if abs_min[-1] < 1e-5: break print("最大偏差{max_dx}".format(max_dx=abs_min[-1])) # 更新变量 @@ -373,9 +421,9 @@ def solve(): if loop >= loop_end: print("不收敛") else: - print(loop) - print(val_delta_li) - print(val_sigma_i) + 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) diff --git a/data.py b/data.py index cd37172..414ec8b 100644 --- a/data.py +++ b/data.py @@ -24,7 +24,7 @@ lambda_i_array = [ ] # 取400m代表档距下 sigma_m = 28517 / area # 架线时，初伸长未释放前的最低点水平应力。单位N/mm2 -span_count = 9 # 几个档距 +span_count = 10 # 几个档距 # n个档距,n-1个直线塔 h_array = [ 0, diff --git a/main.py b/main.py index 16aa0c4..c4c3632 100644 --- a/main.py +++ b/main.py @@ -132,10 +132,10 @@ def cal_loop(): b_i += _delta_l_i length_i = string_length g_i = string_g / data.conductor_n - h_i1 = h_array[i + 1] - l_i1 = l_array[i + 1] if i < span_count - 1: lambda_i1 = lambda_i_array[i + 1] + h_i1 = h_array[i + 1] + l_i1 = l_array[i + 1] try: sigma_i1 = fun_sigma_i1( area, @@ -257,10 +257,10 @@ def verify( print("!!!应力等式不满足") return print("等式满足。") - print('sigma') - print(sigma_array) - print('delta_li') - print(delta_l_i_array) + # print('sigma') + # print(sigma_array) + # print('delta_li') + # print(delta_l_i_array) if __name__ == "__main__":