如何调试 LSQ 子问题中的失败

Question

我正在使用“SLSQP”和 SCIPY Optimizer 驱动程序解决 dymos 问题，由于“LSQ 子问题中超过 3^n 次迭代”，我的问题失败了。我一直在尝试了解此错误的确切含义，以确定我应该在模型中的何处开始调试，但我没有找到有关该主题的任何有用资源。你们都提供了关于如何调试其他优化器错误的有用反馈，所以我想也许你们对这个错误也有一些建议。

请注意，我的变量缩放比例非常好。另外，如果我关闭优化变量opt=False，让优化器收敛轨迹，优化就完成了。然后，如果我启用优化变量 opt=True 并热启动问题，优化也会收敛。虽然这可行，但我必须解决许多不同的问题，并且为每个问题找到一个热启动解决方案是不现实的。我还觉得奇怪的是 LSQ 子问题如此依赖于初始猜测。

Answer 1

没有更多细节，很难说得很具体。然而，一般优化和基于伪谱方法的最优控制优化对初始值高度敏感的情况并不少见。以下不是作为解决方案提出的，甚至也不是针对所有 dymos 案例的一般推荐最佳实践。但是对于您发现需要更好的初始猜测的问题，这是一个很好的起点：

import openmdao.api as om
import dymos as dm
from dymos.examples.plotting import plot_results
from dymos.examples.brachistochrone import BrachistochroneODE
import matplotlib.pyplot as plt

#
# Initialize the Problem and the optimization driver
#
p = om.Problem(model=om.Group())
p.driver = om.ScipyOptimizeDriver()
# p.driver.options['tol'] = 1e-9
p.driver.declare_coloring()

#
# Create a trajectory and add a phase to it
#
traj = p.model.add_subsystem('traj', dm.Trajectory())

phase = traj.add_phase('phase0',
                       dm.Phase(ode_class=BrachistochroneODE,
                                transcription=dm.GaussLobatto(num_segments=10)))

#
# Set the variables
#
phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10))

phase.add_state('x', fix_initial=True, fix_final=True)

phase.add_state('y', fix_initial=True, fix_final=True)

phase.add_state('v', fix_initial=True, fix_final=False)

phase.add_control('theta', continuity=True, rate_continuity=True,
                  units='deg', lower=0.01, upper=179.9)

phase.add_parameter('g', units='m/s**2', val=9.80665)

#
# Minimize time at the end of the phase
#
phase.add_objective('time', loc='final', scaler=10)

p.model.linear_solver = om.DirectSolver()

#
# Setup the Problem
#
p.setup()

# intial guesses for times and controls are important before running a sim
p['traj.phase0.t_initial'] = 0.0
p['traj.phase0.t_duration'] = 2.0
p.set_val('traj.phase0.controls:theta', phase.interp('theta', ys=[5, 100.5]))


# need this to set initial conditions for x,y,v so simulate starts from the right condition
p.set_val('traj.phase0.states:x', phase.interp('x', ys=[0, 10]))
p.set_val('traj.phase0.states:y', phase.interp('y', ys=[10, 5]))
p.set_val('traj.phase0.states:v', phase.interp('v', ys=[0, 9.9]))

# Run a simulation to get physically valid initial guesses
p.run_model()
exp_out = traj.simulate()

# extract the initial values out of the sim result
x_sim = exp_out.get_val('traj.phase0.timeseries.states:x')
y_sim = exp_out.get_val('traj.phase0.timeseries.states:y')
v_sim = exp_out.get_val('traj.phase0.timeseries.states:v')
t_sim = exp_out.get_val('traj.phase0.timeseries.time')

# need to manually reset the last value to the correct BCs for states where fix_final=True
# but the initial values will be correct because we set them above before sim
x_sim[-1] = 10
y_sim[-1] = 5



# set initial guesses based on initial sim into the model
# (Dymos will re-interpolate them onto the collocation grid)
# comment/uncomment these three lines to see the effect of the initial guess
p.set_val('traj.phase0.states:x', phase.interp('x', xs=t_sim, ys=x_sim))
p.set_val('traj.phase0.states:y', phase.interp('y', xs=t_sim, ys=y_sim))
p.set_val('traj.phase0.states:v', phase.interp('v', xs=t_sim, ys=v_sim))

#
# Solve for the optimal trajectory
#
dm.run_problem(p)

# Check the results
print(p.get_val('traj.phase0.timeseries.time')[-1])

在没有初步猜测的情况下运行时，我得到：

Full total jacobian was computed 3 times, taking 0.019553 seconds.
Total jacobian shape: (40, 50)


Jacobian shape: (40, 50)  (19.95% nonzero)
FWD solves: 13   REV solves: 0
Total colors vs. total size: 13 vs 50  (74.0% improvement)

Sparsity computed using tolerance: 1e-25
Time to compute sparsity: 0.019553 sec.
Time to compute coloring: 0.028856 sec.
Memory to compute coloring: 0.000000 MB.
Optimization terminated successfully    (Exit mode 0)
            Current function value: [18.0161673]
            Iterations: 24
            Function evaluations: 24
            Gradient evaluations: 24
Optimization Complete
-----------------------------------
[1.80161673]

当以初始猜测运行时，我得到：

Full total jacobian was computed 3 times, taking 0.020724 seconds.
Total jacobian shape: (40, 50)


Jacobian shape: (40, 50)  (19.95% nonzero)
FWD solves: 13   REV solves: 0
Total colors vs. total size: 13 vs 50  (74.0% improvement)

Sparsity computed using tolerance: 1e-25
Time to compute sparsity: 0.020724 sec.
Time to compute coloring: 0.029557 sec.
Memory to compute coloring: 0.000000 MB.
Optimization terminated successfully    (Exit mode 0)
            Current function value: [18.01657396]
            Iterations: 11
            Function evaluations: 11
            Gradient evaluations: 11
Optimization Complete
-----------------------------------
[1.8016574]

所以你可以看到这里的初始猜测有助于更快地收敛。我想再次强调，情况并非总是如此。如果我将最终时间的初始猜测（作为优化的一部分解决---因为目标是最小化旅行时间！）从 2 更改为 10，那么初始模拟是最终解决方案的更差的近似值.

p['traj.phase0.t_duration'] = 10.0

然后我得到这个标准的初始猜测：

Optimization terminated successfully    (Exit mode 0)
            Current function value: [18.01612053]
            Iterations: 39
            Function evaluations: 39
            Gradient evaluations: 39
Optimization Complete
-----------------------------------
[1.80161205]

这是模拟的初始猜测：

Optimization terminated successfully    (Exit mode 0)
            Current function value: [18.01677452]
            Iterations: 43
            Function evaluations: 44
            Gradient evaluations: 43
Optimization Complete
-----------------------------------
[1.80167745]

所以这次模拟的猜测又进行了几次迭代（虽然还是找到了正确答案）。您可以直观地看到为什么会出现这种情况，因为我将我对持续时间的初始估计从 2 秒（接近最佳 1.801）更改为 10 秒（远离最佳值），因此初始模拟的估计值要差得多最优解。 因此，更一般地说，您为非状态值（即时间、控件）提供合理起始值的能力将影响您从模拟中得到的初始猜测的好坏。

即使使用“错误”的模拟初始猜测，您仍然可以提高优化的稳定性。

还有一点要注意：虽然出于教学目的，我在这里手动输入了猜测，但 dymos 确实有一个辅助函数，可以从案例数据库中重新加载猜测。请参阅 run_problem method 的 restart 参数。