【问题标题】:Lagrange multipliers numerical approach拉格朗日乘数数值方法
【发布时间】:2023-10-20 09:29:02
【问题描述】:

我正在尝试优化以下受约束的功能:

def damage(a, e, cr, cd):
    return 100*(1+a)*(1+e)*(1+cr*cd)

def constraint(a, e, cr, cd):
    return (100/0.466)*(a+e)+(100/0.622)*(2*cr+cd)

当手动求解拉格朗日时,我得到以下输出:

import numpy as np
import sympy as smp
a, e, c, d, l = smp.symbols('a e c d l')
eq1 = smp.Eq(1/(1+a), (100/46.6)*l)
eq2 = smp.Eq(1/(1+e), (100/46.6)*l)
eq3 = smp.Eq(d/(1+c*d), (100/62.2)*2*l)
eq4 = smp.Eq(c/(1+c*d), (100/62.2)*l)
eq5 = smp.Eq((100/46.6)*(a+e)+(100/62.2)*(2*c+d) - 300, 0)

solution = np.array(smp.solve([eq1, eq2, eq3, eq4, eq5], [a, e, c, d, l]))
print(solution[0]/100)
print('Constraint', '{:,.0f}'.format(constraint(*(solution/100)[0][:-1])))
print('Max damage', '{:,.0f}'.format(float(round(damage(*(solution/100)[0][:-1])))))

[0.344658405015485 0.344658405015485 0.236481193219279 0.472962386438559 0.000131394038153319] 约束 300 最大伤害201

为了能够通过数值方法解决这个问题,我修改了问题的公式,明确地分别说明了约束(将主要约束分离成更小的约束)。我明确说明了变量之间所需的关系,并且只约束了其中一个变量,然后确定了所有其他变量的状态。

# We first convert this into a minimization problem.
from scipy import optimize
def damage_min(x):
    return -100*(1+x[0])*(1+x[1])*(1+x[2]*x[3])

# next we define the constrains (equal to 0)
def constraints(x):
    c1 = x[0] - x[1]
    c2 = 2*x[2] - x[3]
    c3 = x[0]/x[3] - 0.466/0.622
    c4 = x[3] - 0.466
    return np.array([c1, c2, c3, c4])
cons = ({'type': 'eq',
         'fun' : constraints})

# We solve the minimization problem
x_initial = np.array([34.4658405015485, 34.4658405015485, 23.6481193219279, 47.2962386438559])
solution = optimize.minimize(damage_min, x_initial, constraints=cons)
print(solution.x)
print('Constraint',  '{:,.0f}'.format(constraint(*(solution.x))))
print('Max damage', '{:,.0f}'.format(float(round(damage(*(solution.x))))))

[0.3491254 0.3491254 0.233 0.466] 约束 300 最大伤害202

我的问题如下。如何通过数值优化单个函数(例如拉格朗日乘数)重新创建上述最佳结果?当我尝试将这两个函数放在一个函数中时,我得到了这个输出。

const = 300
def lagrangian(a, e, cr, cd, lam):
    return -damage(a, e, cr, cd) + lam*(round(constraint(a, e, cr, cd)) - const)

def vector_lagrangian(x):
    return lagrangian(x[0], x[1], x[2], x[3], x[4])

x_initial = np.array([32.4658405015485, 34.4658405015485, 23.6481193219279, 47.2962386438559, 1])
solution = optimize.minimize(vector_lagrangian, x_initial)

fun: -2.140132414183526e+37
 hess_inv: array([[1, 0, 0, 0, 0],
       [0, 1, 0, 0, 0],
       [0, 0, 1, 0, 0],
       [0, 0, 0, 1, 0],
       [0, 0, 0, 0, 1]])
      jac: array([0., 0., 0., 0., 0.])
  message: 'Optimization terminated successfully.'
     nfev: 119
      nit: 1
     njev: 17
   status: 0
  success: True
        x: array([ 6.90178344e+08,  6.51257507e+08,  9.75839219e+08,  4.87919645e+08,
       -5.08835272e+06])
