优化类问题

优化类问题

跟据目标函数与约束问题分类

线性规划

线性规划问题一般是最简单的，经常直接看，比如用 单纯型法 等即可求解，在此不多做赘述。举个例子：

\[min \qquad z = 2x_1 + 3x_2 + x_3 \\ s.t. \begin{cases} x_1 + 4x_2 + 2x_3\ge 8 \\ 3x_1 + 2x_2 \ge 6 \\ x_1 + x_2 + x_3 = 7 \end{cases} \qquad x_1,x_2,x_3 \ge 0 \]

在利用程序 linprog 求解前需要做出一定变形：

• 必须求的是最小值
• 不等式约束必须改成小于等于

import numpy as np
from scipy import optimize as op

def solution():
    # 定义决策变量范围
    x1 = (0, None)
    x2 = (0, None)
    x3 = (0, None)

    # 定义目标函数系数
    c = np.array([2, 3, 1])

    # 定义约束条件系数
    A_ub = np.array([[-1, -4, -2], [-3, -2, 0]])    # 不等式系数
    B_ub = np.array([-8, -6])                       # 不等式常数项
    A_eq = np.array([1, 1, 1])  # 等式的系数
    B_eq = np.array([7])        # 等式的常数项

    res = op.linprog(c=c, a_ub=A_ub, b_ub=B_ub, a_eq=A_eq, b_eq=B_eq, bounds=(x1, x2, x3))
    return res

非线性规划

非线性问题没有一个定式解，只能通过自己写方法来实现求解。例如：

\[min \qquad z = x_1^2 + x_2^2 + x_3^2 + 8 \\ s.t. \begin{cases} x_1^2 - x_2 + x_3^2 \ge 0 \\ -x_1 + x_2^2 + x_3^2 \le 0 \\ -x_1 - x_2^2 + 2 = 0 \\ x_2 + 2x_3^2 \ge 0 \end{cases} \qquad x_1, x_2, x_3 \ge 0 \]

必须自己实现调用方法 minimize 来进行求解：

import numpy as np
from scipy.optimize import  minimize

# 定义目标函数 min z = x1^2 + x2^2 + x^3 + 8
def objective(x):
    return x[0] ** 2 + x[1] ** 2 + x[2] ** 2 + 8

# 定义约束条件
def constraint1(x):
    return x[0] ** 2 - x[1] + x[2] ** 2
def constraint2(x):
    return -(x[0] + x[1] ** 2 + x[2] ** 2 - 20)
def constraint3(x):
    return -x[0] - x[1] ** 2 + 2
def constraint4(x):
    return x[1] + 2 * x[2] ** 2 - 3

# 定义约束条件 ineq -> 大于等于  eq -> 等于
con1 = {'type': 'ineq', 'fun': constraint1}
con2 = {'type': 'ineq', 'fun': constraint2}
con3 = {'type': 'eq', 'fun': constraint3}
con4 = {'type': 'eq', 'fun': constraint4}
cons = ([con1, con2, con3, con4])

# 决策变量的符号约束
bounds = ((0.0, None), (0.0, None), (0.0, None))

def solution():
    x0 = np.array([0, 0, 0])  # 初始值
    ans = minimize(objective, x0, method='SLSQP', bounds=bounds, constraints=cons)
    return ans

其他无约束问题求极值方法

另外还有一些题目要求在某个函数上求其最小值。常用的方法有： 斐波那契探测法，牛顿迭代法，模拟退火 和 插值法 等。

斐波那契探测法

def fibonacciSearch(begin, end, eps=1e-5):
    rho = 1 - (sqrt(5) - 1) / 2
    x1 = begin + rho * (end - begin)
    x2 = end - rho * (end - begin)
    while x2 - x1 > eps:
        if f(x1) > f(x2):
            begin = x1
            x1 = x2
            x2 = end - rho * (end - begin)
        elif f(x1) < f(x2):
            end = x2
            x2 = x1
            x1 = begin + rho * (end - begin)
        else:
            begin = x1
            end = x2
            x1 = begin + rho * (end - begin)
            x2 = end - rho * (end - begin)
    return round((x1 + x2) / 2, 4)

