拉格朗日乘数数值方法

Question

我正在尝试优化以下受约束约束的函数：

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] Constraint 300 Max damage 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 ] Constraint 300 Max damage 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'

在这种情况下，约束没有被保持并且收敛于局部最小值。为什么会这样？问题出在求解器、正在优化的具体函数还是其他原因？

Answer 1

正如评论中已经提到的，您的数学是错误的，因为最小化拉格朗日量不会产生相应优化问题的局部最小值。假设 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 是 objective 梯度，∇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方法在幕后。

拉格朗日乘数数值方法

Lagrange multipliers numerical approach

python

optimization

scipy

constraints

numerical-methods