Scipy.optimize.minimize 关于轨迹优化（成本函数最小化）
Scipy.optimize.minimize on trajectory optimization (Cost Functional Minimization)
我在这里发布了一个问题，询问一般如何解决这个问题：
Trajectory Optimization for "Rocket" using scipy.optimize.minimize
理想情况下，我只想最小化最后的时间，但我无法让优化器将时间附加到可以适当调整的变量上，所以我决定只尝试最小化 u^2现在。
代码如下：
# Code
t_f = 1.0
t = np.linspace(0., t_f, num = 10) # Time array for 1 second into the future with 0.01 increment
u = np.zeros(t.size) + 650
print(u)
g = -650
initial_position = 0
initial_velocity = 0
final_position = 100
final_velocity = 100
def car_dynamics(x):
# Create time vector
# t = np.linspace(0., t_f, num = 100) # Time array for 1 second into the future with 0.01 increment
# Integrate over entire time to find v as a function of t
a = x + g
v = int.cumtrapz(a, t, initial = 0) + initial_velocity
# Integrate v(t) to get s(t)
s = int.cumtrapz(v, t, initial = 0) + initial_position
return s, v
def constraint1(x): # Final state constraints (Boundary conditions)
s, v = car_dynamics(x)
print('c1', s[0] - initial_position)
return s[0] - initial_position
def constraint2(x): # Initial state constraints (initial conditions)
s, v = car_dynamics(x)
print('c2', v[0] - initial_velocity)
return v[0] - initial_velocity
def constraint3(x):
s, v = car_dynamics(x)
print('c3', s[-1] - final_position)
return s[-1] - final_position
def constraint4(x):
s, v = car_dynamics(x)
print('c4', v[-1] - final_velocity)
return v[-1] - final_velocity
def constraint5(x):
return x - 1000
def objective(x):
u2 = np.square(x)
return np.sum(u2)
cons = [{'type':'eq', 'fun':constraint1},
{'type':'eq', 'fun':constraint2},
{'type':'eq', 'fun':constraint3},
{'type':'eq', 'fun':constraint4}]
# {'type':'ineq', 'fun':constraint5}]
result = minimize(objective, u, constraints = cons, method = 'SLSQP', options={'eps':500, 'maxiter':1000, 'ftol':0.001, 'disp':True})
print(result)
代码运行但优化器失败。这是输出的错误。
c1 0.0
c2 0.0
c3 -100.0
c4 -100.0
c1 0.0
c2 0.0
c3 -100.0
c4 -100.0
c1 0.0
c1 0.0
c1 0.0
c1 0.0
c1 0.0
c1 0.0
c1 0.0
c1 0.0
c1 0.0
c1 0.0
c1 0.0
c2 0.0
c2 0.0
c2 0.0
c2 0.0
c2 0.0
c2 0.0
c2 0.0
c2 0.0
c2 0.0
c2 0.0
c2 0.0
c3 -100.0
c3 -73.76543209876543
c3 -50.617283950617285
c3 -56.79012345679013
c3 -62.96296296296296
c3 -69.1358024691358
c3 -75.30864197530863
c3 -81.4814814814815
c3 -87.65432098765432
c3 -93.82716049382715
c3 -98.45679012345678
c4 -100.0
c4 -72.22222222222223
c4 -44.44444444444445
c4 -44.44444444444445
c4 -44.44444444444445
c4 -44.444444444444436
c4 -44.44444444444445
c4 -44.44444444444448
c4 -44.44444444444445
c4 -44.44444444444442
c4 -72.22222222222221
Singular matrix C in LSQ subproblem (Exit mode 6)
Current function value: 4225000.0
Iterations: 1
Function evaluations: 12
Gradient evaluations: 1
fun: 4225000.0
jac: array([1800., 1800., 1800., 1800., 1800., 1800., 1800., 1800., 1800.,
1800.])
message: 'Singular matrix C in LSQ subproblem'
nfev: 12
nit: 1
njev: 1
status: 6
success: False
x: array([650., 650., 650., 650., 650., 650., 650., 650., 650., 650.])
似乎在一定数量的迭代中没有满足约束条件。我应该切换 objective 函数以包含最终速度和最终位置吗？我尝试了不同的步长大小，并且尝试了不同的退出代码。
有没有更好的方法来使用这个函数来实现我想要得到的东西？我试图在从 t0 到 t_f 的整个时间间隔内获取控制向量 u(t)，这样我就可以将这些命令发送到火箭以获得最佳控制。现在我已经将优化简化为单轴，只是为了学习如何使用该功能。但是如你所见，我没有成功。
类似的例子会非常有帮助，我对其他优化方法持开放态度，只要它们是数字的，并且相对较快，因为我计划最终将其实现为实时模型预测控制器。
您的模型同时具有代数方程和微分方程。您需要 DAE 求解器来求解上述隐式 ODE 函数。我知道的一个这样的包是gekko。 (https://github.com/BYU-PRISM/GEKKO)
Gekko 专注于线性、混合整数和非线性优化问题的动态优化。
下面是一个最小化最终时间的火箭发射问题示例。可在 http://apmonitor.com/wiki/index.php/Apps/RocketLaunch
import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO
# create GEKKO model
m = GEKKO()
# scale 0-1 time with tf
m.time = np.linspace(0,1,101)
# options
m.options.NODES = 6
m.options.SOLVER = 3
m.options.IMODE = 6
m.options.MAX_ITER = 500
m.options.MV_TYPE = 0
m.options.DIAGLEVEL = 0
# final time
tf = m.FV(value=1.0,lb=0.1,ub=100)
tf.STATUS = 1
# force
u = m.MV(value=0,lb=-1.1,ub=1.1)
u.STATUS = 1
u.DCOST = 1e-5
# variables
s = m.Var(value=0)
v = m.Var(value=0,lb=0,ub=1.7)
mass = m.Var(value=1,lb=0.2)
# differential equations scaled by tf
m.Equation(s.dt()==tf*v)
m.Equation(mass*v.dt()==tf*(u-0.2*v**2))
m.Equation(mass.dt()==tf*(-0.01*u**2))
# specify endpoint conditions
m.fix(s, pos=len(m.time)-1,val=10.0)
m.fix(v, pos=len(m.time)-1,val=0.0)
# minimize final time
m.Obj(tf)
# Optimize launch
m.solve()
print('Optimal Solution (final time): ' + str(tf.value[0]))
# scaled time
ts = m.time * tf.value[0]
# plot results
plt.figure(1)
plt.subplot(4,1,1)
plt.plot(ts,s.value,'r-',linewidth=2)
plt.ylabel('Position')
plt.legend(['s (Position)'])
plt.subplot(4,1,2)
plt.plot(ts,v.value,'b-',linewidth=2)
plt.ylabel('Velocity')
plt.legend(['v (Velocity)'])
plt.subplot(4,1,3)
plt.plot(ts,mass.value,'k-',linewidth=2)
plt.ylabel('Mass')
plt.legend(['m (Mass)'])
plt.subplot(4,1,4)
plt.plot(ts,u.value,'g-',linewidth=2)
plt.ylabel('Force')
plt.legend(['u (Force)'])
plt.xlabel('Time')
plt.show()
我在这里发布了一个问题，询问一般如何解决这个问题：
Trajectory Optimization for "Rocket" using scipy.optimize.minimize
理想情况下，我只想最小化最后的时间，但我无法让优化器将时间附加到可以适当调整的变量上，所以我决定只尝试最小化 u^2现在。
代码如下：
# Code
t_f = 1.0
t = np.linspace(0., t_f, num = 10) # Time array for 1 second into the future with 0.01 increment
u = np.zeros(t.size) + 650
print(u)
g = -650
initial_position = 0
initial_velocity = 0
final_position = 100
final_velocity = 100
def car_dynamics(x):
# Create time vector
# t = np.linspace(0., t_f, num = 100) # Time array for 1 second into the future with 0.01 increment
# Integrate over entire time to find v as a function of t
a = x + g
v = int.cumtrapz(a, t, initial = 0) + initial_velocity
# Integrate v(t) to get s(t)
s = int.cumtrapz(v, t, initial = 0) + initial_position
return s, v
def constraint1(x): # Final state constraints (Boundary conditions)
s, v = car_dynamics(x)
print('c1', s[0] - initial_position)
return s[0] - initial_position
def constraint2(x): # Initial state constraints (initial conditions)
s, v = car_dynamics(x)
print('c2', v[0] - initial_velocity)
return v[0] - initial_velocity
def constraint3(x):
s, v = car_dynamics(x)
print('c3', s[-1] - final_position)
return s[-1] - final_position
def constraint4(x):
s, v = car_dynamics(x)
print('c4', v[-1] - final_velocity)
return v[-1] - final_velocity
def constraint5(x):
return x - 1000
def objective(x):
u2 = np.square(x)
return np.sum(u2)
cons = [{'type':'eq', 'fun':constraint1},
{'type':'eq', 'fun':constraint2},
{'type':'eq', 'fun':constraint3},
{'type':'eq', 'fun':constraint4}]
# {'type':'ineq', 'fun':constraint5}]
result = minimize(objective, u, constraints = cons, method = 'SLSQP', options={'eps':500, 'maxiter':1000, 'ftol':0.001, 'disp':True})
print(result)
代码运行但优化器失败。这是输出的错误。
c1 0.0
c2 0.0
c3 -100.0
c4 -100.0
c1 0.0
c2 0.0
c3 -100.0
c4 -100.0
c1 0.0
c1 0.0
c1 0.0
c1 0.0
c1 0.0
c1 0.0
c1 0.0
c1 0.0
c1 0.0
c1 0.0
c1 0.0
c2 0.0
c2 0.0
c2 0.0
c2 0.0
c2 0.0
c2 0.0
c2 0.0
c2 0.0
c2 0.0
c2 0.0
c2 0.0
c3 -100.0
c3 -73.76543209876543
c3 -50.617283950617285
c3 -56.79012345679013
c3 -62.96296296296296
c3 -69.1358024691358
c3 -75.30864197530863
c3 -81.4814814814815
c3 -87.65432098765432
c3 -93.82716049382715
c3 -98.45679012345678
c4 -100.0
c4 -72.22222222222223
c4 -44.44444444444445
c4 -44.44444444444445
c4 -44.44444444444445
c4 -44.444444444444436
c4 -44.44444444444445
c4 -44.44444444444448
c4 -44.44444444444445
c4 -44.44444444444442
c4 -72.22222222222221
Singular matrix C in LSQ subproblem (Exit mode 6)
Current function value: 4225000.0
Iterations: 1
Function evaluations: 12
Gradient evaluations: 1
fun: 4225000.0
jac: array([1800., 1800., 1800., 1800., 1800., 1800., 1800., 1800., 1800.,
1800.])
message: 'Singular matrix C in LSQ subproblem'
nfev: 12
nit: 1
njev: 1
status: 6
success: False
x: array([650., 650., 650., 650., 650., 650., 650., 650., 650., 650.])
似乎在一定数量的迭代中没有满足约束条件。我应该切换 objective 函数以包含最终速度和最终位置吗？我尝试了不同的步长大小，并且尝试了不同的退出代码。
有没有更好的方法来使用这个函数来实现我想要得到的东西？我试图在从 t0 到 t_f 的整个时间间隔内获取控制向量 u(t)，这样我就可以将这些命令发送到火箭以获得最佳控制。现在我已经将优化简化为单轴，只是为了学习如何使用该功能。但是如你所见，我没有成功。
类似的例子会非常有帮助，我对其他优化方法持开放态度，只要它们是数字的，并且相对较快，因为我计划最终将其实现为实时模型预测控制器。
您的模型同时具有代数方程和微分方程。您需要 DAE 求解器来求解上述隐式 ODE 函数。我知道的一个这样的包是gekko。 (https://github.com/BYU-PRISM/GEKKO) Gekko 专注于线性、混合整数和非线性优化问题的动态优化。
下面是一个最小化最终时间的火箭发射问题示例。可在 http://apmonitor.com/wiki/index.php/Apps/RocketLaunch
import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO
# create GEKKO model
m = GEKKO()
# scale 0-1 time with tf
m.time = np.linspace(0,1,101)
# options
m.options.NODES = 6
m.options.SOLVER = 3
m.options.IMODE = 6
m.options.MAX_ITER = 500
m.options.MV_TYPE = 0
m.options.DIAGLEVEL = 0
# final time
tf = m.FV(value=1.0,lb=0.1,ub=100)
tf.STATUS = 1
# force
u = m.MV(value=0,lb=-1.1,ub=1.1)
u.STATUS = 1
u.DCOST = 1e-5
# variables
s = m.Var(value=0)
v = m.Var(value=0,lb=0,ub=1.7)
mass = m.Var(value=1,lb=0.2)
# differential equations scaled by tf
m.Equation(s.dt()==tf*v)
m.Equation(mass*v.dt()==tf*(u-0.2*v**2))
m.Equation(mass.dt()==tf*(-0.01*u**2))
# specify endpoint conditions
m.fix(s, pos=len(m.time)-1,val=10.0)
m.fix(v, pos=len(m.time)-1,val=0.0)
# minimize final time
m.Obj(tf)
# Optimize launch
m.solve()
print('Optimal Solution (final time): ' + str(tf.value[0]))
# scaled time
ts = m.time * tf.value[0]
# plot results
plt.figure(1)
plt.subplot(4,1,1)
plt.plot(ts,s.value,'r-',linewidth=2)
plt.ylabel('Position')
plt.legend(['s (Position)'])
plt.subplot(4,1,2)
plt.plot(ts,v.value,'b-',linewidth=2)
plt.ylabel('Velocity')
plt.legend(['v (Velocity)'])
plt.subplot(4,1,3)
plt.plot(ts,mass.value,'k-',linewidth=2)
plt.ylabel('Mass')
plt.legend(['m (Mass)'])
plt.subplot(4,1,4)
plt.plot(ts,u.value,'g-',linewidth=2)
plt.ylabel('Force')
plt.legend(['u (Force)'])
plt.xlabel('Time')
plt.show()