注意事项：初始给定的区间内必须有且只有单峰极小值。在这里偷懒使用近似的 0.618 代替了斐波那契数的比值

牛顿迭代法

def Newton(begin, eps=1e-6):
    x = torch.tensor([begin * 1.0], requires_grad=True)
    while abs(f(x)) > eps:
        f_val = f(x)
        fp = torch.autograd.grad(f(x), x, create_graph=True)[0]
        x = x - f_val / fp
    return x.data[0]

由于求导机制的复杂性，某些问题不一定可以求出显式的导数式，因此在此采用 torch 库下的自动求导机制，转化为张量后求解。

注意事项：初始值必须给定比较靠近零点，否则可能会失真

除了最简单的牛顿迭代以外，由于多元函数求梯度时的 \(Hesse\) 矩阵求解时间复杂度较高，因此有改进版本的 拟牛顿法 以及 BFGS 等在此不做赘述

模拟退火

class SA:
    def __init__(self, begin, end):
        self.begin = begin  # 边界左端点
        self.end = end      # 边界右端点
        self.iter = 100     # 内循环迭代次数，即 L = 100
        self.alpha = 0.99   # 降温系数， alpha = 0.99
        self.T0 = 200       # 初始温度为 T0 = 200
        self.Tf = 0.01      # 结束温度为 Tf = 0.01
        self.T = 100        # 当前温度
        self.x = np.linspace(self.begin, self.end, self.iter)

    def generate_new(self, x):       # 扰动产生新解
        while True:
            x_new = x + self.T * (random() - random())
            if self.begin <= x_new <= self.end:
                break
        return x_new

    def Metropolis(self, f, f_new):     # Metropolis 准则
        if f_new <= f:
            return True
        return random() < math.exp((f - f_new) / self.T)

    def best(self):     # 计算当前状态下最优解
        f_min = f(self.x[0])
        pos = self.x[0]
        for xp in self.x:
            if f(xp) < f_min:
                f_min = f(xp)
                pos = xp
        return f_min, pos

    def run(self):
        count = 1
        while self.T > self.Tf:
            for i in range(self.iter):
                f_prime = f(self.x[i])
                x_new = self.generate_new(self.x[i])
                f_new = f(x_new)
                if self.Metropolis(f_prime, f_new):		# 判断是否接受新解
                    self.x[i] = x_new
            self.T = self.T * self.alpha
            # f_min, pos = self.best()
            # print('condition {0}: f_val = {1}, x_pos = {2}'.format(count, f_min, pos))
            count += 1
        f_min, pos = self.best()
        return round(f_min, 4), round(pos, 4)

相比于上述算法，模拟退火实现的是对全局最优解的查找。上面的算法寻找的是初始值附近的局部最优解，因此使用之前还需要判断是否为凸函数或者确定定义域等。模拟退火没有这些事，直接使用即可。不过由于其不稳定性最好多算几次求个平均。

辛普森自适应积分

连续性问题不可避免涉及到积分运算，补充一个数值积分内常用的积分

def simpson(left, right, eps=1e-4):
    def cal(left, right):
        mid = (left + right) / 2.0
        return (f(left) + 4.0 * f(mid) + f(right)) * (right - left) / 6.0
    mid = (left + right) / 2.0
    ST = cal(left, right)
    SL = cal(left, mid)
    SR = cal(mid, right)
    if math.fabs(SL + SR - ST) < 15.0 * eps:
        return SL + SR + (SL + SR - ST) / 15.0
    return simpson(left, mid, eps=eps/2.0) + simpson(mid, right, eps=eps/2.0)

二次规划

二次规划解决的是二次型下的最优解问题，其形式如：

