Python GEKKO：Objective 函数未显示正确的结果
Python GEKKO: Objective function not showing correct results
我正在尝试优化推力推进系统的轨迹。控制变量是质量流率，最终的objective是最大化机器人的质量，最小化使用的推进量。轨迹类似于弹道，具有初始上升阶段和最终下降阶段。
我想我得到了一个很好的初始猜测，但是算法没有收敛。我检查了控制台中的输出，似乎 objective 函数没有正常工作，我认为这就是它没有收敛的原因。
这是我的代码
~# -*- coding: utf-8 -*-
from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
# create GEKKO model
m = GEKKO(remote=False)
# m = GEKKO()
print(m.path)
g_0 = 9.81
Isp = 46.28
m_fb=1
m_prop=0.1*m_fb
m_dry=value=m_fb-m_prop
c = 0.8
a = 0.4
g = 1.62
R_moon = 1737.4e3
J_y = 0.06666
m_dot_up=2.5*m_fb*g/(Isp*g_0)
m_dot_low=0*m_fb*g/(Isp*g_0)
theta_0=70*np.pi/180
v0= 3
x_f=np.sin(2*theta_0)*v0*v0/g
v0_x= v0*np.cos(theta_0)
y_max=v0*v0/(2*g)
m0=1
m1=(np.e**(-v0/(Isp*g_0)))
time_burn= m_fb*(m0-m1)/m_dot_up
tf=2*v0*np.sin(theta_0)/g + 2*time_burn
t0=0
nr_intervals=30
step=tf/nr_intervals
t=np.linspace(t0, tf, nr_intervals)
time_burn_node= min(range(len(t)), key=lambda i: abs(t[i]-time_burn))
m.time=t
m0=1
m1=(np.e**(-v0/(Isp*g_0)))
time_burn= m_fb*(m0-m1)/m_dot_up
time_burn_node= min(range(len(t)), key=lambda i: abs(t[i]-time_burn))
y_f=1
v_f_x = 1
v_f_y=3
gamma_f = -90*np.pi/180
alpha_f= 180*np.pi/180
alpha_dot_f=5*np.pi/180
parabola_profile= lambda x: -4*y_max/(x_f*x_f)*x*x+4*y_max*x/x_f
velocity_profile_y=lambda t: (v0*np.sin(theta_0)-g*t)
velocity_profile_x= lambda t: v0*np.cos(theta_0)
velocity_profile=lambda t: np.sqrt(velocity_profile_x(t)*velocity_profile_x(t)+velocity_profile_y(t)*velocity_profile_y(t))
x2 = m.Var(value=[np.linspace(0,v0,time_burn_node+1)[int(i)] for i in np.linspace(0,time_burn_node,time_burn_node+1)]+[velocity_profile(i) for i in t[1:int(nr_intervals)-time_burn_node]],lb=0,ub=1e3)
x3 = m.Var(value=[np.linspace(np.pi/2,theta_0,time_burn_node+1)[int(i)] for i in np.linspace(0,time_burn_node,time_burn_node+1)]+[np.arctan2(velocity_profile_y(i),velocity_profile_x(i)) for i in t[1:int(nr_intervals)-2*time_burn_node-1]]+[np.linspace(np.arctan2(velocity_profile_y(t[int(nr_intervals)-2*time_burn_node]),velocity_profile_x(t[int(nr_intervals)-2*time_burn_node])),-np.pi/2,time_burn_node+1)[int(i)] for i in np.linspace(0,time_burn_node,time_burn_node+1)],lb=-np.pi*2,ub=np.pi*2)
x0 = m.Var(value=[np.trapz(np.multiply(x2.value[0:i],np.cos(x3.value[0:i])),x=t[0:i]) for i in range(0,int(nr_intervals))],lb=0,ub=x_f+1)
x1 = m.Var(value=[np.trapz(np.multiply(x2.value[0:i],np.sin(x3.value[0:i])),x=t[0:i]) for i in range(0,int(nr_intervals))],lb=0,ub=x_f+15)
x4 = m.Var(value=np.concatenate((np.zeros(int(nr_intervals/2)),np.pi*np.ones(int(nr_intervals)-int(nr_intervals/2)))),lb=-np.pi*2,ub=np.pi*2)
cc= [x*180/3.1415 for x in x4.value]
x5 = m.Var(value=[(x4.value[i+1]-x4.value[i])/step for i in range(0,len(x4.value)-1)]+[(x4[-1]-x4[-2])/step],lb=-step/J_y,ub=step/J_y)
x6 = m.Var(value=np.concatenate((np.linspace(m0,m1,time_burn_node),m1*np.ones(len(t)-time_burn_node))),lb=0,ub=m0)
x2 = m.Var(value=[np.linspace(0,v0,time_burn_node+1)[int(i)] for i in np.linspace(0,time_burn_node,time_burn_node+1)]+[velocity_profile(i) for i in t[1:int(nr_intervals)-time_burn_node]],lb=0,ub=1e3)
x3 = m.Var(value=[np.linspace(np.pi/2,theta_0,time_burn_node+1)[int(i)] for i in np.linspace(0,time_burn_node,time_burn_node+1)]+[np.arctan2(velocity_profile_y(i),velocity_profile_x(i)) for i in t[1:int(nr_intervals)-2*time_burn_node-1]]+[np.linspace(np.arctan2(velocity_profile_y(t[int(nr_intervals)-2*time_burn_node]),velocity_profile_x(t[int(nr_intervals)-2*time_burn_node])),-np.pi/2,time_burn_node+1)[int(i)] for i in np.linspace(0,time_burn_node,time_burn_node+1)],lb=-np.pi*2,ub=np.pi*2)
x0 = m.Var(value=[np.trapz(np.multiply(x2.value[0:i],np.cos(x3.value[0:i])),x=t[0:i]) for i in range(0,int(nr_intervals))],lb=0,ub=x_f+1)
x1 = m.Var(value=[np.trapz(np.multiply(x2.value[0:i],np.sin(x3.value[0:i])),x=t[0:i]) for i in range(0,int(nr_intervals))],lb=0,ub=x_f+15)
x4 = m.Var(value=np.concatenate((np.zeros(int(nr_intervals/2)),np.pi*np.ones(int(nr_intervals)-int(nr_intervals/2)))),lb=-np.pi*2,ub=np.pi*2)
x5 = m.Var(value=[(x4.value[i+1]-x4.value[i])/step for i in range(0,len(x4.value)-1)]+[(x4[-1]-x4[-2])/step],lb=-step/J_y,ub=step/J_y)
x6 = m.Var(value=np.concatenate((np.linspace(m0,m1,time_burn_node),m1*np.ones(len(t)-time_burn_node))),lb=0,ub=m0)
m_dot= m.MV(value=np.concatenate((m_dot_up*np.ones(time_burn_node),np.zeros(len(t)-time_burn_node))),lb=m_dot_low,ub=m_dot_up)
m_dot.value=m_dot.value[0:len(x1.value)]
p = np.zeros(len(x1.value))
p[-1]=1
final = m.Param(value=p,lb=0,ub=1)
m.Equation(x0.dt()==x2*m.cos(x3))
m.Equation(x1.dt()==x2*m.sin(x3))
m.Equation(x6*x2.dt()==(Isp*m_dot*g_0)*m.cos(x4)-x6*g*m.sin(x3))
m.Equation(x6*x2*x3.dt()==x2**2*m.cos(x3)*x6/R_moon+(Isp*m_dot*g_0)*m.sin(x4)-x6*g*m.cos(x3))
m.Equation(x6.dt()==-m_dot)
m.fix_final(x0,x0.value[-1])
m.Equation(x1*final<=1)
m.Minimize((-x6*final))
m.options.MAX_ITER = 1000 # adjust maximum iterations
m.options.SOLVER = 3
m.options.IMODE = 6
m.options.NODES = 3
m.solve()
print(final.value)
print(f"Final Mass: {x6.value[-1]:.3f} s")
这是控制台的输出，显示 objective 函数从负到正，这没有意义，因为最终质量始终为正
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
90r-8.8685190e-001 2.32e+000 9.25e+003 0.9 3.41e+001 - 1.18e-001 2.17e-002f 1
91 -8.8659478e-001 3.03e+000 1.79e+005 0.3 5.68e+001 - 2.78e-001 1.27e-001f 2
92 -8.8283316e-001 6.47e+000 2.86e+007 1.6 1.16e+002 - 9.91e-001 4.76e-002f 1
93 1.8507003e+000 1.87e+001 1.55e+007 1.6 5.71e+003 - 2.41e-001 3.44e-001f 1
94 1.8506996e+000 1.86e+001 1.16e+011 1.6 1.25e+001 10.6 6.17e-003 2.28e-003h 1
95 1.8506996e+000 1.86e+001 1.16e+011 1.6 1.25e+001 10.1 5.75e-002 2.28e-005h 1
我还输出了最终矢量和最终质量的最终值，它们显示了正确的结果：
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0]
Final Mass: 0.993 s
非常感谢您的任何建议
使用 m.open_folder() 检查模型文件并在文本编辑器中查看 gk0_model.apm 文件。
Model
Parameters
p1<= 0.008920571233734825, >= 0.0
p2<= 1, >= 0
End Parameters
Variables
v1<= 1000.0, >= 0
v2<= 6.283185307179586, >= -6.283185307179586
v3<= 4.57104227603633, >= 0
v4<= 18.57104227603633, >= 0
v5<= 6.283185307179586, >= -6.283185307179586
v6<= 2.478718169538892, >= -2.478718169538892
v7<= 1, >= 0
v8<= 1000.0, >= 0
v9<= 6.283185307179586, >= -6.283185307179586
v10<= 4.57104227603633, >= 0
v11<= 18.57104227603633, >= 0
v12<= 6.283185307179586, >= -6.283185307179586
v13<= 2.478718169538892, >= -2.478718169538892
v14<= 1, >= 0
End Variables
Equations
$v10=((v8)*(cos(v9)))
$v11=((v8)*(sin(v9)))
((v14)*($v8))=(((((((46.28)*(p1)))*(9.81)))*(cos(v12)))-((((v14)*(1.62)))*(sin(v9))))
((((v14)*(v8)))*($v9))=((((((((((v8)^(2)))*(cos(v9))))*(v14)))/(1737400.0))+((((((46.28)*(p1)))*(9.81)))*(sin(v12))))-((((v14)*(1.62)))*(cos(v9))))
$v14=(-p1)
((v11)*(p2))<=1
minimize (((-v14))*(p2))
End Equations
Connections
p(end).n(end).v10=4.25390890270255
p(end).n(end).v10=fixed
End Connections
End Model
objective 是 minimize (((-v14))*(p2)) 并且求解器不受限于采用可行的搜索路径达到最优。 v14 的值可能在正值和负值之间波动。
初始化和软化约束有助于找到解决方案。尝试用软约束替换硬终端约束，例如：
#m.fix_final(x0,x0.value[-1])
#m.Equation(x1*final<=1)
m.Minimize(final*(x0-x0.value[-1])**2)
m.Minimize(final*(x1-1)**2)
另一种策略是首先求解一个可行的解决方案，其中决策变量固定在合理的值。设置 m.options.COLDSTART=1 关闭 MV 以找到可行的解决方案。
# -*- coding: utf-8 -*-
from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
# create GEKKO model
m = GEKKO(remote=False)
# m = GEKKO()
print(m.path)
g_0 = 9.81
Isp = 46.28
m_fb=1
m_prop=0.1*m_fb
m_dry=value=m_fb-m_prop
c = 0.8
a = 0.4
g = 1.62
R_moon = 1737.4e3
J_y = 0.06666
m_dot_up=2.5*m_fb*g/(Isp*g_0)
m_dot_low=0*m_fb*g/(Isp*g_0)
theta_0=70*np.pi/180
v0= 3
x_f=np.sin(2*theta_0)*v0*v0/g
v0_x= v0*np.cos(theta_0)
y_max=v0*v0/(2*g)
m0=1
m1=(np.e**(-v0/(Isp*g_0)))
time_burn= m_fb*(m0-m1)/m_dot_up
tf=2*v0*np.sin(theta_0)/g + 2*time_burn
t0=0
nr_intervals=30
step=tf/nr_intervals
t=np.linspace(t0, tf, nr_intervals)
time_burn_node= min(range(len(t)), key=lambda i: abs(t[i]-time_burn))
m.time=t
m0=1
m1=(np.e**(-v0/(Isp*g_0)))
time_burn= m_fb*(m0-m1)/m_dot_up
time_burn_node= min(range(len(t)), key=lambda i: abs(t[i]-time_burn))
y_f=1
v_f_x = 1
v_f_y=3
gamma_f = -90*np.pi/180
alpha_f= 180*np.pi/180
alpha_dot_f=5*np.pi/180
parabola_profile= lambda x: -4*y_max/(x_f*x_f)*x*x+4*y_max*x/x_f
velocity_profile_y=lambda t: (v0*np.sin(theta_0)-g*t)
velocity_profile_x= lambda t: v0*np.cos(theta_0)
velocity_profile=lambda t: np.sqrt(velocity_profile_x(t)*velocity_profile_x(t)+velocity_profile_y(t)*velocity_profile_y(t))
x2 = m.Var(value=[np.linspace(0,v0,time_burn_node+1)[int(i)] for i in np.linspace(0,time_burn_node,time_burn_node+1)]+[velocity_profile(i) for i in t[1:int(nr_intervals)-time_burn_node]],lb=0,ub=1e3)
x3 = m.Var(value=[np.linspace(np.pi/2,theta_0,time_burn_node+1)[int(i)] for i in np.linspace(0,time_burn_node,time_burn_node+1)]+[np.arctan2(velocity_profile_y(i),velocity_profile_x(i)) for i in t[1:int(nr_intervals)-2*time_burn_node-1]]+[np.linspace(np.arctan2(velocity_profile_y(t[int(nr_intervals)-2*time_burn_node]),velocity_profile_x(t[int(nr_intervals)-2*time_burn_node])),-np.pi/2,time_burn_node+1)[int(i)] for i in np.linspace(0,time_burn_node,time_burn_node+1)],lb=-np.pi*2,ub=np.pi*2)
x0 = m.Var(value=[np.trapz(np.multiply(x2.value[0:i],np.cos(x3.value[0:i])),x=t[0:i]) for i in range(0,int(nr_intervals))],lb=0,ub=x_f+1)
x1 = m.Var(value=[np.trapz(np.multiply(x2.value[0:i],np.sin(x3.value[0:i])),x=t[0:i]) for i in range(0,int(nr_intervals))],lb=0,ub=x_f+15)
x4 = m.Var(value=np.concatenate((np.zeros(int(nr_intervals/2)),np.pi*np.ones(int(nr_intervals)-int(nr_intervals/2)))),lb=-np.pi*2,ub=np.pi*2)
cc= [x*180/3.1415 for x in x4.value]
x5 = m.Var(value=[(x4.value[i+1]-x4.value[i])/step for i in range(0,len(x4.value)-1)]+[(x4[-1]-x4[-2])/step],lb=-step/J_y,ub=step/J_y)
x6 = m.Var(value=np.concatenate((np.linspace(m0,m1,time_burn_node),m1*np.ones(len(t)-time_burn_node))),lb=0,ub=m0)
x2 = m.Var(value=[np.linspace(0,v0,time_burn_node+1)[int(i)] for i in np.linspace(0,time_burn_node,time_burn_node+1)]+[velocity_profile(i) for i in t[1:int(nr_intervals)-time_burn_node]],lb=0,ub=1e3)
x3 = m.Var(value=[np.linspace(np.pi/2,theta_0,time_burn_node+1)[int(i)] for i in np.linspace(0,time_burn_node,time_burn_node+1)]+[np.arctan2(velocity_profile_y(i),velocity_profile_x(i)) for i in t[1:int(nr_intervals)-2*time_burn_node-1]]+[np.linspace(np.arctan2(velocity_profile_y(t[int(nr_intervals)-2*time_burn_node]),velocity_profile_x(t[int(nr_intervals)-2*time_burn_node])),-np.pi/2,time_burn_node+1)[int(i)] for i in np.linspace(0,time_burn_node,time_burn_node+1)],lb=-np.pi*2,ub=np.pi*2)
x0 = m.Var(value=[np.trapz(np.multiply(x2.value[0:i],np.cos(x3.value[0:i])),x=t[0:i]) for i in range(0,int(nr_intervals))],lb=0,ub=x_f+1)
x1 = m.Var(value=[np.trapz(np.multiply(x2.value[0:i],np.sin(x3.value[0:i])),x=t[0:i]) for i in range(0,int(nr_intervals))],lb=0,ub=x_f+15)
x4 = m.Var(value=np.concatenate((np.zeros(int(nr_intervals/2)),np.pi*np.ones(int(nr_intervals)-int(nr_intervals/2)))),lb=-np.pi*2,ub=np.pi*2)
x5 = m.Var(value=[(x4.value[i+1]-x4.value[i])/step for i in range(0,len(x4.value)-1)]+[(x4[-1]-x4[-2])/step],lb=-step/J_y,ub=step/J_y)
x6 = m.Var(value=np.concatenate((np.linspace(m0,m1,time_burn_node),m1*np.ones(len(t)-time_burn_node))),lb=0,ub=m0)
m_dot= m.MV(value=np.concatenate((m_dot_up*np.ones(time_burn_node),np.zeros(len(t)-time_burn_node))),lb=m_dot_low,ub=m_dot_up)
m_dot.value=m_dot.value[0:len(x1.value)]
p = np.zeros(len(x1.value))
p[-1]=1
final = m.Param(value=p,lb=0,ub=1)
m.Equation(x0.dt()==x2*m.cos(x3))
m.Equation(x1.dt()==x2*m.sin(x3))
m.Equation(x6*x2.dt()==(Isp*m_dot*g_0)*m.cos(x4)-x6*g*m.sin(x3))
m.Equation(x6*x2*x3.dt()==x2**2*m.cos(x3)*x6/R_moon+(Isp*m_dot*g_0)*m.sin(x4)-x6*g*m.cos(x3))
m.Equation(x6.dt()==-m_dot)
#m.fix_final(x0,x0.value[-1])
#m.Equation(x1*final<=1)
m.Minimize(final*(x0-x0.value[-1])**2)
m.Minimize(final*(x1-1)**2)
m.Minimize((-x6*final))
m.options.MAX_ITER = 2000 # adjust maximum iterations
m.options.SOLVER = 2
m.options.IMODE = 6
m.options.NODES = 3
m.open_folder()
m.options.COLDSTART = 1
m.solve()
print(final.value)
print(f"Final Mass: {x6.value[-1]:.3f} s")
这给出了一个错误：
Number of state variables: 1247
Number of total equations: - 725
Number of slack variables: - 0
---------------------------------------
Degrees of freedom : 522
@error: Degrees of Freedom
* Error: DOF must be zero for this mode
STOPPING...
Traceback (most recent call last):
File "C:\Users\johnh\Desktop\test.py", line 111, in <module>
m.solve()
File "C:\Users\johnh\Python39\lib\site-packages\gekko\gekko.py", line 2185, in solve
raise Exception(response)
Exception: @error: Degrees of Freedom
* Error: DOF must be zero for this mode
STOPPING...
尝试删除任何约束并确定每个变量与约束的配对。如果自由度为零（相同数量的方程和变量），则求解器有时可以找到可以随后优化的初始解。
我正在尝试优化推力推进系统的轨迹。控制变量是质量流率，最终的objective是最大化机器人的质量，最小化使用的推进量。轨迹类似于弹道，具有初始上升阶段和最终下降阶段。
我想我得到了一个很好的初始猜测，但是算法没有收敛。我检查了控制台中的输出，似乎 objective 函数没有正常工作，我认为这就是它没有收敛的原因。
这是我的代码
~# -*- coding: utf-8 -*-
from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
# create GEKKO model
m = GEKKO(remote=False)
# m = GEKKO()
print(m.path)
g_0 = 9.81
Isp = 46.28
m_fb=1
m_prop=0.1*m_fb
m_dry=value=m_fb-m_prop
c = 0.8
a = 0.4
g = 1.62
R_moon = 1737.4e3
J_y = 0.06666
m_dot_up=2.5*m_fb*g/(Isp*g_0)
m_dot_low=0*m_fb*g/(Isp*g_0)
theta_0=70*np.pi/180
v0= 3
x_f=np.sin(2*theta_0)*v0*v0/g
v0_x= v0*np.cos(theta_0)
y_max=v0*v0/(2*g)
m0=1
m1=(np.e**(-v0/(Isp*g_0)))
time_burn= m_fb*(m0-m1)/m_dot_up
tf=2*v0*np.sin(theta_0)/g + 2*time_burn
t0=0
nr_intervals=30
step=tf/nr_intervals
t=np.linspace(t0, tf, nr_intervals)
time_burn_node= min(range(len(t)), key=lambda i: abs(t[i]-time_burn))
m.time=t
m0=1
m1=(np.e**(-v0/(Isp*g_0)))
time_burn= m_fb*(m0-m1)/m_dot_up
time_burn_node= min(range(len(t)), key=lambda i: abs(t[i]-time_burn))
y_f=1
v_f_x = 1
v_f_y=3
gamma_f = -90*np.pi/180
alpha_f= 180*np.pi/180
alpha_dot_f=5*np.pi/180
parabola_profile= lambda x: -4*y_max/(x_f*x_f)*x*x+4*y_max*x/x_f
velocity_profile_y=lambda t: (v0*np.sin(theta_0)-g*t)
velocity_profile_x= lambda t: v0*np.cos(theta_0)
velocity_profile=lambda t: np.sqrt(velocity_profile_x(t)*velocity_profile_x(t)+velocity_profile_y(t)*velocity_profile_y(t))
x2 = m.Var(value=[np.linspace(0,v0,time_burn_node+1)[int(i)] for i in np.linspace(0,time_burn_node,time_burn_node+1)]+[velocity_profile(i) for i in t[1:int(nr_intervals)-time_burn_node]],lb=0,ub=1e3)
x3 = m.Var(value=[np.linspace(np.pi/2,theta_0,time_burn_node+1)[int(i)] for i in np.linspace(0,time_burn_node,time_burn_node+1)]+[np.arctan2(velocity_profile_y(i),velocity_profile_x(i)) for i in t[1:int(nr_intervals)-2*time_burn_node-1]]+[np.linspace(np.arctan2(velocity_profile_y(t[int(nr_intervals)-2*time_burn_node]),velocity_profile_x(t[int(nr_intervals)-2*time_burn_node])),-np.pi/2,time_burn_node+1)[int(i)] for i in np.linspace(0,time_burn_node,time_burn_node+1)],lb=-np.pi*2,ub=np.pi*2)
x0 = m.Var(value=[np.trapz(np.multiply(x2.value[0:i],np.cos(x3.value[0:i])),x=t[0:i]) for i in range(0,int(nr_intervals))],lb=0,ub=x_f+1)
x1 = m.Var(value=[np.trapz(np.multiply(x2.value[0:i],np.sin(x3.value[0:i])),x=t[0:i]) for i in range(0,int(nr_intervals))],lb=0,ub=x_f+15)
x4 = m.Var(value=np.concatenate((np.zeros(int(nr_intervals/2)),np.pi*np.ones(int(nr_intervals)-int(nr_intervals/2)))),lb=-np.pi*2,ub=np.pi*2)
cc= [x*180/3.1415 for x in x4.value]
x5 = m.Var(value=[(x4.value[i+1]-x4.value[i])/step for i in range(0,len(x4.value)-1)]+[(x4[-1]-x4[-2])/step],lb=-step/J_y,ub=step/J_y)
x6 = m.Var(value=np.concatenate((np.linspace(m0,m1,time_burn_node),m1*np.ones(len(t)-time_burn_node))),lb=0,ub=m0)
x2 = m.Var(value=[np.linspace(0,v0,time_burn_node+1)[int(i)] for i in np.linspace(0,time_burn_node,time_burn_node+1)]+[velocity_profile(i) for i in t[1:int(nr_intervals)-time_burn_node]],lb=0,ub=1e3)
x3 = m.Var(value=[np.linspace(np.pi/2,theta_0,time_burn_node+1)[int(i)] for i in np.linspace(0,time_burn_node,time_burn_node+1)]+[np.arctan2(velocity_profile_y(i),velocity_profile_x(i)) for i in t[1:int(nr_intervals)-2*time_burn_node-1]]+[np.linspace(np.arctan2(velocity_profile_y(t[int(nr_intervals)-2*time_burn_node]),velocity_profile_x(t[int(nr_intervals)-2*time_burn_node])),-np.pi/2,time_burn_node+1)[int(i)] for i in np.linspace(0,time_burn_node,time_burn_node+1)],lb=-np.pi*2,ub=np.pi*2)
x0 = m.Var(value=[np.trapz(np.multiply(x2.value[0:i],np.cos(x3.value[0:i])),x=t[0:i]) for i in range(0,int(nr_intervals))],lb=0,ub=x_f+1)
x1 = m.Var(value=[np.trapz(np.multiply(x2.value[0:i],np.sin(x3.value[0:i])),x=t[0:i]) for i in range(0,int(nr_intervals))],lb=0,ub=x_f+15)
x4 = m.Var(value=np.concatenate((np.zeros(int(nr_intervals/2)),np.pi*np.ones(int(nr_intervals)-int(nr_intervals/2)))),lb=-np.pi*2,ub=np.pi*2)
x5 = m.Var(value=[(x4.value[i+1]-x4.value[i])/step for i in range(0,len(x4.value)-1)]+[(x4[-1]-x4[-2])/step],lb=-step/J_y,ub=step/J_y)
x6 = m.Var(value=np.concatenate((np.linspace(m0,m1,time_burn_node),m1*np.ones(len(t)-time_burn_node))),lb=0,ub=m0)
m_dot= m.MV(value=np.concatenate((m_dot_up*np.ones(time_burn_node),np.zeros(len(t)-time_burn_node))),lb=m_dot_low,ub=m_dot_up)
m_dot.value=m_dot.value[0:len(x1.value)]
p = np.zeros(len(x1.value))
p[-1]=1
final = m.Param(value=p,lb=0,ub=1)
m.Equation(x0.dt()==x2*m.cos(x3))
m.Equation(x1.dt()==x2*m.sin(x3))
m.Equation(x6*x2.dt()==(Isp*m_dot*g_0)*m.cos(x4)-x6*g*m.sin(x3))
m.Equation(x6*x2*x3.dt()==x2**2*m.cos(x3)*x6/R_moon+(Isp*m_dot*g_0)*m.sin(x4)-x6*g*m.cos(x3))
m.Equation(x6.dt()==-m_dot)
m.fix_final(x0,x0.value[-1])
m.Equation(x1*final<=1)
m.Minimize((-x6*final))
m.options.MAX_ITER = 1000 # adjust maximum iterations
m.options.SOLVER = 3
m.options.IMODE = 6
m.options.NODES = 3
m.solve()
print(final.value)
print(f"Final Mass: {x6.value[-1]:.3f} s")
这是控制台的输出，显示 objective 函数从负到正，这没有意义，因为最终质量始终为正
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
90r-8.8685190e-001 2.32e+000 9.25e+003 0.9 3.41e+001 - 1.18e-001 2.17e-002f 1
91 -8.8659478e-001 3.03e+000 1.79e+005 0.3 5.68e+001 - 2.78e-001 1.27e-001f 2
92 -8.8283316e-001 6.47e+000 2.86e+007 1.6 1.16e+002 - 9.91e-001 4.76e-002f 1
93 1.8507003e+000 1.87e+001 1.55e+007 1.6 5.71e+003 - 2.41e-001 3.44e-001f 1
94 1.8506996e+000 1.86e+001 1.16e+011 1.6 1.25e+001 10.6 6.17e-003 2.28e-003h 1
95 1.8506996e+000 1.86e+001 1.16e+011 1.6 1.25e+001 10.1 5.75e-002 2.28e-005h 1
我还输出了最终矢量和最终质量的最终值，它们显示了正确的结果：
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0]
Final Mass: 0.993 s
非常感谢您的任何建议
使用 m.open_folder() 检查模型文件并在文本编辑器中查看 gk0_model.apm 文件。
Model
Parameters
p1<= 0.008920571233734825, >= 0.0
p2<= 1, >= 0
End Parameters
Variables
v1<= 1000.0, >= 0
v2<= 6.283185307179586, >= -6.283185307179586
v3<= 4.57104227603633, >= 0
v4<= 18.57104227603633, >= 0
v5<= 6.283185307179586, >= -6.283185307179586
v6<= 2.478718169538892, >= -2.478718169538892
v7<= 1, >= 0
v8<= 1000.0, >= 0
v9<= 6.283185307179586, >= -6.283185307179586
v10<= 4.57104227603633, >= 0
v11<= 18.57104227603633, >= 0
v12<= 6.283185307179586, >= -6.283185307179586
v13<= 2.478718169538892, >= -2.478718169538892
v14<= 1, >= 0
End Variables
Equations
$v10=((v8)*(cos(v9)))
$v11=((v8)*(sin(v9)))
((v14)*($v8))=(((((((46.28)*(p1)))*(9.81)))*(cos(v12)))-((((v14)*(1.62)))*(sin(v9))))
((((v14)*(v8)))*($v9))=((((((((((v8)^(2)))*(cos(v9))))*(v14)))/(1737400.0))+((((((46.28)*(p1)))*(9.81)))*(sin(v12))))-((((v14)*(1.62)))*(cos(v9))))
$v14=(-p1)
((v11)*(p2))<=1
minimize (((-v14))*(p2))
End Equations
Connections
p(end).n(end).v10=4.25390890270255
p(end).n(end).v10=fixed
End Connections
End Model
objective 是 minimize (((-v14))*(p2)) 并且求解器不受限于采用可行的搜索路径达到最优。 v14 的值可能在正值和负值之间波动。
初始化和软化约束有助于找到解决方案。尝试用软约束替换硬终端约束，例如：
#m.fix_final(x0,x0.value[-1])
#m.Equation(x1*final<=1)
m.Minimize(final*(x0-x0.value[-1])**2)
m.Minimize(final*(x1-1)**2)
另一种策略是首先求解一个可行的解决方案，其中决策变量固定在合理的值。设置 m.options.COLDSTART=1 关闭 MV 以找到可行的解决方案。
# -*- coding: utf-8 -*-
from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
# create GEKKO model
m = GEKKO(remote=False)
# m = GEKKO()
print(m.path)
g_0 = 9.81
Isp = 46.28
m_fb=1
m_prop=0.1*m_fb
m_dry=value=m_fb-m_prop
c = 0.8
a = 0.4
g = 1.62
R_moon = 1737.4e3
J_y = 0.06666
m_dot_up=2.5*m_fb*g/(Isp*g_0)
m_dot_low=0*m_fb*g/(Isp*g_0)
theta_0=70*np.pi/180
v0= 3
x_f=np.sin(2*theta_0)*v0*v0/g
v0_x= v0*np.cos(theta_0)
y_max=v0*v0/(2*g)
m0=1
m1=(np.e**(-v0/(Isp*g_0)))
time_burn= m_fb*(m0-m1)/m_dot_up
tf=2*v0*np.sin(theta_0)/g + 2*time_burn
t0=0
nr_intervals=30
step=tf/nr_intervals
t=np.linspace(t0, tf, nr_intervals)
time_burn_node= min(range(len(t)), key=lambda i: abs(t[i]-time_burn))
m.time=t
m0=1
m1=(np.e**(-v0/(Isp*g_0)))
time_burn= m_fb*(m0-m1)/m_dot_up
time_burn_node= min(range(len(t)), key=lambda i: abs(t[i]-time_burn))
y_f=1
v_f_x = 1
v_f_y=3
gamma_f = -90*np.pi/180
alpha_f= 180*np.pi/180
alpha_dot_f=5*np.pi/180
parabola_profile= lambda x: -4*y_max/(x_f*x_f)*x*x+4*y_max*x/x_f
velocity_profile_y=lambda t: (v0*np.sin(theta_0)-g*t)
velocity_profile_x= lambda t: v0*np.cos(theta_0)
velocity_profile=lambda t: np.sqrt(velocity_profile_x(t)*velocity_profile_x(t)+velocity_profile_y(t)*velocity_profile_y(t))
x2 = m.Var(value=[np.linspace(0,v0,time_burn_node+1)[int(i)] for i in np.linspace(0,time_burn_node,time_burn_node+1)]+[velocity_profile(i) for i in t[1:int(nr_intervals)-time_burn_node]],lb=0,ub=1e3)
x3 = m.Var(value=[np.linspace(np.pi/2,theta_0,time_burn_node+1)[int(i)] for i in np.linspace(0,time_burn_node,time_burn_node+1)]+[np.arctan2(velocity_profile_y(i),velocity_profile_x(i)) for i in t[1:int(nr_intervals)-2*time_burn_node-1]]+[np.linspace(np.arctan2(velocity_profile_y(t[int(nr_intervals)-2*time_burn_node]),velocity_profile_x(t[int(nr_intervals)-2*time_burn_node])),-np.pi/2,time_burn_node+1)[int(i)] for i in np.linspace(0,time_burn_node,time_burn_node+1)],lb=-np.pi*2,ub=np.pi*2)
x0 = m.Var(value=[np.trapz(np.multiply(x2.value[0:i],np.cos(x3.value[0:i])),x=t[0:i]) for i in range(0,int(nr_intervals))],lb=0,ub=x_f+1)
x1 = m.Var(value=[np.trapz(np.multiply(x2.value[0:i],np.sin(x3.value[0:i])),x=t[0:i]) for i in range(0,int(nr_intervals))],lb=0,ub=x_f+15)
x4 = m.Var(value=np.concatenate((np.zeros(int(nr_intervals/2)),np.pi*np.ones(int(nr_intervals)-int(nr_intervals/2)))),lb=-np.pi*2,ub=np.pi*2)
cc= [x*180/3.1415 for x in x4.value]
x5 = m.Var(value=[(x4.value[i+1]-x4.value[i])/step for i in range(0,len(x4.value)-1)]+[(x4[-1]-x4[-2])/step],lb=-step/J_y,ub=step/J_y)
x6 = m.Var(value=np.concatenate((np.linspace(m0,m1,time_burn_node),m1*np.ones(len(t)-time_burn_node))),lb=0,ub=m0)
x2 = m.Var(value=[np.linspace(0,v0,time_burn_node+1)[int(i)] for i in np.linspace(0,time_burn_node,time_burn_node+1)]+[velocity_profile(i) for i in t[1:int(nr_intervals)-time_burn_node]],lb=0,ub=1e3)
x3 = m.Var(value=[np.linspace(np.pi/2,theta_0,time_burn_node+1)[int(i)] for i in np.linspace(0,time_burn_node,time_burn_node+1)]+[np.arctan2(velocity_profile_y(i),velocity_profile_x(i)) for i in t[1:int(nr_intervals)-2*time_burn_node-1]]+[np.linspace(np.arctan2(velocity_profile_y(t[int(nr_intervals)-2*time_burn_node]),velocity_profile_x(t[int(nr_intervals)-2*time_burn_node])),-np.pi/2,time_burn_node+1)[int(i)] for i in np.linspace(0,time_burn_node,time_burn_node+1)],lb=-np.pi*2,ub=np.pi*2)
x0 = m.Var(value=[np.trapz(np.multiply(x2.value[0:i],np.cos(x3.value[0:i])),x=t[0:i]) for i in range(0,int(nr_intervals))],lb=0,ub=x_f+1)
x1 = m.Var(value=[np.trapz(np.multiply(x2.value[0:i],np.sin(x3.value[0:i])),x=t[0:i]) for i in range(0,int(nr_intervals))],lb=0,ub=x_f+15)
x4 = m.Var(value=np.concatenate((np.zeros(int(nr_intervals/2)),np.pi*np.ones(int(nr_intervals)-int(nr_intervals/2)))),lb=-np.pi*2,ub=np.pi*2)
x5 = m.Var(value=[(x4.value[i+1]-x4.value[i])/step for i in range(0,len(x4.value)-1)]+[(x4[-1]-x4[-2])/step],lb=-step/J_y,ub=step/J_y)
x6 = m.Var(value=np.concatenate((np.linspace(m0,m1,time_burn_node),m1*np.ones(len(t)-time_burn_node))),lb=0,ub=m0)
m_dot= m.MV(value=np.concatenate((m_dot_up*np.ones(time_burn_node),np.zeros(len(t)-time_burn_node))),lb=m_dot_low,ub=m_dot_up)
m_dot.value=m_dot.value[0:len(x1.value)]
p = np.zeros(len(x1.value))
p[-1]=1
final = m.Param(value=p,lb=0,ub=1)
m.Equation(x0.dt()==x2*m.cos(x3))
m.Equation(x1.dt()==x2*m.sin(x3))
m.Equation(x6*x2.dt()==(Isp*m_dot*g_0)*m.cos(x4)-x6*g*m.sin(x3))
m.Equation(x6*x2*x3.dt()==x2**2*m.cos(x3)*x6/R_moon+(Isp*m_dot*g_0)*m.sin(x4)-x6*g*m.cos(x3))
m.Equation(x6.dt()==-m_dot)
#m.fix_final(x0,x0.value[-1])
#m.Equation(x1*final<=1)
m.Minimize(final*(x0-x0.value[-1])**2)
m.Minimize(final*(x1-1)**2)
m.Minimize((-x6*final))
m.options.MAX_ITER = 2000 # adjust maximum iterations
m.options.SOLVER = 2
m.options.IMODE = 6
m.options.NODES = 3
m.open_folder()
m.options.COLDSTART = 1
m.solve()
print(final.value)
print(f"Final Mass: {x6.value[-1]:.3f} s")
这给出了一个错误：
Number of state variables: 1247
Number of total equations: - 725
Number of slack variables: - 0
---------------------------------------
Degrees of freedom : 522
@error: Degrees of Freedom
* Error: DOF must be zero for this mode
STOPPING...
Traceback (most recent call last):
File "C:\Users\johnh\Desktop\test.py", line 111, in <module>
m.solve()
File "C:\Users\johnh\Python39\lib\site-packages\gekko\gekko.py", line 2185, in solve
raise Exception(response)
Exception: @error: Degrees of Freedom
* Error: DOF must be zero for this mode
STOPPING...
尝试删除任何约束并确定每个变量与约束的配对。如果自由度为零（相同数量的方程和变量），则求解器有时可以找到可以随后优化的初始解。