<h3>首先,一些警告:</h3>
<ul>
<li>数值优化很难做对</li>
<li><a href="https://www.coin-or.org/Ipopt/" rel="nofollow noreferrer">ipopt</a>是非常复杂的软件
<ul>
<li>将ipopt与数值微分相结合听起来像是在自找麻烦,但这当然取决于你的问题</li>
<li>ipopt几乎总是基于<a href="https://en.wikipedia.org/wiki/Automatic_differentiation" rel="nofollow noreferrer">automatic-differentiation tools</a>,而不是<a href="https://en.wikipedia.org/wiki/Numerical_differentiation" rel="nofollow noreferrer">numerical-differentiation</a>!在</li>
</ul></li>
</ul>
<h3>还有更多:</h3>
<ul>
<li>由于这是一项复杂的任务,而且python+ipopt的状态没有其他语言中的好(例如<a href="https://julialang.org/" rel="nofollow noreferrer">julia</a>+<a href="https://jump.readthedocs.io/en/latest/" rel="nofollow noreferrer">JuMP</a>),所以需要做一些工作</li>
</ul>
<h3>还有一些替代品:</h3>
<ul>
<li>使用<a href="http://www.pyomo.org/" rel="nofollow noreferrer">pyomo</a>,它包装了ipopt并具有自动区分功能</li>
<li>使用<a href="https://github.com/casadi/casadi" rel="nofollow noreferrer">casadi</a>,它也包装ipopt并具有自动区分功能</li>
<li>使用<a href="https://github.com/HIPS/autograd" rel="nofollow noreferrer">autograd</a>自动计算numpy代码子集上的渐变
<ul>
<li>然后使用cyipot添加这些</li>
</ul></li>
<li><a href="https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html" rel="nofollow noreferrer">scipy.minimize with solvers SLSQP or COBYLA</a>,<em>可以为您做任何事情</em>(SLSQP可以使用等式和不等式约束;只有COBYLA的不等式约束,其中通过<code>x >= y</code>+<code>x <= y</code><em>可以</em>来模拟等式约束)</li>
</ul>
<h3>使用工具完成任务</h3>
<p>完整的示例问题在<a href="http://www.ai7.uni-bayreuth.de/test_problem_coll.pdf" rel="nofollow noreferrer">Test Examples for Nonlinear Programming Codes</a>中定义:</p>
<p><a href="https://i.imgur.com/8ogR0gE.png" rel="nofollow noreferrer"><img src="https://i.imgur.com/8ogR0gE.png" alt="enter image description here"/></a></p>
<p>这里有一些代码,基于数值微分,解决你的测试问题,包括官方设置(函数,梯度,起点,边界…)</p>
<pre><code>import numpy as np
import scipy.sparse as sps
import ipopt
from scipy.optimize import approx_fprime
class Problem40(object):
""" # Hock & Schittkowski test problem #40
Basic structure follows:
- cyipopt example from https://pythonhosted.org/ipopt/tutorial.html#defining-the-problem
- which follows ipopt's docs from: https://www.coin-or.org/Ipopt/documentation/node22.html
Changes:
- numerical-diff using scipy for function & constraints
- removal of hessian-calculation
- we will use limited-memory approximation
- ipopt docs: https://www.coin-or.org/Ipopt/documentation/node31.html
- (because i'm too lazy to reason about the math; lagrange and co.)
"""
def __init__(self):
self.num_diff_eps = 1e-8 # maybe tuning needed!
def objective(self, x):
# callback for objective
return -np.prod(x) # -x1 x2 x3 x4
def constraint_0(self, x):
return np.array([x[0]**3 + x[1]**2 -1])
def constraint_1(self, x):
return np.array([x[0]**2 * x[3] - x[2]])
def constraint_2(self, x):
return np.array([x[3]**2 - x[1]])
def constraints(self, x):
# callback for constraints
return np.concatenate([self.constraint_0(x),
self.constraint_1(x),
self.constraint_2(x)])
def gradient(self, x):
# callback for gradient
return approx_fprime(x, self.objective, self.num_diff_eps)
def jacobian(self, x):
# callback for jacobian
return np.concatenate([
approx_fprime(x, self.constraint_0, self.num_diff_eps),
approx_fprime(x, self.constraint_1, self.num_diff_eps),
approx_fprime(x, self.constraint_2, self.num_diff_eps)])
def hessian(self, x, lagrange, obj_factor):
return False # we will use quasi-newton approaches to use hessian-info
# progress callback
def intermediate(
self,
alg_mod,
iter_count,
obj_value,
inf_pr,
inf_du,
mu,
d_norm,
regularization_size,
alpha_du,
alpha_pr,
ls_trials
):
print("Objective value at iteration #%d is - %g" % (iter_count, obj_value))
# Remaining problem definition; still following official source:
# http://www.ai7.uni-bayreuth.de/test_problem_coll.pdf
# start-point -> infeasible
x0 = [0.8, 0.8, 0.8, 0.8]
# variable-bounds -> empty => np.inf-approach deviates from cyipopt docs!
lb = [-np.inf, -np.inf, -np.inf, -np.inf]
ub = [np.inf, np.inf, np.inf, np.inf]
# constraint bounds -> c == 0 needed -> both bounds = 0
cl = [0, 0, 0]
cu = [0, 0, 0]
nlp = ipopt.problem(
n=len(x0),
m=len(cl),
problem_obj=Problem40(),
lb=lb,
ub=ub,
cl=cl,
cu=cu
)
# IMPORTANT: need to use limited-memory / lbfgs here as we didn't give a valid hessian-callback
nlp.addOption(b'hessian_approximation', b'limited-memory')
x, info = nlp.solve(x0)
print(x)
print(info)
# CORRECT RESULT & SUCCESSFUL STATE
</code></pre>
<p>输出:</p>
^{pr2}$
<h3>关于代码</h3>
<ul>
<li>我们使用<a href="https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.approx_fprime.html" rel="nofollow noreferrer">scipy's approx_fprime</a>,它基本上是为<a href="https://docs.scipy.org/doc/scipy/reference/optimize.html" rel="nofollow noreferrer">scipy.optimize</a>中所有基于梯度的优化器添加的</li>
<li>正如源代码中所述;我没有考虑ipopt对hessian的需求,我们使用ipopts <a href="https://www.coin-or.org/Ipopt/documentation/node31.html" rel="nofollow noreferrer">hessian-approximation</a>
<ul>
<li>基本思想在<a href="https://en.wikipedia.org/wiki/Limited-memory_BFGS" rel="nofollow noreferrer">wiki: LBFGS</a>中描述</li>
</ul></li>
<li>我忽略了ipopts对约束Jacobian的稀疏结构的需要<em>
<ul>
<li>一个默认假设:<em>默认的hessian结构是一个较低的三角形矩阵,使用了</em>,我不会对这里可能发生的情况给出任何保证(性能差vs.破坏一切)</li>
</ul></li>
</ul>