\[\underset{x}{min} \qquad \frac{1}{2}x^TPx + q^Tx \\ subject \quad to \begin{cases} Gx \preceq h \\ Ax = b \end{cases} \]

其中我们需要对矩阵 \(P\) 进行讨论：

• 当 \(P\) 是正定矩阵时存在全局最优解。当 \(P\) 是半正定矩阵时目标函数为凸函数，存在最优解(不一定唯一)
• 当 \(P\) 是不定矩阵时目标函数非凸，存在多个局部最小值和稳定点

以题目为例：

\[min \qquad 2x_1^2 + x_2^2 + x_1x_2 + x_1 + x_2 \\ s.t. \begin{cases} -x_1 \le 0 \\ -x_2 \le 0 \\ x_1 + x_2 = 1 \end{cases} \]

有求解代码：

from cvxopt import matrix, solvers

# cvxopt.matrix 为列优先
P = matrix([[4.0, 1.0], [1.0, 2.0]])        # 二次项对应系数矩阵
q = matrix([1.0, 1.0])                      # 一次项系数
G = matrix([[-1.0, 0.0], [0.0, -1.0]])      # 不等式约束系数
h = matrix([0.0, 0.0])                      # 不等式约束常数
A = matrix([1.0, 1.0], (1, 2))              # 等式约束系数
b = matrix([1.0])                           # 等式约束常数

def solution():
    result = solvers.qp(P, q, G, h, A, b)
    return result['x']

按照控制变量类型分类

整数规划

常用的方法有： 分枝定界法，割平面法 和 隐枚举法 等，匈牙利算法放在二分图下的指派问题，在此不做解释。

另外必须明确的一点是，目前还没有完整的解决整数规划的算法，在此的讨论全部默认是 整数线性规划。

给出分支定界法代码：

'''
# param c: 目标函数系数
# param A_ub: 不等式系数
# param B_ub: 不等式常数
# param A_eq: 等式系数
# param B_eq: 等式常数
'''
def branch_and_bound(c, A_ub=None, B_ub=None, A_eq=None, B_eq=None, eps=1e-5):
    res = linprog(c, A_ub, B_ub, A_eq, B_eq)
    bestVal = sys.maxsize
    bestX = res.x
    if not(type(res.x) is float or res.status != 0):
        bestVal = sum([x*y for x, y in zip(c, bestX)])
    if all((x-math.floor(x)) < eps or (math.ceil(x)-x) < eps for x in bestX):
        return round(bestVal, 5), [round(x, 5) for x in bestX]
    else:
        ind = [i for i, x in enumerate(bestX) if (x-math.floor(x)) > eps and (math.ceil(x)-x) > eps][0]
        newCon1 = [0] * len(A_ub[0])
        newCon2 = [0] * len(A_ub[0])
        newCon1[ind] = -1
        newCon2[ind] = 1
        newA1 = A_ub.copy()
        newA2 = A_ub.copy()
        newA1.append(newCon1)
        newA2.append(newCon2)
        newB1 = B_ub.copy()
        newB1.append(-math.ceil(bestX[ind]))
        newB2 = B_ub.copy()
        newB2.append(math.floor(bestX[ind]))
        res1 = branch_and_bound(c, newA1, newB1, A_eq, B_eq, eps)
        res2 = branch_and_bound(c, newA2, newB2, A_eq, B_eq, eps)
        if res1[0] < res2[0]:
            return res1
        return res2

割平面法整体思路与分支定界法类似，不作单独的代码整理

混合整数规划

对于这种可以采用 MILP 算法进行求解，使用 \(\text{MATLAB}\) 中的函数 intlingrog 即可。

01整数规划

对于这类二分图带权最大匹配问题，使用的是匈牙利算法(或者在 ACM 里实际用的是 KM 算法？)

匈牙利算法

