ロケット垂直上昇問題 SciPy

Scipyとは

SciPyは、Pythonで科学技術計算を行うためのオープンソースのライブラリです。

NumPyを基盤としており、高度な数学関数、線形代数、最適化、統計、信号処理、画像処理などの機能を提供しています。

SciPyは、科学技術計算において広く使用されており、データ解析や機械学習の分野でも重要な役割を果たしています。

SciPyの主な機能は以下の通りです。

1. 線形代数（scipy.linalg）:

行列の分解や逆行列の計算、線形方程式の解法など、線形代数の基本的な操作を提供します。

2. 最適化（scipy.optimize）:

非線形方程式の解法、最小化問題、最適化問題、制約付き最適化など、最適化アルゴリズムを提供します。

3. 統計（scipy.stats）:

確率分布、統計的検定、相関、回帰分析など、統計分析に関する機能を提供します。

4. 信号処理（scipy.signal）:

フィルタリング、ウィンドウ関数、スペクトル解析など、信号処理に関する機能を提供します。

5. 画像処理（scipy.ndimage）:

画像のフィルタリング、形態学的操作、変換など、画像処理に関する機能を提供します。

6. 積分（scipy.integrate）:

数値積分、常微分方程式の解法など、積分に関する機能を提供します。

7. 補間（scipy.interpolate）:

線形補間、スプライン補間、多次元補間など、データ間の補間に関する機能を提供します。

8. 特殊関数（scipy.special）:

ガンマ関数、ベータ関数、誤差関数など、特殊関数の計算を提供します。

これらの機能を利用することで、Pythonで科学技術計算を効率的に行うことができます。

SciPyは、研究や開発において幅広く活用されており、Pythonのエコシステムにおいて重要な位置を占めています。

ロケット垂直上昇問題

最適制御問題の一例として、簡単なロケットの垂直上昇問題を考えます。

ロケットの質量は$ m $、重力加速度は$ g $、推力は$ u $とします。

目標は、ロケットの高さを最大化することです。

制約条件として、燃料の消費量に制限があります。

この問題をSciPyを使って解きます。

まず、必要なライブラリをインポートします。

1
2
3
import numpy as np
from scipy.integrate import solve_ivp
from scipy.optimize import minimize

次に、ロケットの運動方程式を定義します。

1
2
3
4
5
6
7
def rocket_dynamics(t, y, u):
    m = 1.0  # mass of the rocket
    g = 9.81  # gravitational acceleration

    # y[0] is the height, y[1] is the velocity
    dydt = [y[1], u / m - g]
    return dydt

制御変数uを決定するために、最適化問題を定義します。

1
2
3
4
5
6
7
8
9
def objective(u):
    # Initial conditions
    y0 = [0, 0]  # initial height and velocity

    # Integrate the rocket dynamics with the given control input u
    sol = solve_ivp(rocket_dynamics, [0, 1], y0, args=(u,), t_eval=[1])

    # The objective is to maximize the height, so we minimize the negative height
    return -sol.y[0][-1]

燃料の消費量に関する制約条件を定義します。

1
2
3
def fuel_constraint(u):
    max_fuel = 10  # maximum fuel consumption
    return max_fuel - abs(u)

制約条件を設定し、最適化問題を解きます。

1
2
constraints = {'type': 'ineq', 'fun': fuel_constraint}
result = minimize(objective, 0, constraints=constraints)

最適な制御入力と最大高度を表示します。

1
2
3
4
optimal_u = result.x[0]
max_height = -result.fun
print(f"Optimal control input: {optimal_u}")
print(f"Maximum height: {max_height}")

このコードは、ロケットの垂直上昇問題を最適制御問題として定式化し、SciPyを使って解いています。

最適な制御入力と最大高度が得られます。

[実行結果]
Optimal control input: 10.0
Maximum height: 0.0949999999999997

グラフ化

ロケットの高さと速度を時間の関数としてプロットします。

最適制御入力optimal_uを使用して、ロケットの運動方程式を解き、結果をグラフ化します。

まず、必要なライブラリをインポートします。

1
import matplotlib.pyplot as plt

次に、最適制御入力を使用してロケットの運動方程式を解きます。

1
2
3
t_eval = np.linspace(0, 1, 100)
y0 = [0, 0]  # initial height and velocity
sol = solve_ivp(rocket_dynamics, [0, 1], y0, args=(optimal_u,), t_eval=t_eval)

最後に、高さと速度を時間の関数としてプロットします。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
plt.figure(figsize=(10, 5))

plt.subplot(1, 2, 1)
plt.plot(t_eval, sol.y[0])
plt.xlabel('Time (s)')
plt.ylabel('Height (m)')
plt.title('Height vs Time')

plt.subplot(1, 2, 2)
plt.plot(t_eval, sol.y[1])
plt.xlabel('Time (s)')
plt.ylabel('Velocity (m/s)')
plt.title('Velocity vs Time')

plt.tight_layout()
plt.show()

このコードは、最適制御入力optimal_uを使用してロケットの運動方程式を解き、高さと速度を時間の関数としてプロットします。

[実行結果]

左のグラフは、ロケットの高さが時間とともにどのように変化するかを示しています。

右のグラフは、ロケットの速度が時間とともにどのように変化するかを示しています。

最適制御入力を使用すると、ロケットの高さが最大化されることがわかります。