'constraint': '680,080,111,963'

在这种情况下,约束没有被保持,它收敛于局部最小值。为什么会这样?问题是由求解器引起的,是正在优化的具体功能,还是有其他原因?

【问题讨论】:

  • 欢迎来到Stack Overflow.!要求有关问题方法的一般指导的问题通常过于宽泛,不适合本网站。人们有自己解决问题的方法,因此不可能有正确的答案。仔细阅读 Where to StartMinimal Reproducible Example,然后编辑您的帖子。
  • 我想说这个问题非常具体。我正在通过数值方法运行优化并且得到不正确的结果。我想知道为什么这个特定情况下没有约束。同时,可以通过运行上面的函数来重现结果,这里没有额外的链接。最后,欢迎提出进一步的问题。如果需要,我不介意详细说明。
  • 您所说的问题是“有谁知道这个求解器没有做什么?”。没有具体的“正确”答案可以解决这个问题。因此,这个问题只能通过意见来回答。
  • @itprorh66 我同意这个问题应该更加集中和减少。但是,我不同意这个问题只能通过意见来回答。求解器没有产生正确解的原因很简单:OP 的数学是错误的。
  • 感谢您的反馈。我将编辑问题。

标签: python optimization scipy constraints numerical-methods


【解决方案1】:

正如 cmets 中已经提到的,您的数学是错误的,因为最小化拉格朗日不会产生相应优化问题的局部最小值。假设 f : R^n -> R 和 g : R^n -> R^m 都是可微函数,你想解决优化问题

min f(x) s.t. g(x) = 0

那么一阶必要最优性条件(FOC)为

∇L(x, λ) = ∇f(x) + ∇g(x)^T * λ = 0
                          g(x) = 0

其中 L 是拉格朗日,∇f 是目标梯度,∇g 是函数 g 的转置雅可比。因此,您需要找到函数 H(x,λ) = (∇f(x) + λ^T * ∇g(x), g(x))^T 的根来求解 FOC,它可以是通过scipy.optimize.root完成:

from scipy.optimize import minimize, root
from scipy.optimize._numdiff import approx_derivative

def damage_min(x):
    return -100*(1+x[0])*(1+x[1])*(1+x[2]*x[3])

def constraints(x):
    c1 = x[0] - x[1]
    c2 = 2*x[2] - x[3]
    c3 = x[0]/x[3] - 0.466/0.622
    c4 = x[3] - 0.466
    return np.array([c1, c2, c3, c4])

def f_grad(x):
    return approx_derivative(damage_min, x)

def g_jac(x):
    return approx_derivative(constraints, x)

def H(z, f_grad, g, g_jac):
    g_evaluated = g(z)
    x, λ   = np.split(z, (-g_evaluated.size, ))
    eq1 = f_grad(x) + g_jac(x).T @ λ
    eq2 = g_evaluated
    return np.array([*eq1, *eq2])

# res.x contains the solution
res = root(lambda z: H(z, f_grad, constraints, g_jac), x0=np.ones(8))

产生解(由 x 和拉格朗日乘数 λ 组成):

array([ 3.49125402e-01,  3.49125402e-01,  2.33000000e-01,  4.66000000e-01,
       -1.49561074e+02,  4.24092469e+01,  1.39390921e+02,  3.08919653e+02])

几点说明:

  • 一般来说,强烈建议提供精确的导数,而不是通过 approx_derivate 的有限差分来逼近它们。
  • 如果您真的想要解决最小化问题,可以通过最小化函数 H 的欧几里得范数来解决 FOC。这正是 root 方法在幕后所做的。

【讨论】:

  • 非常感谢您的清晰解释。使用函数 H 是结合拉格朗日条件的绝佳方式。 approx_derivate 也很强大。就我而言,我计划从目标函数的自然对数中提供梯度的精确导数。据我了解,单调变换不会产生梯度的新问题。同时,鉴于大多数约束在本质上都是线性的,希望 approx_derivate 在后台使用的有限差分近似就足够了。