这种算法对应是 \(\text{KM}\) 算法的弱化版本，因为美赛中大部分指派问题都是完全二分图 (二分图两侧任意两点之间均有边)，因此简化出这种简便算法。给出代码：

from scipy.optimize import linear_sum_assignment
def Hungarian(cost):
    # cost[i][j] 表示第 i 个工人做第 j 项任务的工作成本，利用匈牙利算法以实现工作成本的最小化
    row_ind, col_ind = linear_sum_assignment(cost)
    tot_cost = cost[row_ind, col_ind].sum()
    for i in range(len(row_ind)):
        print('第 {0} 个人做第 {1} 项工作，其成本为 {2}'.format(row_ind[i]+1, col_ind[i]+1, cost[row_ind[i]][col_ind[i]]))
    print('总成本为：{0}'.format(tot_cost))

KM 算法

人不可能全能，故有些人是无法做某一部分工作的，这样导致了其不是一个完全二分图，就应用到了这种算法。具体算法思想是贪心，不断找当前权值最大的策略，如果找不到就下降一个单位重新寻找直至匹配结束。附上代码链接：KM算法代码实现

智能算法

除去上述的模拟退火以外，还需要补充一个遗传算法。以求解该函数的最大值为例：

\[F(x,y) = 3(1-x)^2e^{-x^2-(y+1)^2}-10(\frac{x}{5} - x^3 - y^5)e^{-x^2-y^2}-\frac{1}{3}e^{-(x+1)^2 - y^2} \\ x \in [-3, 3], y \in [-3, 3] \]

给出代码有:

def F(x, y):
    return  3*(1-x)**2*np.exp(-(x**2)-(y+1)**2) - \
            10*(x/5 - x**3 - y**5)*np.exp(-x**2-y**2) - \
            1/3**np.exp(-(x+1)**2 - y**2)

class GA:
    def __init__(self, X_BOUND, Y_BOUND):
        self.DNA_SIZE = 24              # 编码长度
        self.POP_SIZE = 500             # 种群数量
        self.CROSSOVER_RATE = 0.8       # 交叉概率
        self.MUTATION_RATE = 0.05       # 变异概率
        self.N_GENERATIONS = 500        # 遗传代数
        self.X_BOUND = X_BOUND         # x 取值范围
        self.Y_BOUND = Y_BOUND          # y 取值范围
        self.pop = np.random.randint(2, size=(self.POP_SIZE, self.DNA_SIZE*2))    # 种群矩阵

    def translationDNA(self):           # 解码为具体的坐标值
        x_pop = self.pop[:, 1::2]
        y_pop = self.pop[:, ::2]
        x = x_pop.dot(2 ** np.arange(self.DNA_SIZE)[::-1]) / float(2 ** self.DNA_SIZE - 1) \
            * (self.X_BOUND[1] - self.X_BOUND[0]) + self.X_BOUND[0]
        y = y_pop.dot(2 ** np.arange(self.DNA_SIZE)[::-1]) / float(2 ** self.DNA_SIZE - 1) \
            * (self.Y_BOUND[1] - self.Y_BOUND[0]) + self.Y_BOUND[0]
        return x, y

    def get_fitness(self):          # 计算适应度
        x, y = self.translationDNA()
        pred = F(x, y)
        return pred - np.min(pred)

    def select(self, fitness):      # 根据不同适应度选择
        idx = np.random.choice(np.arange(self.POP_SIZE), size=self.POP_SIZE, replace=True,
                               p=fitness/(fitness.sum()))
        return self.pop[idx]

    def mutation(self, child):      # 变异
        if np.random.rand() < self.MUTATION_RATE:
            mutate_point = np.random.randint(0, self.DNA_SIZE * 2)
            child[mutate_point] = child[mutate_point] ^ 1
        return child

    def crossover_and_mutation(self):       # 遗传交叉与变异
        new_pop = []
        for father in self.pop:
            child = father
            if np.random.rand() < self.CROSSOVER_RATE:
                mother = self.pop[np.random.randint(self.POP_SIZE)]
                cross_point = np.random.randint(self.DNA_SIZE*2)
                child[cross_point:] = mother[cross_point:]
            child = self.mutation(child)
            new_pop.append(child)
        return new_pop

    def run(self):
        for _ in range(self.N_GENERATIONS):
            self.pop = np.array(self.crossover_and_mutation())
            fitness = self.get_fitness()
            self.pop = self.select(fitness)
        self.printInfo()

    def printInfo(self):
        fitness = self.get_fitness()
        max_fitness_idx = np.argmax(fitness)
        print('max_fitness: ', fitness[max_fitness_idx])
        x, y = self.translationDNA()
        print('最优的基因型为: ', self.pop[max_fitness_idx])
        print('(x, y) = ', (x[max_fitness_idx], y[max_fitness_idx]))
        print('F(x, y) = ', F(x[max_fitness_idx], y[max_fitness_idx]))

