Python数模笔记-模拟退火算法（3）整数规划问题

1、整数规划问题

　　整数规划问题在工业、经济、国防、医疗等各行各业应用十分广泛，是指规划中的变量（全部或部分）限制为整数，属于离散优化问题（Discrete Optimization）。
　　线性规划问题的最优解可能是分数或小数。但很多实际问题常常要求某些变量必须是整数解，例如：机器的台数、工作的人数或装货的车数。根据对决策变量的不同要求，整数规划又可以分为：纯整数规划、混合整数规划、0-1整数规划、混合0－1规划。
　　整数规划与线性规划的差别只在于增加了整数约束。初看起来似乎只要把线性规划得到的非整数解舍入化整就可以得到整数解，但是这样化整后的整数解不一定是最优解，甚至可能不是可行解。因此，通常需要采用特殊的方法来求解整数规划，这比求解线性规划问题复杂的多，以至于至今还没有一般的多项式解法。因此，整数规划问题被看作数学规划中、甚至是数学中最困难的问题之一。
　　求解整数规划比较成功又流行的方法是分支定界法和割平面法。核心思想是把整数规划问题分解为一系列线性规划问题，并追踪整数规划问题的上界（最优可行解）和下界（最优线性松弛解），逐步迭代收敛到最优解。由于精确算法为指数复杂度，因此在有限时间内也不能获得全局最优解，只能获得近似最优解。
　　目前整数规划问题的优化求解器主要有：IBM Cplex，Gurobi，FICO Xpress，SCIP，2018年中科院发布了CMIP混合整数规划求解器。使用 Lingo 可以求解整数规划问题，使用 Matlab 也可以用intlinprog 函数求解整数规划问题，实际上都是使用软件中内建的求解器。Python 也可以使用第三方库求解整数规划问题，例如 Cvxpy、PuLp 都可以求解整数规划问题，Cplex、Gurobi也有自己的python API。　　

欢迎关注 Youcans 原创系列，每周更新数模笔记

Python数模笔记-PuLP库
 Python数模笔记-StatsModels统计回归
 Python数模笔记-Sklearn
Python数模笔记-NetworkX
Python数模笔记-模拟退火算法

2、模拟退火算法处理整数约束

　　由于整数规划问题在有限时间内不能获得全局最优解，启发式算法就有了用武之地。下面我们讨论模拟退火算法处理整数约束，求解整数规划问题。
　　上一篇文章中我们讨论模拟退火算法处理线性规划的约束条件时，方法比其它常用算法复杂的多。但是，模拟退火算法在处理整数约束时，方法却极其简单：
　　对于决策变量为连续变量的一般优化问题，基本的模拟退火算法在决策变量的取值范围随机产生初始解，新解则是在现有解的邻域施加扰动产生，算法上通过均匀分布或正态分布的随机数来实现：

xInitial = random.uniform(xMin, xMax)
# random.uniform(min,max) 在 [min,max] 范围内随机生成一个实数

xNew = xNow + scale * (xMax-xMin) * random.normalvariate(0, 1)
# random.normalvariate(0, 1)：产生服从均值为0、标准差为 1 的正态分布随机实数
xNew = max(min(xNew, xMax), xMin) # 保证新解在 [min,max] 范围内

　　对于整数规划问题，只要将产生初值/新解的随机实数发生器 random.uniform、random.normalvariate 改为随机整数发生器 random.randint即可：

xInitial = random.randint(xMin, xMax)
# random.randint(xMin, xMax) 产生 [min,max]之间的随机整数

　　由于模拟退火算法与问题无关（Problem-independent），所以通常来说这样处理并不会影响算法的性能：既不会引起不可行解，也不用担心得不到最优解——近似算法只能得到近似最优解的，而且可以得到近似最优解。
　　既然如此，更简单的处理方法，连随机整数发生器都不需要，直接把线性规划得到的非整数解舍入化整就可以了：
　　

xNew = round(xNow + scale * (xMax-xMin) * random.normalvariate(0, 1))
# random.normalvariate(0, 1)：产生服从均值为0、标准差为 1 的正态分布随机实数
xNew = max(min(xNew, xMax), xMin) # 保证新解在 [min,max] 范围内

