有约束的优化问题之 Lagrangian 乘子法、投影梯度算法和罚函数法¶
使用 Lagrangian 乘子法、投影梯度算法、罚函数法求解有约束的优化问题。
Lagrangian 乘子法¶
Lagrangian 乘子的具体形式¶
\[
L(x, y, \lambda) = x^2 + 5y^2 + xy + \lambda(x^2 + y^2 - 1)
\]
对偶问题¶
\[
\max_{\lambda} g(\lambda)\\\
\text{ s.t. } \lambda \geq 0
\]
\[
\text{where } g(\lambda) = \min_{x, y} L(x, y, \lambda)
\]
问题描述：¶
投影梯度算法求原问题的解¶
算法步骤¶
先按照普通的无约束优化问题进行梯度下降，如果迭代后得到的自变量取值\(y_{k+1}\)不满足约束条件，则在满足约束条件的可行集中找一个距离\(y_{k+1}\)最近的点，作为\(x_{k+1}\)。如果迭代后得到的自变量取值\(y_{k+1}\)已经满足约束条件，则\(x_{k+1}\)就等于\(y_{k+1}\)。
Python
# 导入包
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
# 支持中文
matplotlib.rcParams["font.sans-serif"] = ["SimSun"]
matplotlib.rcParams["axes.unicode_minus"] = False
# 渲染公式
plt.rcParams["text.usetex"] = True
定义原函数、约束函数和梯度函数等¶
Python
# 定义原函数
def f(X):
x = X[0]
y = X[1]
return x**2 + 5 * y**2 + x * y
# 定义约束函数
def h(X):
x = X[0]
y = X[1]
return np.array([x**2 + y**2 - 1])
# 定义梯度函数
def grad_f(X):
x = X[0]
y = X[1]
return np.array([2 * x + y, 10 * y + x])
# 定义模长函数
def norm(vector):
return np.sqrt(sum([i**2 for i in vector]))
初始化参数，迭代实现投影梯度算法¶
Python
# 回溯线性搜索实现梯度下降
alpha = 0.4
beta = 0.5
epsilon = 1e-6
k = 0
X = np.array([2, 2])
t = 1
k_list = [k]
X_all = X
# 当梯度的模仍然高于精度要求时，需要一直迭代
while norm(grad_f(X)) > epsilon:
# 回溯线性搜索，检查目前的步长是否满足要求，若步长太大，则需要缩小步长
while f(np.subtract(X, t * grad_f(X))) > f(X) - alpha * t * grad_f(X).T @ grad_f(X):
t = beta * t
# 检查 X 是否在约束范围内
if h(X) > 0:
# 对于不满足约束的 X，需要将其投影到约束范围内
# 本题的约束函数是一个单位圆，因此投影的方法是将 X 向量的模缩放到 1
X = X / norm(X)
# 记录此次的自变量值
k = k + 1
k_list.append(k)
X = X - t * grad_f(X)
X_all = np.vstack((X_all, X))
# print("k = %d, X = %s, f(X) = %f" % (k, X, f(X)))
# 退出循环，说明已经得到了满足精度要求的梯度。此时的自变量值就是极小值点
print("极小值点为：x* = {:.2f}, y* = {:.2f}".format(X[0], X[1]))
绘制迭代变化图¶
Python
# 绘制迭代次数与自变量取值
fig = plt.figure(figsize=(5, 4), dpi=200)
plt.plot(
X_all[:, 0], "r-o", markerfacecolor="none", label=r"$x$"
) # 空心圆用 markerfacecolor='none',
plt.plot(X_all[:, 1], "b-o", markerfacecolor="none", label=r"$y$")
plt.xlabel("迭代次数", usetex=False) # 中文不用 tex 渲染，否则会报错
plt.ylabel("自变量的取值", usetex=False)
plt.legend()
罚函数方法求原问题的解¶
算法步骤¶
把不等式的约束函数作为惩罚，在前期先允许自变量不满足约束函数，后期随着\(u\)变得非常大，则不满足约束条件带来的惩罚也会非常大。这会迫使加入惩罚函数之后的最小值是满足约束条件的（相当于逃避了惩罚）。
定义罚函数和罚函数的梯度函数¶
Python
# 定义罚函数
def g(X):
return max(0, h(X)) ** 2
# 定义罚函数的梯度
def grad_g(X):
if h(X) > 0:
return np.array([2 * X[0], 2 * X[1]])
else:
return np.array([0, 0])
# 定义加入罚函数之后的目标函数
def penalty(X, u):
return f(X) + u * g(X)
# 定义加入罚函数之后的目标函数的梯度
def grad_penalty(X, u):
return grad_f(X) + u * grad_g(X)
初始化参数，迭代实现罚函数法¶
Python
# 初始化参数
u = 1
k = 0
X = np.array([2, 2])
t = 1
k_list = [k]
X_all = X
# 当梯度的模仍然高于精度要求，或者约束函数不被满足时，需要一直迭代
# 通过无约束优化算法求解加入罚函数之后的目标函数的极小值点
while norm(grad_f(X)) > epsilon or h(X) > epsilon:
while norm(grad_penalty(X, u)) > epsilon:
# 回溯线性搜索，检查目前的步长是否满足要求，若步长太大，则需要缩小步长
while penalty(np.subtract(X, t * grad_penalty(X, u)), u) > penalty(
X, u
) - alpha * t * grad_penalty(X, u).T @ grad_penalty(X, u):
t = beta * t
# 更新自变量值
X = X - t * grad_penalty(X, u)
# 记录此次的自变量值，它是加入罚函数之后的目标函数的极小值点
k = k + 1
k_list.append(k)
X_all = np.vstack((X_all, X))
# 增大惩罚参数
u = 10 * u
# 退出循环，说明已经得到了满足精度要求的梯度，并且也满足约束条件。此时的自变量值就是极小值点
print("极小值点为：x* = {:.2f}, y* = {:.2f}".format(X[0], X[1]))
绘制迭代变化图¶
Python
# 绘制迭代次数与自变量取值
fig = plt.figure(figsize=(5, 4), dpi=200)
plt.plot(
X_all[:, 0], "r-o", markerfacecolor="none", label=r"$x$"
) # 空心圆用 markerfacecolor='none',
plt.plot(X_all[:, 1], "b-o", markerfacecolor="none", label=r"$y$")
plt.xlabel("迭代次数", usetex=False) # 中文不用 tex 渲染，否则会报错
plt.ylabel("自变量的取值", usetex=False)
plt.legend()