其中由于遗传算法需要映射到单独的染色体上，因此解决的问题往往是离散型的。而模拟退火可以随机"袋鼠跳"到任意一个可能的取值上，常常用来解决连续型全局最优解的问题。

其他分类方法

多目标规划模型

有时候所需要求解的模型中，往往我们想要让多个目标取到最值，尽可能优化，这时候不同解下的不同目标空间则张成了一个解空间。在这个解空间内，我们无法直接去比较两个解的好坏(除非存在偏序关系，即多个目标之间均存在一致的某个解优于另一个解)，由此产生出了劣解，非劣解，绝对最优解，有效解和弱有效解等概念。为了解决这个问题，我们就需要认为去确定一个比较函数。

最直接的想法可以采用欧几里得距离，通过比较这个点与原点之间的距离以确定解的好坏。或者我们也可以对不同的目标之间给予加权，针对现实问题的轻重具体问题具体分析。但是总的思路均是转化为单目标规划去求解。

• 主要目标函数法：选择其中一个目标最优化，其他目标放宽条件，只需满足一定条件即可
• 线性加权和目标函数法：按照一定的规则分别给不同的目标赋予权重，并线性组合为一个新的评价目标，以该目标为最优化目标
• 分层序列法：把规划中的各个目标按照重要程度排序，以此求单目标规划下的最优解，最后得到的解序列依次进行闭区间套操作。

动态规划

哈哈哈这个老生常谈的话题。作为 ACM 老手和菜鸟的分水岭，做这类 dp 问题的熟练程度可以很明显分出一个人算法水平的好坏。其基本思想是：将一个问题划分为多个子问题的求解，如果子问题达到最优解的话，那么这个问题也会达到最优解，称之为 最优子结构 的性质。具体的入门题目必须掌握 01背包问题，当然进阶的优化的滚动数组之类的先不谈。之后必须了解的是动态规划与贪心的区别，两者在数学证明上有根本性的区分，如何选择两种解法十分关键。具体可以比较一下完全背包和01背包解法的异同，即可简单了解。另外分类还有如 线性dp，区间dp，数位dp，插头dp 等，优化方式有例如四边形优化，斜率优化，状压(转移方程为一个矩阵)等问题，可以进一步学习研究。

随机规划

相比于上述的线性规划或者非线性规划，其每一项前面的系数均给出确定值，而在随机规划中，会对这个确定的系数动手脚。作为《运筹学》中较为高级的优化方法，随机规划关注的是不确定情形下的目标函数的期望，而鲁棒优化致力于将最坏情况变得更好(相对更保守)。在随机规划中，会给出各个系数的一个随机变量，给出其均值或者方差，或者遵循的分布等，然后去求目标函数的最优解。针对这种问题，我们可以根据不同的分布去采样，转化为确定的规划之后求解，然后多次求解后取平均即可得到目标函数的期望值。这种随机规划往往可以放在对一个最优化模型的灵敏性分析中。