　　这样处理的好处是：（1）简单、直接，（2）便于实现所需的概率分布。

3、数模案例

　　为了便于理解，本文仍使用之前的案例。

3.1 问题描述：

　　某厂生产甲乙两种饮料，每百箱甲饮料需用原料6千克、工人10名，获利10万元；每百箱乙饮料需用原料5千克、工人20名，获利9万元。
　　今工厂共有原料60千克、工人150名，又由于其他条件所限甲饮料产量不超过8百箱。
　　（5）若不允许散箱（按整百箱生产），如何安排生产计划，即两种饮料各生产多少使获利最大？

3.2 问题分析：

　　问题（5）要求按整百箱生产，即要求决策变量为整数，是整数规划问题。
　　对于模拟退火算法，基本算法中的初值/新解都是随机生成的浮点实数（均匀分布或正态分布）。对于整数规划问题，只要将产生初值/新解的随机实数发生器改为随机整数发生器即可，或者把线性规划得到的非整数解舍入化整。
　　

3.3 问题建模：

　　决策变量：
　　　　x1：甲饮料产量，正整数（单位：百箱）
　　　　x2：乙饮料产量，正整数（单位：百箱）
　　目标函数：
　　　　max fx = 10x1 + 9x2
　　约束条件：
　　　　6x1 + 5x2 <= 60
　　　　10x1 + 20x2 <= 150
　　取值范围：
　　　　给定条件：x1, x2 >= 0，x1 <= 8
　　　　推导条件：由 x1,x2>=0 和 10*x1+20*x2<=150 可知：0<=x1<=15；0<=x2<=7.5
　　　　因此，0 <= x1<=8，0 <= x2<=7.5

3.4 惩罚函数法求解约束优化问题：

　　构造惩罚函数：
　　　　p1 = (max(0, 6*x1+5*x2-60))**2
　　　　p2 = (max(0, 10*x1+20*x2-150))**2
　　说明：如存在等式约束，例如：x1 + 2*x2 = m，也可以转化为惩罚函数：
　　　　p3 = (x1+2*x2-m)**2
　　　　P(x) = p1 + p2 + ...
　　构造增广目标函数：
　　　　L(x,m(k)) = min(fx) + m(k)*P(x)
　　　　m(k)：惩罚因子，随迭代次数 k 逐渐增大

　　在模拟退火算法中，m(k) 随外循环迭代次数逐渐增大，但在内循环中应保持不变。

4、模拟退火算法 Python 程序：求解整数规划问题

# 模拟退火算法 程序：求解线性规划问题（整数规划）
# Program: SimulatedAnnealing_v4.py
# Purpose: Simulated annealing algorithm for function optimization
# v4.0: 整数规划：满足决策变量的取值为整数（初值和新解都是随机生成的整数）
# Copyright 2021 YouCans, XUPT
# Crated：2021-05-01

#  -*- coding: utf-8 -*-
import math                         # 导入模块
import random                       # 导入模块
import pandas as pd                 # 导入模块 YouCans, XUPT
import numpy as np                  # 导入模块 numpy，并简写成 np
import matplotlib.pyplot as plt     
from datetime import datetime

#=== 关注 Youcans，分享更多原创系列 https://www.cnblogs.com/youcans/ ===
# 子程序：定义优化问题的目标函数
def cal_Energy(X, nVar, mk): 	# m(k)：惩罚因子，随迭代次数 k 逐渐增大
    p1 = (max(0, 6*X[0]+5*X[1]-60))**2
    p2 = (max(0, 10*X[0]+20*X[1]-150))**2
    fx = -(10*X[0]+9*X[1])
    return fx+mk*(p1+p2)


# 子程序：模拟退火算法的参数设置
def ParameterSetting():
    cName = "funcOpt"           # 定义问题名称 YouCans, XUPT
    nVar = 2                    # 给定自变量数量，y=f(x1,..xn)
    xMin = [0, 0]               # 给定搜索空间的下限，x1_min,..xn_min
    xMax = [8, 8]               # 给定搜索空间的上限，x1_max,..xn_max

    tInitial = 100.0            # 设定初始退火温度(initial temperature)
    tFinal  = 1                 # 设定终止退火温度(stop temperature)
    alfa    = 0.98              # 设定降温参数，T(k)=alfa*T(k-1)
    meanMarkov = 100            # Markov链长度，也即内循环运行次数
    scale   = 0.5               # 定义搜索步长，可以设为固定值或逐渐缩小
    return cName, nVar, xMin, xMax, tInitial, tFinal, alfa, meanMarkov, scale


