我编写了一个不收敛的非线性轨道推力和角度优化问题。下面是 Python 代码，下面是日志的片段。基本上，每个时间步长都有一个推力和角度控制变量，需要对其进行优化以最小化燃料消耗，该变量由推力随时间的总和定义。然而，这并没有按预期收敛。我尝试过使用 scipy.optimise.minimize 但这似乎也无法最小化目标函数。是不是我设置错了？为什么它不能最小化目标函数？我的限制和界限看起来很合理。
from gekko import GEKKO
import numpy as np
# Initialize Model
m = GEKKO(remote=True)
m.time = np.linspace(0,140,801)
# State Variables
r = m.Var(value = 1.0)
theta = m.Var(value = 0.0)
vr = m.Var(value = 0.0)
vt = m.Var(value = 1.0)
# Control Variables
thrust = m.MV(value = 0.005, lb = 0.0, ub = 0.01)
thrust.STATUS = 1
angle = m.MV(value = 0.0, lb = -np.pi/2, ub = np.pi/2)
angle.STATUS = 1
# Optimise fuel consumption
fuel_consumption = m.Var(value = 0.0)
# Dynamics
m.Equation(r.dt() == vr)
m.Equation(theta.dt() == vt/r)
m.Equation(vr.dt() == (vt**2)/r - 1/(r**2) + thrust*m.sin(angle))
m.Equation(vt.dt() == -(vr*vt)/r + thrust*m.cos(angle))
m.Equation(fuel_consumption.dt() == thrust) # Accumulate thrust over time
# Boundary Conditions
m.fix(r, pos=len(m.time)-1,val = 4.0)
m.fix(theta, pos=len(m.time)-1,val = 0.0)
m.fix(vr, pos=len(m.time)-1,val = 0)
m.fix(vt, pos=len(m.time)-1,val = 0.5)
# Objective function
m.Obj(fuel_consumption)
#Set global options
m.options.IMODE = 6 # Dynamic Optimization (Simultaneous ODEs)
m.options.NODES = 3 # Collocation nodes
m.options.SOLVER = 3 # IPOPT solver
#Solve simulation
m.solve(disp=True) # solve on public server
#Results
print('')
print('Results')
print('Objective: ',m.options.OBJFCNVAL)
print('Solution: ', thrust.value)
日志是：
----------------------------------------------------------------
APMonitor, Version 1.0.1
APMonitor Optimization Suite
----------------------------------------------------------------
--------- APM Model Size ------------
Each time step contains
Objects : 0
Constants : 0
Variables : 7
Intermediates: 0
Connections : 8
Equations : 6
Residuals : 6
Number of state variables: 25592
Number of total equations: - 24000
Number of slack variables: - 0
---------------------------------------
Degrees of freedom : 1592
**********************************************
Dynamic Control with Interior Point Solver
**********************************************
Info: Exact Hessian
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit http://projects.coin-or.org/Ipopt
******************************************************************************
This is Ipopt version 3.12.10, running with linear solver ma57.
Number of nonzeros in equality constraint Jacobian...: 67966
Number of nonzeros in inequality constraint Jacobian.: 6400
Number of nonzeros in Lagrangian Hessian.............: 11195
Total number of variables............................: 25592
variables with only lower bounds: 3200
variables with lower and upper bounds: 3200
variables with only upper bounds: 0
Total number of equality constraints.................: 20800
Total number of inequality constraints...............: 3200
inequality constraints with only lower bounds: 3200
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 1.3103987e-02 3.00e+00 1.00e+00 0.0 0.00e+00 - 0.00e+00 0.00e+00 0
Reallocating memory for MA57: lfact (641403)
1 3.8902526e+01 2.77e+00 4.18e+00 -1.3 1.47e+02 - 2.34e-01 7.51e-02f 1
......
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
250r 2.8607519e+02 1.73e+02 9.99e+02 2.2 0.00e+00 3.1 0.00e+00 2.58e-07R 5
Number of Iterations....: 250
(scaled) (unscaled)
Objective...............: 2.8607519277584993e+02 2.8607519277584993e+02
Dual infeasibility......: 2.1703502621114161e+00 2.1703502621114161e+00
Constraint violation....: 1.7251865116637478e+02 1.7251865116637478e+02
Complementarity.........: 3.1072980675208779e+00 3.1072980675208779e+00
Overall NLP error.......: 1.7251865116637478e+02 1.7251865116637478e+02
Number of objective function evaluations = 323
Number of objective gradient evaluations = 252
Number of equality constraint evaluations = 323
Number of inequality constraint evaluations = 323
Number of equality constraint Jacobian evaluations = 254
Number of inequality constraint Jacobian evaluations = 254
Number of Lagrangian Hessian evaluations = 250
Total CPU secs in IPOPT (w/o function evaluations) = 26.390
Total CPU secs in NLP function evaluations = 45.994
EXIT: Maximum Number of Iterations Exceeded.
我尝试过更改时间尺度，并使用 scipy minimize，但似乎无法使其最小化。
尝试通过软化边界条件分两个阶段解决问题：
# Soft Boundary Conditions
p = np.zeros_like(m.time); p[-1]=1
final = m.Param(p)
m.Minimize(final*(r-4)**2)
m.Minimize(final*(theta-0)**2)
m.Minimize(final*(vr-0)**2)
m.Minimize(final*(vt-0.5)**2)
第一个优化是收敛方程的模拟：
Number of state variables: 24000
Number of total equations: - 24000
Number of slack variables: - 0
---------------------------------------
Degrees of freedom : 0
第二次优化使方程收敛，并使用
angle.STATUS=1thrust.STATUS=1 Number of state variables: 25592
Number of total equations: - 24000
Number of slack variables: - 0
---------------------------------------
Degrees of freedom : 1592
这给出了一个成功的解决方案：
from gekko import GEKKO
import numpy as np
# Initialize Model
m = GEKKO(remote=True)
m.time = np.linspace(0,140,801)
# State Variables
r = m.Var(value = 1.0)
theta = m.Var(value = 0.0)
vr = m.Var(value = 0.0)
vt = m.Var(value = 1.0)
# Control Variables
thrust = m.MV(value = 0.005, lb = 0.0, ub = 0.01)
thrust.STATUS = 0
angle = m.MV(value = 0.0, lb = -np.pi/4, ub = np.pi/4)
angle.STATUS = 0
# Optimise fuel consumption
fuel_consumption = m.Var(value = 0.0)
# Dynamics
m.Equation(r.dt() == vr)
m.Equation(theta.dt() == vt/r)
m.Equation(vr.dt() == (vt**2)/r - 1/(r**2) + thrust*m.sin(angle))
m.Equation(vt.dt() == -(vr*vt)/r + thrust*m.cos(angle))
m.Equation(fuel_consumption.dt() == thrust) # Accumulate thrust over time
# Soft Boundary Conditions
p = np.zeros_like(m.time); p[-1]=1
final = m.Param(p)
m.Minimize(final*(r-4)**2)
m.Minimize(final*(theta-0)**2)
m.Minimize(final*(vr-0)**2)
m.Minimize(final*(vt-0.5)**2)
# Objective function
m.Obj(fuel_consumption)
#Set global options
m.options.IMODE = 6 # Dynamic Optimization (Simultaneous ODEs)
m.options.NODES = 3 # Collocation nodes
m.options.SOLVER = 3 # IPOPT solver
#Solve simulation
m.solve(disp=True) # solve on public server
thrust.STATUS = 1
angle.STATUS = 1
# Hard Boundary Conditions
m.fix(theta, pos=len(m.time)-1,val = 0.0)
m.fix(vr, pos=len(m.time)-1,val = 0)
m.fix(r, pos=len(m.time)-1,val = 4.0)
m.fix(vt, pos=len(m.time)-1,val = 0.5)
m.solve(disp=True)
#Results
print('')
print('Results')
print('Objective: ',m.options.OBJFCNVAL)
print('Solution: ', thrust.value)
两阶段求解的主要思想是帮助求解器在优化之前找到初始可行解。硬终端约束也可能不可行。推荐的方法是尝试软终端约束（基于目标，可能会被违反）以利于获得解决方案。