Python超人-宇宙模拟器

bd901895 · 三月三net · 2342b5c6 · bd901895 · bd901895 · bd901895
6 changed file
--- a/sim_lab/halley_comet_gradient_descent.py
+++ b/sim_lab/halley_comet_gradient_descent.py
+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}")
--- a/sim_lab/halley_comet_log.txt
+++ b/sim_lab/halley_comet_log.txt
--- a/sim_lab/halley_comet_minimize.py
+++ b/sim_lab/halley_comet_minimize.py
+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}")
--- a/sim_lab/halley_comet_plt_draw.py
+++ b/sim_lab/halley_comet_plt_draw.py
+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()
--- 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()
--- 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]
    )