# 模拟退火算法
def OptimizationSSA(nVar,xMin,xMax,tInitial,tFinal,alfa,meanMarkov,scale):
    # ====== 初始化随机数发生器 ======
    randseed = random.randint(1, 100)
    random.seed(randseed)  # 随机数发生器设置种子，也可以设为指定整数

    # ====== 随机产生优化问题的初始解 ======
    xInitial = np.zeros((nVar))   # 初始化，创建数组
    for v in range(nVar):
        # xInitial[v] = random.uniform(xMin[v], xMax[v]) # 产生 [xMin, xMax] 范围的随机实数
        xInitial[v] = random.randint(xMin[v], xMax[v]) # 产生 [xMin, xMax] 范围的随机整数
    # 调用子函数 cal_Energy 计算当前解的目标函数值
    fxInitial = cal_Energy(xInitial, nVar, 1) # m(k)：惩罚因子，初值为 1

    # ====== 模拟退火算法初始化 ======
    xNew = np.zeros((nVar))         # 初始化，创建数组
    xNow = np.zeros((nVar))         # 初始化，创建数组
    xBest = np.zeros((nVar))        # 初始化，创建数组
    xNow[:]  = xInitial[:]          # 初始化当前解，将初始解置为当前解
    xBest[:] = xInitial[:]          # 初始化最优解，将当前解置为最优解
    fxNow  = fxInitial              # 将初始解的目标函数置为当前值
    fxBest = fxInitial              # 将当前解的目标函数置为最优值
    print('x_Initial:{:.6f},{:.6f},\tf(x_Initial):{:.6f}'.format(xInitial[0], xInitial[1], fxInitial))

    recordIter = []                 # 初始化，外循环次数
    recordFxNow = []                # 初始化，当前解的目标函数值
    recordFxBest = []               # 初始化，最佳解的目标函数值
    recordPBad = []                 # 初始化，劣质解的接受概率
    kIter = 0                       # 外循环迭代次数，温度状态数
    totalMar = 0                    # 总计 Markov 链长度
    totalImprove = 0                # fxBest 改善次数
    nMarkov = meanMarkov            # 固定长度 Markov链

    # ====== 开始模拟退火优化 ======
    # 外循环，直到当前温度达到终止温度时结束
    tNow = tInitial                 # 初始化当前温度(current temperature)
    while tNow >= tFinal:           # 外循环，直到当前温度达到终止温度时结束
        # 在当前温度下，进行充分次数(nMarkov)的状态转移以达到热平衡
        kBetter = 0                 # 获得优质解的次数
        kBadAccept = 0              # 接受劣质解的次数
        kBadRefuse = 0              # 拒绝劣质解的次数

        # ---内循环，循环次数为Markov链长度
        for k in range(nMarkov):    # 内循环，循环次数为Markov链长度
            totalMar += 1           # 总 Markov链长度计数器

            # ---产生新解
            # 产生新解：通过在当前解附近随机扰动而产生新解，新解必须在 [min,max] 范围内
            # 方案 1：只对 n元变量中的一个进行扰动，其它 n-1个变量保持不变
            xNew[:] = xNow[:]
            v = random.randint(0, nVar-1)   # 产生 [0,nVar-1]之间的随机数
            xNew[v] = round(xNow[v] + scale * (xMax[v]-xMin[v]) * random.normalvariate(0, 1))
            # 满足决策变量为整数，采用最简单的方案：产生的新解按照四舍五入取整
            xNew[v] = max(min(xNew[v], xMax[v]), xMin[v])  # 保证新解在 [min,max] 范围内

            # ---计算目标函数和能量差
            # 调用子函数 cal_Energy 计算新解的目标函数值
            fxNew = cal_Energy(xNew, nVar, kIter)
            deltaE = fxNew - fxNow

            # ---按 Metropolis 准则接受新解
            # 接受判别：按照 Metropolis 准则决定是否接受新解
            if fxNew < fxNow:  # 更优解：如果新解的目标函数好于当前解，则接受新解
                accept = True
                kBetter += 1
            else:  # 容忍解：如果新解的目标函数比当前解差，则以一定概率接受新解
                pAccept = math.exp(-deltaE / tNow)  # 计算容忍解的状态迁移概率
                if pAccept > random.random():
                    accept = True  # 接受劣质解
                    kBadAccept += 1
                else:
                    accept = False  # 拒绝劣质解
                    kBadRefuse += 1

            # 保存新解
            if accept == True:  # 如果接受新解，则将新解保存为当前解
                xNow[:] = xNew[:]
                fxNow = fxNew
                if fxNew < fxBest:  # 如果新解的目标函数好于最优解，则将新解保存为最优解
                    fxBest = fxNew
                    xBest[:] = xNew[:]
                    totalImprove += 1
                    scale = scale*0.99  # 可变搜索步长，逐步减小搜索范围，提高搜索精度
                    
        # ---内循环结束后的数据整理
        # 完成当前温度的搜索，保存数据和输出
        pBadAccept = kBadAccept / (kBadAccept + kBadRefuse)  # 劣质解的接受概率
        recordIter.append(kIter)  # 当前外循环次数
        recordFxNow.append(round(fxNow, 4))  # 当前解的目标函数值
        recordFxBest.append(round(fxBest, 4))  # 最佳解的目标函数值
        recordPBad.append(round(pBadAccept, 4))  # 最佳解的目标函数值

        if kIter%10 == 0:                           # 模运算，商的余数
            print('i:{},t(i):{:.2f}, badAccept:{:.6f}, f(x)_best:{:.6f}'.\
                format(kIter, tNow, pBadAccept, fxBest))

        # 缓慢降温至新的温度，降温曲线：T(k)=alfa*T(k-1)
        tNow = tNow * alfa
        kIter = kIter + 1
        fxBest = cal_Energy(xBest, nVar, kIter)  # 由于迭代后惩罚因子增大，需随之重构增广目标函数
        # ====== 结束模拟退火过程 ======

    print('improve:{:d}'.format(totalImprove))
    return kIter,xBest,fxBest,fxNow,recordIter,recordFxNow,recordFxBest,recordPBad

