跳转至

有约束的优化问题之 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) \]

迭代次数与自变量取值 - 投影梯度算法

问题描述:

problem-statement

投影梯度算法求原问题的解

算法步骤

先按照普通的无约束优化问题进行梯度下降,如果迭代后得到的自变量取值\(y_{k+1}\)不满足约束条件,则在满足约束条件的可行集中找一个距离\(y_{k+1}\)最近的点,作为\(x_{k+1}\)。如果迭代后得到的自变量取值\(y_{k+1}\)已经满足约束条件,则\(x_{k+1}\)就等于\(y_{k+1}\)

projected-gradient-method

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]))
Text Only
## 极小值点为:x* = 0.00, y* = -0.00
Python
print("极小值为:f(x*, y*) = {:.2f}".format(f(X)))
Text Only
## 极小值为:f(x*, y*) = 0.00

绘制迭代变化图

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()

unnamed-chunk-4-1

罚函数方法求原问题的解

算法步骤

把不等式的约束函数作为惩罚,在前期先允许自变量不满足约束函数,后期随着\(u\)变得非常大,则不满足约束条件带来的惩罚也会非常大。这会迫使加入惩罚函数之后的最小值是满足约束条件的(相当于逃避了惩罚)。

penalty-function-method-1

penalty-function-method-2

定义罚函数和罚函数的梯度函数

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]))
Text Only
## 极小值点为:x* = 0.00, y* = -0.00
Python
print("极小值为:f(x*, y*) = {:.2f}".format(f(X)))
Text Only
## 极小值为:f(x*, y*) = 0.00

绘制迭代变化图

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()

unnamed-chunk-7-3

评论