我正在尝试尽可能快地执行以下脚本:
- 从文本文件中读取浮点列表
- 求解微分方程
- 将一些结果输出到
stdout
。
应该认为它是测试优化工具的“黑匣子”,该优化工具与脚本作为系统命令交互。因此,该工具调用python blackbox.py x.txt
(其中x.txt
是评估黑盒的点,此处有0到3.5之间的500个元素),并且输出是目标函数以及对违反约束的度量。 (例如-1.007 -0.067 -4.471
)。
在这种情况下,我正在尝试解决Goddard rocket problem。
这是我所拥有的代码,以及我所知道的加快技巧的所有技巧:
import sys
from numpy import max,sqrt,exp,errstate # Avoid resolving NumPy methods (as np.max) within the code
from scipy.integrate import solve_ivp
TF = 0.2 # Final time
# Read point
with open(sys.argv[1],'r') as f:
x_u = [float(xi) for xi in f.read().split()]
i_last_x_u = len(x_u)-1 # Index of last element
# Function to solve
def f(t,y):
h,v,m = y # Decompose by rows
D_hv = 0.5*620*v**2*exp(500*(1-h)) # Drag
u_t = x_u[int(t/TF*i_last_x_u)] # "Interpolation" on input vector
return [
v,# h'
(u_t - D_hv)/m - 1/h**2,# v'
-u_t/0.5,# m'
]
# Solver call
with errstate(all='ignore'): # Errors are catched by looking at sol.success
sol = solve_ivp(f,t_span=(0,TF),y0=[1,1],vectorized=True)
if sol.success:
# Solution
h_sol,v_sol,m_sol = sol.y # Decompose by rows
# Objective function: maximum height
# Minus sign since NOMAD tries to minimize obj
obj_h_max = -max(h_sol)
# Constraint: Final mass must be >0.6
pb_m = 0.6 - m_sol[-1]
# Constraint: Dynamic pressure
pb_q = max( v_sol - sqrt(2*10/exp(500*(1-h_sol))) )
# Output solution
print(f'{obj_h_max} {pb_m} {pb_q}')
else:
# Solver failed
# Non-zero exit status is flagged as a failed evaluation
print('0 0 0')
sys.exit(1)
下面是Spyder探查器的输出。看来罪魁祸首是_find_and_load
。
如果不使用C编写整个代码,我该怎么做才能加快速度?
谢谢!