# 结果校验与输出
def ResultOutput(cName,nVar,xBest,fxBest,kIter,recordFxNow,recordFxBest,recordPBad,recordIter):
    # ====== 优化结果校验与输出 ======
    fxCheck = cal_Energy(xBest, nVar, kIter)
    if abs(fxBest - fxCheck)>1e-3:   # 检验目标函数
        print("Error 2: Wrong total millage!")
        return
    else:
        print("\nOptimization by simulated annealing algorithm:")
        for i in range(nVar):
            print('\tx[{}] = {:.1f}'.format(i,xBest[i]))
        print('\n\tf(x) = {:.1f}'.format(cal_Energy(xBest,nVar,0)))

    return


# 主程序
def main():

    # 参数设置，优化问题参数定义，模拟退火算法参数设置
    [cName, nVar, xMin, xMax, tInitial, tFinal, alfa, meanMarkov, scale] = ParameterSetting()
    # print([nVar, xMin, xMax, tInitial, tFinal, alfa, meanMarkov, scale])

    # 模拟退火算法    [kIter,xBest,fxBest,fxNow,recordIter,recordFxNow,recordFxBest,recordPBad] \
        = OptimizationSSA(nVar,xMin,xMax,tInitial,tFinal,alfa,meanMarkov,scale)
    # print(kIter, fxNow, fxBest, pBadAccept)

    # 结果校验与输出
    ResultOutput(cName, nVar,xBest,fxBest,kIter,recordFxNow,recordFxBest,recordPBad,recordIter)


if __name__ == '__main__':
    main()

5、运行结果

Optimization by simulated annealing algorithm:
    x[0] = 8.0
    x[1] = 2.0
    f(x) = -98.0

参考文献：

（1）田澎，杨自厚，张嗣瀛，一类非线性整数规划的模拟退火求解，1993年控制理论及其应用年会论文集，海洋出版社，1993，533-537.