diff --git a/sim_lab/halley_comet_gradient_descent.py b/sim_lab/halley_comet_gradient_descent.py new file mode 100644 index 0000000000000000000000000000000000000000..58a854db5dca250211093a352d3388599cbad291 --- /dev/null +++ b/sim_lab/halley_comet_gradient_descent.py @@ -0,0 +1,81 @@ +import numpy as np + +from common.consts import AU +from sim_scenes.solar_system.halley_comet_lib import HalleyCometParams +from sim_lab.halley_comet_research_calc import calc_simulator + + +def target_function(x, y, z): + params = HalleyCometParams( + start_time='1982-09-24 00:00:00', + # init_velocity=[-2.836 + (x / 1000), 4.705 + (y / 1000), 8.85 + (z / 1000)], + init_velocity=[x, y, z], + init_position=[0, -5 * AU, -10 * AU] + ) + return calc_simulator(params) + + +def loss_function(a, b, c): + # 0.586AU, 35.1AU, 13817天 + target_a = 0.586 + target_b = 35.1 + target_c = 13817.0 + + loss_a = (a - target_a) ** 2 + loss_b = (b - target_b) ** 2 + loss_c = (c - target_c) ** 2 + + total_loss = loss_a + loss_b + loss_c + + return total_loss + + +def compute_gradients(x, y, z, epsilon=1e-5): + # 数值方法计算梯度 + grad_x = (loss_function(*target_function(x + epsilon, y, z)) - loss_function(*target_function(x, y, z))) / epsilon + grad_y = (loss_function(*target_function(x, y + epsilon, z)) - loss_function(*target_function(x, y, z))) / epsilon + grad_z = (loss_function(*target_function(x, y, z + epsilon)) - loss_function(*target_function(x, y, z))) / epsilon + + return grad_x, grad_y, grad_z + + +def adaptive_learning_rate(iteration): + # 自适应学习率调整 + return 0.01 / np.sqrt(iteration + 1) + + +def gradient_descent(init_velocity, iterations): + x, y, z = init_velocity + + for i in range(1, iterations + 1): # 从1开始，避免除以0 + a, b, c = target_function(x, y, z) + loss = loss_function(a, b, c) + + # 计算梯度 + grad_x, grad_y, grad_z = compute_gradients(x, y, z) + + # 自适应学习率调整 + learning_rate = adaptive_learning_rate(i) + + # 更新参数 + x = x - learning_rate * grad_x + y = y - learning_rate * grad_y + z = z - learning_rate * grad_z + + # 打印当前迭代的结果 + if i % 100 == 0: + print(f"Iteration {i}: Loss = {loss}, a = {a}, b = {b}, c = {c}") + + return x, y, z + + +# 设置初始值 +init_velocity = [-2.816, 4.725, 8.87] # 请替换为实际的初始值 +iterations = 1000 + +# 运行梯度下降算法 +final_x, final_y, final_z = gradient_descent(init_velocity, iterations) + +# 打印最终结果 +final_a, final_b, final_c = target_function(final_x, final_y, final_z) +print(f"Final Result: a = {final_a}, b = {final_b}, c = {final_c}") diff --git a/sim_lab/halley_comet_log.txt b/sim_lab/halley_comet_log.txt new file mode 100644 index 0000000000000000000000000000000000000000..d4025475786e6b5f68dac594be82d008d9e01642 --- /dev/null +++ b/sim_lab/halley_comet_log.txt @@ -0,0 +1,505 @@ +Index:1 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +WARNING: TimeDeltaMissingUnitWarning: Numerical value without unit or explicit format passed to TimeDelta, assuming days [astropy.time.core] + + +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.856, 4.6850000000000005, 8.83], init_position=[0, -747989353.5, -1495978707.0]) +远日点:34.775 AU, 近日点:0.591 AU, 天数: 13542 +Index:2 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.856, 4.6850000000000005, 8.84], init_position=[0, -747989353.5, -1495978707.0]) +远日点:34.900 AU, 近日点:0.591 AU, 天数: 13609 +Index:3 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.856, 4.6850000000000005, 8.85], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.025 AU, 近日点:0.592 AU, 天数: 13682 +Index:4 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.856, 4.6850000000000005, 8.86], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.152 AU, 近日点:0.593 AU, 天数: 13754 +Index:5 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.856, 4.6850000000000005, 8.87], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.279 AU, 近日点:0.594 AU, 天数: 13827 +Index:6 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.856, 4.695, 8.83], init_position=[0, -747989353.5, -1495978707.0]) +远日点:34.840 AU, 近日点:0.591 AU, 天数: 13578 +Index:7 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.856, 4.695, 8.84], init_position=[0, -747989353.5, -1495978707.0]) +远日点:34.966 AU, 近日点:0.592 AU, 天数: 13645 +Index:8 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.856, 4.695, 8.85], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.091 AU, 近日点:0.592 AU, 天数: 13718 +Index:9 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.856, 4.695, 8.86], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.218 AU, 近日点:0.593 AU, 天数: 13791 +Index:10 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.856, 4.695, 8.87], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.346 AU, 近日点:0.594 AU, 天数: 13876 +Index:11 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.856, 4.705, 8.83], init_position=[0, -747989353.5, -1495978707.0]) +远日点:34.907 AU, 近日点:0.592 AU, 天数: 13615 +Index:12 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.856, 4.705, 8.84], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.032 AU, 近日点:0.592 AU, 天数: 13688 +Index:13 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.856, 4.705, 8.85], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.158 AU, 近日点:0.593 AU, 天数: 13761 +Index:14 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.856, 4.705, 8.86], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.285 AU, 近日点:0.594 AU, 天数: 13834 +Index:15 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.856, 4.705, 8.87], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.414 AU, 近日点:0.593 AU, 天数: 13913 +Index:16 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.856, 4.715, 8.83], init_position=[0, -747989353.5, -1495978707.0]) +远日点:34.973 AU, 近日点:0.592 AU, 天数: 13651 +Index:17 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.856, 4.715, 8.84], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.099 AU, 近日点:0.593 AU, 天数: 13724 +Index:18 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.856, 4.715, 8.85], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.225 AU, 近日点:0.594 AU, 天数: 13797 +Index:19 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.856, 4.715, 8.86], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.353 AU, 近日点:0.595 AU, 天数: 13876 +Index:20 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.856, 4.715, 8.87], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.482 AU, 近日点:0.593 AU, 天数: 13955 +Index:21 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.856, 4.725, 8.83], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.040 AU, 近日点:0.593 AU, 天数: 13694 +Index:22 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.856, 4.725, 8.84], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.165 AU, 近日点:0.594 AU, 天数: 13761 +Index:23 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.856, 4.725, 8.85], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.293 AU, 近日点:0.595 AU, 天数: 13840 +Index:24 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.856, 4.725, 8.86], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.421 AU, 近日点:0.594 AU, 天数: 13919 +Index:25 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.856, 4.725, 8.87], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.551 AU, 近日点:0.593 AU, 天数: 13998 +Index:26 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.8459999999999996, 4.6850000000000005, 8.83], init_position=[0, -747989353.5, -1495978707.0]) +远日点:34.740 AU, 近日点:0.587 AU, 天数: 13517 +Index:27 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.8459999999999996, 4.6850000000000005, 8.84], init_position=[0, -747989353.5, -1495978707.0]) +远日点:34.864 AU, 近日点:0.588 AU, 天数: 13584 +Index:28 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.8459999999999996, 4.6850000000000005, 8.85], init_position=[0, -747989353.5, -1495978707.0]) +远日点:34.990 AU, 近日点:0.588 AU, 天数: 13657 +Index:29 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.8459999999999996, 4.6850000000000005, 8.86], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.115 AU, 近日点:0.590 AU, 天数: 13730 +Index:30 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.8459999999999996, 4.6850000000000005, 8.87], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.243 AU, 近日点:0.589 AU, 天数: 13809 +Index:31 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.8459999999999996, 4.695, 8.83], init_position=[0, -747989353.5, -1495978707.0]) +远日点:34.805 AU, 近日点:0.588 AU, 天数: 13554 +Index:32 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.8459999999999996, 4.695, 8.84], init_position=[0, -747989353.5, -1495978707.0]) +远日点:34.930 AU, 近日点:0.588 AU, 天数: 13627 +Index:33 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.8459999999999996, 4.695, 8.85], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.055 AU, 近日点:0.589 AU, 天数: 13694 +Index:34 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.8459999999999996, 4.695, 8.86], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.182 AU, 近日点:0.590 AU, 天数: 13773 +Index:35 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.8459999999999996, 4.695, 8.87], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.309 AU, 近日点:0.589 AU, 天数: 13852 +Index:36 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.8459999999999996, 4.705, 8.83], init_position=[0, -747989353.5, -1495978707.0]) +远日点:34.871 AU, 近日点:0.588 AU, 天数: 13590 +Index:37 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.8459999999999996, 4.705, 8.84], init_position=[0, -747989353.5, -1495978707.0]) +远日点:34.996 AU, 近日点:0.589 AU, 天数: 13663 +Index:38 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.8459999999999996, 4.705, 8.85], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.122 AU, 近日点:0.590 AU, 天数: 13736 +Index:39 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.8459999999999996, 4.705, 8.86], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.249 AU, 近日点:0.590 AU, 天数: 13815 +Index:40 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.8459999999999996, 4.705, 8.87], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.378 AU, 近日点:0.588 AU, 天数: 13888 +Index:41 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.8459999999999996, 4.715, 8.83], init_position=[0, -747989353.5, -1495978707.0]) +远日点:34.937 AU, 近日点:0.589 AU, 天数: 13627 +Index:42 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.8459999999999996, 4.715, 8.84], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.063 AU, 近日点:0.590 AU, 天数: 13700 +Index:43 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.8459999999999996, 4.715, 8.85], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.189 AU, 近日点:0.591 AU, 天数: 13773 +Index:44 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.8459999999999996, 4.715, 8.86], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.317 AU, 近日点:0.590 AU, 天数: 13858 +Index:45 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.8459999999999996, 4.715, 8.87], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.445 AU, 近日点:0.588 AU, 天数: 13931 +Index:46 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.8459999999999996, 4.725, 8.83], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.004 AU, 近日点:0.590 AU, 天数: 13669 +Index:47 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.8459999999999996, 4.725, 8.84], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.129 AU, 近日点:0.591 AU, 天数: 13742 +Index:48 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.8459999999999996, 4.725, 8.85], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.256 AU, 近日点:0.591 AU, 天数: 13821 +Index:49 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.8459999999999996, 4.725, 8.86], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.384 AU, 近日点:0.589 AU, 天数: 13894 +Index:50 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.8459999999999996, 4.725, 8.87], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.514 AU, 近日点:0.588 AU, 天数: 13967 +Index:51 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.836, 4.6850000000000005, 8.83], init_position=[0, -747989353.5, -1495978707.0]) +远日点:34.705 AU, 近日点:0.584 AU, 天数: 13493 +Index:52 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.836, 4.6850000000000005, 8.84], init_position=[0, -747989353.5, -1495978707.0]) +远日点:34.829 AU, 近日点:0.584 AU, 天数: 13566 +Index:53 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.836, 4.6850000000000005, 8.85], init_position=[0, -747989353.5, -1495978707.0]) +远日点:34.953 AU, 近日点:0.585 AU, 天数: 13633 +Index:54 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.836, 4.6850000000000005, 8.86], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.079 AU, 近日点:0.585 AU, 天数: 13712 +Index:55 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.836, 4.6850000000000005, 8.87], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.206 AU, 近日点:0.584 AU, 天数: 13785 +Index:56 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.836, 4.695, 8.83], init_position=[0, -747989353.5, -1495978707.0]) +远日点:34.770 AU, 近日点:0.584 AU, 天数: 13530 +Index:57 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.836, 4.695, 8.84], init_position=[0, -747989353.5, -1495978707.0]) +远日点:34.894 AU, 近日点:0.585 AU, 天数: 13602 +Index:58 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.836, 4.695, 8.85], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.020 AU, 近日点:0.586 AU, 天数: 13675 +Index:59 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.836, 4.695, 8.86], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.146 AU, 近日点:0.585 AU, 天数: 13748 +Index:60 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.836, 4.695, 8.87], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.274 AU, 近日点:0.584 AU, 天数: 13827 +Index:61 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.836, 4.705, 8.83], init_position=[0, -747989353.5, -1495978707.0]) +远日点:34.836 AU, 近日点:0.585 AU, 天数: 13566 +Index:62 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.836, 4.705, 8.84], init_position=[0, -747989353.5, -1495978707.0]) +远日点:34.960 AU, 近日点:0.586 AU, 天数: 13639 +Index:63 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.836, 4.705, 8.85], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.086 AU, 近日点:0.586 AU, 天数: 13718 +Index:64 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.836, 4.705, 8.86], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.213 AU, 近日点:0.585 AU, 天数: 13791 +Index:65 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.836, 4.705, 8.87], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.341 AU, 近日点:0.584 AU, 天数: 13864 +Index:66 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.836, 4.715, 8.83], init_position=[0, -747989353.5, -1495978707.0]) +远日点:34.901 AU, 近日点:0.586 AU, 天数: 13609 +Index:67 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.836, 4.715, 8.84], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.027 AU, 近日点:0.587 AU, 天数: 13675 +Index:68 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.836, 4.715, 8.85], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.153 AU, 近日点:0.586 AU, 天数: 13754 +Index:69 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.836, 4.715, 8.86], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.281 AU, 近日点:0.585 AU, 天数: 13827 +Index:70 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.836, 4.715, 8.87], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.409 AU, 近日点:0.584 AU, 天数: 13907 +Index:71 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.836, 4.725, 8.83], init_position=[0, -747989353.5, -1495978707.0]) +远日点:34.969 AU, 近日点:0.587 AU, 天数: 13645 +Index:72 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.836, 4.725, 8.84], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.094 AU, 近日点:0.588 AU, 天数: 13724 +Index:73 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.836, 4.725, 8.85], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.221 AU, 近日点:0.586 AU, 天数: 13797 +Index:74 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.836, 4.725, 8.86], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.348 AU, 近日点:0.585 AU, 天数: 13870 +Index:75 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.836, 4.725, 8.87], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.478 AU, 近日点:0.584 AU, 天数: 13949 +Index:76 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.826, 4.6850000000000005, 8.83], init_position=[0, -747989353.5, -1495978707.0]) +远日点:34.670 AU, 近日点:0.580 AU, 天数: 13469 +Index:77 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.826, 4.6850000000000005, 8.84], init_position=[0, -747989353.5, -1495978707.0]) +远日点:34.794 AU, 近日点:0.581 AU, 天数: 13542 +Index:78 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.826, 4.6850000000000005, 8.85], init_position=[0, -747989353.5, -1495978707.0]) +远日点:34.918 AU, 近日点:0.582 AU, 天数: 13615 +Index:79 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.826, 4.6850000000000005, 8.86], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.044 AU, 近日点:0.580 AU, 天数: 13688 +Index:80 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.826, 4.6850000000000005, 8.87], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.171 AU, 近日点:0.579 AU, 天数: 13761 +Index:81 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.826, 4.695, 8.83], init_position=[0, -747989353.5, -1495978707.0]) +远日点:34.735 AU, 近日点:0.581 AU, 天数: 13505 +Index:82 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.826, 4.695, 8.84], init_position=[0, -747989353.5, -1495978707.0]) +远日点:34.859 AU, 近日点:0.582 AU, 天数: 13578 +Index:83 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.826, 4.695, 8.85], init_position=[0, -747989353.5, -1495978707.0]) +远日点:34.984 AU, 近日点:0.582 AU, 天数: 13657 +Index:84 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.826, 4.695, 8.86], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.110 AU, 近日点:0.580 AU, 天数: 13724 +Index:85 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.826, 4.695, 8.87], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.237 AU, 近日点:0.579 AU, 天数: 13803 +Index:86 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.826, 4.705, 8.83], init_position=[0, -747989353.5, -1495978707.0]) +远日点:34.801 AU, 近日点:0.582 AU, 天数: 13548 +Index:87 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.826, 4.705, 8.84], init_position=[0, -747989353.5, -1495978707.0]) +远日点:34.924 AU, 近日点:0.583 AU, 天数: 13615 +Index:88 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.826, 4.705, 8.85], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.051 AU, 近日点:0.581 AU, 天数: 13694 +Index:89 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.826, 4.705, 8.86], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.177 AU, 近日点:0.580 AU, 天数: 13767 +Index:90 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.826, 4.705, 8.87], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.305 AU, 近日点:0.579 AU, 天数: 13840 +Index:91 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.826, 4.715, 8.83], init_position=[0, -747989353.5, -1495978707.0]) +远日点:34.867 AU, 近日点:0.583 AU, 天数: 13584 +Index:92 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.826, 4.715, 8.84], init_position=[0, -747989353.5, -1495978707.0]) +远日点:34.992 AU, 近日点:0.583 AU, 天数: 13657 +Index:93 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.826, 4.715, 8.85], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.117 AU, 近日点:0.581 AU, 天数: 13730 +Index:94 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.826, 4.715, 8.86], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.244 AU, 近日点:0.580 AU, 天数: 13803 +Index:95 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.826, 4.715, 8.87], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.373 AU, 近日点:0.579 AU, 天数: 13882 +Index:96 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.826, 4.725, 8.83], init_position=[0, -747989353.5, -1495978707.0]) +远日点:34.933 AU, 近日点:0.584 AU, 天数: 13621 +Index:97 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.826, 4.725, 8.84], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.058 AU, 近日点:0.582 AU, 天数: 13700 +Index:98 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.826, 4.725, 8.85], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.184 AU, 近日点:0.581 AU, 天数: 13773 +Index:99 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.826, 4.725, 8.86], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.313 AU, 近日点:0.580 AU, 天数: 13846 +Index:100 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.826, 4.725, 8.87], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.440 AU, 近日点:0.579 AU, 天数: 13919 +Index:101 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.816, 4.6850000000000005, 8.83], init_position=[0, -747989353.5, -1495978707.0]) +远日点:34.636 AU, 近日点:0.577 AU, 天数: 13450 +Index:102 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.816, 4.6850000000000005, 8.84], init_position=[0, -747989353.5, -1495978707.0]) +远日点:34.759 AU, 近日点:0.578 AU, 天数: 13517 +Index:103 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.816, 4.6850000000000005, 8.85], init_position=[0, -747989353.5, -1495978707.0]) +远日点:34.883 AU, 近日点:0.577 AU, 天数: 13596 +Index:104 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.816, 4.6850000000000005, 8.86], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.009 AU, 近日点:0.576 AU, 天数: 13663 +Index:105 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.816, 4.6850000000000005, 8.87], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.134 AU, 近日点:0.575 AU, 天数: 13736 +Index:106 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.816, 4.695, 8.83], init_position=[0, -747989353.5, -1495978707.0]) +远日点:34.701 AU, 近日点:0.578 AU, 天数: 13487 +Index:107 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.816, 4.695, 8.84], init_position=[0, -747989353.5, -1495978707.0]) +远日点:34.824 AU, 近日点:0.578 AU, 天数: 13560 +Index:108 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.816, 4.695, 8.85], init_position=[0, -747989353.5, -1495978707.0]) +远日点:34.949 AU, 近日点:0.577 AU, 天数: 13633 +Index:109 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.816, 4.695, 8.86], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.074 AU, 近日点:0.575 AU, 天数: 13706 +Index:110 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.816, 4.695, 8.87], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.202 AU, 近日点:0.575 AU, 天数: 13779 +Index:111 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.816, 4.705, 8.83], init_position=[0, -747989353.5, -1495978707.0]) +远日点:34.765 AU, 近日点:0.579 AU, 天数: 13523 +Index:112 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.816, 4.705, 8.84], init_position=[0, -747989353.5, -1495978707.0]) +远日点:34.890 AU, 近日点:0.578 AU, 天数: 13596 +Index:113 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.816, 4.705, 8.85], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.015 AU, 近日点:0.577 AU, 天数: 13669 +Index:114 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.816, 4.705, 8.86], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.142 AU, 近日点:0.575 AU, 天数: 13742 +Index:115 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.816, 4.705, 8.87], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.269 AU, 近日点:0.575 AU, 天数: 13815 +Index:116 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.816, 4.715, 8.83], init_position=[0, -747989353.5, -1495978707.0]) +远日点:34.831 AU, 近日点:0.579 AU, 天数: 13566 +Index:117 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.816, 4.715, 8.84], init_position=[0, -747989353.5, -1495978707.0]) +远日点:34.956 AU, 近日点:0.578 AU, 天数: 13639 +Index:118 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.816, 4.715, 8.85], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.083 AU, 近日点:0.576 AU, 天数: 13712 +Index:119 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.816, 4.715, 8.86], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.208 AU, 近日点:0.575 AU, 天数: 13785 +Index:120 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.816, 4.715, 8.87], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.337 AU, 近日点:0.575 AU, 天数: 13858 +Index:121 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.816, 4.725, 8.83], init_position=[0, -747989353.5, -1495978707.0]) +远日点:34.898 AU, 近日点:0.579 AU, 天数: 13602 +Index:122 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.816, 4.725, 8.84], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.023 AU, 近日点:0.578 AU, 天数: 13675 +Index:123 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.816, 4.725, 8.85], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.149 AU, 近日点:0.576 AU, 天数: 13748 +Index:124 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.816, 4.725, 8.86], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.276 AU, 近日点:0.576 AU, 天数: 13821 +Index:125 --------------------------- +北京时间： 1982-09-24 00:00:00+08:00 +HalleyCometParams(start_time='1982-09-24 00:00:00', init_velocity=[-2.816, 4.725, 8.87], init_position=[0, -747989353.5, -1495978707.0]) +远日点:35.405 AU, 近日点:0.575 AU, 天数: 13894 + +Process finished with exit code 0 diff --git a/sim_lab/halley_comet_minimize.py b/sim_lab/halley_comet_minimize.py new file mode 100644 index 0000000000000000000000000000000000000000..1b7c0ec8e9412e195c2c23821f02fdad87a704bc --- /dev/null +++ b/sim_lab/halley_comet_minimize.py @@ -0,0 +1,94 @@ +import numpy as np + +from common.consts import AU +from sim_scenes.solar_system.halley_comet_lib import HalleyCometParams +from sim_lab.halley_comet_research_calc import calc_simulator + + +def f(x, y, z): + params = HalleyCometParams( + start_time='1982-09-24 00:00:00', + # init_velocity=[-2.836 + (x / 1000), 4.705 + (y / 1000), 8.85 + (z / 1000)], + init_velocity=[x, y, z], + init_position=[0, -5 * AU, -10 * AU] + ) + return calc_simulator(params) + + +from scipy.optimize import minimize +import numpy as np + + + +# 定义损失函数 +def loss_function(x, target): + # 定义损失函数，这里使用简单的平方差，并调整权重 + weights = np.array([10, 1, 0.1]) # 调整权重，使近日点和远日点的偏差对最终结果的影响更大 + return np.sum(weights * (x - target) ** 2) + + +# 定义目标函数 +def target_function(v): + vx, vy, vz = v + result = f(vx, vy, vz) + # 计算损失函数，这里假设目标是 [0.586, 35.1, 13817] + target = np.array([0.586, 35.1, 13817]) + loss = loss_function(result, target) + return loss + + +# 设置初始速度 +initial_velocity = [-2.816, 4.725, 8.87] # 请替换为实际的初始值 +# 设置初始速度时添加随机扰动 +initial_velocity = [-2.816 + np.random.uniform(-1, 1), 4.725 + np.random.uniform(-1, 1), 8.87 + np.random.uniform(-1, 1)] + + +# 运行优化算法 +result = minimize(target_function, initial_velocity, method='trust-constr', options={'maxiter': 10000, 'disp': True}, + bounds=[(-100, 100), (-100, 100), (-100, 100)]) + +# 打印最终结果 +final_velocity = result.x +print(f"Optimal Velocity: {final_velocity}") + +# +# # # 定义目标函数 +# # def f(vx, vy, vz): +# # # 这里需要根据你的实际情况定义哈雷彗星轨道的计算 +# # # 返回近日点距离（单位：AU）、远日点距离（单位：AU）、以及哈雷彗星近日点到远日点花了多长时间（单位：天） +# # # 这里使用一个简单的示例，你需要根据实际情况修改 +# # return np.array([0.586, 35.1, 13817]) +# +# # 定义损失函数 +# def loss_function(x, target): +# # 定义损失函数，这里使用简单的平方差 +# return (x - target) ** 2 +# +# +# # 定义目标函数 +# def target_function(v): +# vx, vy, vz = v +# result = f(vx, vy, vz) +# # 计算损失函数，这里假设目标是 [0.586, 35.1, 13817] +# target = np.array([0.586, 35.1, 13817]) +# loss = np.sum([loss_function(result[i], target[i]) for i in range(len(result))]) +# return loss +# +# +# # 设置初始速度 +# initial_velocity = [-2.816, 4.725, 8.87] # 请替换为实际的初始值 +# +# # 运行优化算法 +# # result = minimize(target_function, initial_velocity, method='BFGS', +# # options={'maxiter': 1000, 'disp': True}) # 使用BFGS优化算法，你可以尝试其他算法 +# +# # 运行优化算法 +# # result = minimize(target_function, initial_velocity, method='L-BFGS-B', options={'maxiter': 10000, 'disp': True}) +# +# # 运行优化算法 +# result = minimize(target_function, initial_velocity, method='trust-constr', options={'maxiter': 10000, 'disp': True}, +# bounds=[(-100, 100), (-100, 100), (-100, 100)]) +# +# # 打印最终结果 +# final_velocity = result.x +# print(f"Optimal Velocity: {final_velocity}") diff --git a/sim_lab/halley_comet_plt_draw.py b/sim_lab/halley_comet_plt_draw.py new file mode 100644 index 0000000000000000000000000000000000000000..f969566e8829abcf889f31e50cf3f4e064f358c0 --- /dev/null +++ b/sim_lab/halley_comet_plt_draw.py @@ -0,0 +1,142 @@ +import numpy as np +import matplotlib.pyplot as plt +from common.consts import AU +from sim_scenes.solar_system.halley_comet_lib import HalleyCometParams +from sim_lab.halley_comet_research_calc import calc_simulator, target_function + +# 支持中文 +plt.rcParams['font.sans-serif'] = ['SimHei'] # 用来正常显示中文标签 +plt.rcParams['axes.unicode_minus'] = False # 用来正常显示负号 + + +def plt_draw_3d(): + # idx = 0 + r = 1 + xs = [] + ys = [] + zs = [] + a_s = [] + b_s = [] + c_s = [] + for xi in range(-r, r + 1): + for yi in range(-r, r + 1): + for zi in range(-r, r + 1): + x, y, z = -2.47 + (xi / 100), 5.1 + (yi / 100), 8.73 + (zi / 100) + xs.append(x) + ys.append(y) + zs.append(z) + a, b, c = target_function(x, y, z) + # a,b,c = x,y,z + a_s.append(a) + b_s.append(b) + c_s.append(c) + # 把 x y z 和 结果 a b c 进行可视化，完成代码，方便查看 x y z 和 结果 a b c 之间的关系 + fig = plt.figure() + ax = fig.add_subplot(111, projection='3d') + ax.plot(xs, ys, zs, label="IN") + ax.plot(a_s, b_s, c_s, label="OUT") + # ax.scatter(xs, ys, zs, c='r', marker='o') # 绘制三维散点图，使用红色标记，可以根据需要更改颜色和标记类型 + # ax.scatter(a_s, b_s, c_s, c='b', marker='o') # 在同样的图上，使用蓝色标记绘制函数的结果点 + ax.legend() + plt.show() + + +import matplotlib.pyplot as plt + + +def plot_twin(t, _y1, _y2, _y3, xlabel, _ylabel1, _ylabel2, _ylabel3): + fig, ax1 = plt.subplots() + color = 'tab:blue' + ax1.set_xlabel(xlabel) + ax1.set_ylabel(_ylabel1, color=color) + ax1.plot(t, _y1, color=color) + ax1.tick_params(axis='y', labelcolor=color) + + ax2 = ax1.twinx() # 创建共用x轴的第二个y轴 + + color = 'tab:red' + ax2.set_ylabel(_ylabel2, color=color) + ax2.plot(t, _y2, color=color) + ax2.tick_params(axis='y', labelcolor=color) + + ax3 = ax1.twinx() # 创建共用x轴的第二个y轴 + + color = 'tab:green' + ax3.set_ylabel(_ylabel3, color=color) + ax3.plot(t, _y3, color=color) + ax3.tick_params(axis='y', labelcolor=color) + + # fig.tight_layout() + plt.show() + + +def plt_draw_2d(): + import numpy as np + t = np.array([1, 2, 3]) + # 创建模拟数据 + data1 = np.array([35.032, 35.033, 35.034]) + data2 = np.array([0.477, 0.477, 0.477]) + data3 = np.array([13615, 13615, 13615]) + plot_twin(t, data1, data2, data3, 'exp', 'sine', 'year') + + +# target_function(-2.835, 4.71, 8.86) # 远日点:35.243 AU, 近日点:0.584 AU, 天数: 13809 +# target_function(-2.835, 4.73, 8.85) # OK 远日点:35.251 AU, 近日点:0.585 AU, 天数: 13815 + +R = 3 +X, Y, Z = -2.835, 4.73, 8.85 +X, Y, Z = -2.835, 4.72, 8.847 +# [-2.816, 4.725, 8.84] + +def plt_draw_x(): + x_s = [] + a_s = [] + b_s = [] + c_s = [] + for x in range(-R, R + 1): + x, y, z = X + (x / 100), Y, Z + x_s.append(x) + a, b, c = target_function(x, y, z) + a_s.append(a) + b_s.append(b) + c_s.append(c) + + plot_twin(x_s, a_s, b_s, c_s, 'X', '远日点', '近日点', '天数') + + +def plt_draw_y(): + y_s = [] + a_s = [] + b_s = [] + c_s = [] + for y in range(-R, R + 1): + x, y, z = X, Y + (y / 100), Z + y_s.append(y) + a, b, c = target_function(x, y, z) + a_s.append(a) + b_s.append(b) + c_s.append(c) + + plot_twin(y_s, a_s, b_s, c_s, 'Y', '远日点', '近日点', '天数') + + +def plt_draw_z(): + z_s = [] + a_s = [] + b_s = [] + c_s = [] + for z in range(-R, R + 1): + x, y, z = X, Y, Z + (z / 100) + z_s.append(z) + a, b, c = target_function(x, y, z) + a_s.append(a) + b_s.append(b) + c_s.append(c) + + plot_twin(z_s, a_s, b_s, c_s, 'Z', '远日点', '近日点', '天数') + + +if __name__ == '__main__': + plt_draw_x() + plt_draw_y() + plt_draw_z() diff --git a/sim_lab/halley_comet_research_calc.py b/sim_lab/halley_comet_research_calc.py index 6a898297cae62839c65937db246d0fd01e53b188..27335aeb31d810ce926646918343b30f6e02a790 100644 --- a/sim_lab/halley_comet_research_calc.py +++ b/sim_lab/halley_comet_research_calc.py @@ -5,7 +5,18 @@ # date :2023-11-16 # link :https://gitcode.net/pythoncr/ # python_version :3.9 +# ChatGPT # ============================================================================== +""" +我把我的详细需求写给你，你帮我分析用什么方法。我要计算哈雷彗星在太阳系中某个位置上面的初始速度（vx, vy, vz） ， +通过万有引力的计算，哈雷彗星会在一个轨道上有2个点，近日点和远日点。哈雷彗星从近日点和远日点花了多长时间。 +我现在定义了一个函数，f(vx,vy,vz) 参数就是初始速度（vx, vy, vz）， +这个函数f返回3个值，分别是 近日点距离（单位：AU）、远日点距离（单位：AU）、以及哈雷彗星近日点到远日点花了多长时间（单位：天）。 +我现在就是希望如果函数f返回值是 0.586AU, 35.1AU, 13817天，那么 vx,vy,vz的值是多少？ 以上的需求，我用什么算法可以实现。 + +""" +import threading + from bodies import Body, Sun, Earth, Moon from objs import Obj, Satellite, Satellite2 from common.consts import SECONDS_PER_HOUR, SECONDS_PER_DAY, SECONDS_PER_WEEK, SECONDS_PER_MONTH @@ -110,6 +121,8 @@ def calc_simulator(params): sim.total_times = 0 sim.build() + sim.result = None + def on_ready(context: CalcContext): # for body in sim.bodies: # body.reset() @@ -198,8 +211,8 @@ def calc_simulator(params): 2023年12月9日 - 1986年2月9日 = 13817 35.1AU 0.586AU """ - if abs(comet_peri - 0.586) < 0.01 and abs(comet_aphel - 35.1) < 0.1 and (diff_days - 13817) < 2: - print("OK ......",comet_peri,comet_aphel,diff_days) + if abs(comet_peri - 0.586) < 0.01 and abs(comet_aphel - 35.1) < 0.1 and abs(diff_days - 13817) < 2: + print("OK ......", comet_peri, comet_aphel, diff_days) # print("year days:", time_data.years, time_data.days) # print("远日点: %.3f AU(%s)" % (sim.comet_aphel / AU, sim.comet_aphel_dt)) # print("近日点: %.3f AU(%s)" % (sim.comet_peri / AU, sim.comet_peri_dt)) @@ -207,9 +220,14 @@ def calc_simulator(params): # print("远\近日天数: %i 天 " % (sim.comet_aphel_dt[2] - sim.comet_peri_dt[2])) def on_finished(context: CalcContext): + global result # print("on_finished", context) time_data = context.get_param("time_data") log(time_data) + comet_aphel = sim.comet_aphel / AU + comet_peri = sim.comet_peri / AU + diff_days = sim.comet_aphel_dt[2] - sim.comet_peri_dt[2] + sim.result = (comet_aphel, comet_peri, float(diff_days)) CalcSimulator.on_ready_subscription(on_ready) CalcSimulator.on_after_evolve_subscription(after_evolve) @@ -218,13 +236,58 @@ def calc_simulator(params): calc_run(sim.bodies, SECONDS_PER_YEAR, evolve_next=evolve_next) + return sim.result -if __name__ == '__main__': + +def target_function(x, y, z): params = HalleyCometParams( start_time='1982-09-24 00:00:00', - init_velocity=[-2.836, 4.705, 8.85], + # init_velocity=[-2.836 + (x / 1000), 4.705 + (y / 1000), 8.85 + (z / 1000)], + init_velocity=[x, y, z], init_position=[0, -5 * AU, -10 * AU] ) + return calc_simulator(params) + + +if __name__ == '__main__': + pass + # 2023年12月9日 - 1986年2月9日 = 13817 + # 35.1AU 0.586AU + # target_function(-2.836, 4.705, 8.85) + # target_function(-2.816, 4.725, 8.86) # 远日点:35.276 AU, 近日点:0.576 AU, 天数: 13821 + # target_function(-2.816, 4.725, 8.855) # 远日点:35.212 AU, 近日点:0.576 AU, 天数: 13785 + # target_function(-2.816, 4.725, 8.857) # 远日点:35.237 AU, 近日点:0.576 AU, 天数: 13797 + # target_function(-2.816, 4.725, 8.858) # 远日点:35.251 AU, 近日点:0.576 AU, 天数: 13809 + # target_function(-2.81, 4.725, 8.858) # 远日点:35.230 AU, 近日点:0.573 AU, 天数: 13791 + # target_function(-2.81, 4.735, 8.858) # 远日点:35.297 AU, 近日点:0.573 AU, 天数: 13834 + # target_function(-2.81, 4.73, 8.858) # 远日点:35.264 AU, 近日点:0.573 AU, 天数: 13815 + # target_function(-2.820, 4.73, 8.858) # 远日点:35.298 AU, 近日点:0.578 AU, 天数: 13834 + # target_function(-2.820, 4.73, 8.85) # 远日点:35.196 AU, 近日点:0.578 AU, 天数: 13779 + # target_function(-2.825, 4.73, 8.85) # 远日点:35.215 AU, 近日点:0.581 AU, 天数: 13791 + # target_function(-2.825, 4.73, 8.848) # 远日点:35.190 AU, 近日点:0.581 AU, 天数: 13773 + # x, y, z => a, b, c + # x=a,b,c + # y=a,c + # z=a,c + # target_function(-2.835, 4.72, 8.847) # 远日点:35.144 AU, 近日点:0.586 AU, 天数: 13748 + # target_function(-2.835, 4.71, 8.86) # 远日点:35.243 AU, 近日点:0.584 AU, 天数: 13809 + # target_function(-2.835, 4.73, 8.85) # OK 远日点:35.251 AU, 近日点:0.585 AU, 天数: 13815 + a, b, c = target_function(-2.47, 5.13, 8.73) + + + # [-2.5, 5.02, 8.763] 远日点:35.032 AU, 近日点:0.477 AU, 天数: 13615 + # [-2.835, 4.81, 8.8] 远日点:35.161 AU, 近日点:0.591 AU, 天数: 13754 + # [-2.835, 4.8, 8.81] 远日点:35.220 AU, 近日点:0.590 AU, 天数: 13797 + # [-2.835, 4.8, 8.8], 远日点:35.092 AU, 近日点:0.590 AU, 天数: 13718 + # [-2.835, 4.83, 8.8] 远日点:35.299 AU, 近日点:0.592 AU, 天数: 13846 + # [-2.835, 4.75, 8.83] 远日点:35.132 AU, 近日点:0.588 AU, 天数: 13742 + # [-2.835, 4.74, 8.84] 远日点:35.191 AU, 近日点:0.587 AU, 天数: 13779 + # [-2.835, 4.79, 8.8] 远日点:35.025 AU, 近日点:0.589 AU, 天数: 13682 + # params = HalleyCometParams( + # start_time='1982-09-24 00:00:00', + # init_velocity=[-2.836, 4.705, 8.85], + # init_position=[0, -5 * AU, -10 * AU] + # ) # year days: 42 6 # 远日点: 35.086 AU((40, 346, 14946.8522)) # 近日点: 0.586 AU((3, 133, 1228.3418)) @@ -233,16 +296,16 @@ if __name__ == '__main__': # for x in range(-2,2): # for y in range(-2, 2): # for z in range(-2, 2): - idx = 0 - r = 2 - for x in range(-r, r+1): - for y in range(-r, r+1): - for z in range(-r, r+1): - params = HalleyCometParams( - start_time='1982-09-24 00:00:00', - init_velocity=[-2.836 + (x / 100), 4.705 + (y / 100), 8.85 + (z / 100)], - init_position=[0, -5 * AU, -10 * AU] - ) - idx += 1 - print(f"Index:{idx} ---------------------------") - calc_simulator(params) + + # def t(): + # global idx + # params = HalleyCometParams( + # start_time='1982-09-24 00:00:00', + # init_velocity=[-2.836 + (x / 100), 4.705 + (y / 100), 8.85 + (z / 100)], + # init_position=[0, -5 * AU, -10 * AU] + # ) + # idx += 1 + # print(f"Index:{idx} ---------------------------") + # calc_simulator(params) + # # threading.Thread(target=t).start() + # t() diff --git a/sim_scenes/solar_system/halley_comet_sim.py b/sim_scenes/solar_system/halley_comet_sim.py index bd25ba888d83428df0c25ae0fb0e3517327b497e..bd16b807a00a70fef51477c75931922f831ceb2f 100644 --- a/sim_scenes/solar_system/halley_comet_sim.py +++ b/sim_scenes/solar_system/halley_comet_sim.py @@ -227,7 +227,8 @@ if __name__ == '__main__': # 2019年5月6日 34.772 params = HalleyCometParams( start_time='1982-09-24 00:00:00', - init_velocity=[-2.836, 4.705, 8.85], + init_velocity=[-2.835, 4.72, 8.847], + # init_velocity=[-2.836, 4.705, 8.85], init_position=[0, -5 * AU, -10 * AU] )