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


1、整数规划问题

整数规划问题在工业、经济、国防、医疗等各行各业应用十分广泛,是指规划中的变量(全部或部分)限制为整数,属于离散优化问题(Discrete Optimization)。
线性规划问题的最优解可能是分数或小数。但很多实际问题常常要求某些变量必须是整数解,例如:机器的台数、工作的人数或装货的车数。根据对决策变量的不同要求,整数规划又可以分为:纯整数规划、混合整数规划、0-1整数规划、混合0-1规划。
整数规划与线性规划的差别只在于增加了整数约束。初看起来似乎只要把线性规划得到的非整数解舍入化整就可以得到整数解,但是这样化整后的整数解不一定是最优解,甚至可能不是可行解。因此,通常需要采用特殊的方法来求解整数规划,这比求解线性规划问题复杂的多,以至于至今还没有一般的多项式解法。因此,整数规划问题被看作数学规划中、甚至是数学中最困难的问题之一。
求解整数规划比较成功又流行的方法是分支定界法和割平面法。核心思想是把整数规划问题分解为一系列线性规划问题,并追踪整数规划问题的上界(最优可行解)和下界(最优线性松弛解),逐步迭代收敛到最优解。由于精确算法为指数复杂度,因此在有限时间内也不能获得全局最优解,只能获得近似最优解。YouCans
目前整数规划问题的优化求解器主要有: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 = 10*x1 + 9*x2
  约束条件:
    6*x1 + 5*x2 <= 60
    10*x1 + 20*x2 <= 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
# = 关注 Youcans,分享原创系列 https://blog.csdn.net/youcans =
#  -*- 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# 子程序:定义优化问题的目标函数
def cal_Energy(X, nVar, mk): 	# m(k):惩罚因子,随迭代次数 k 逐渐增大p1 = (max(0, 6*X[0]+5*X[1]-60))**2p2 = (max(0, 10*X[0]+20*X[1]-150))**2fx = -(10*X[0]+9*X[1])return fx+mk*(p1+p2)# 子程序:模拟退火算法的参数设置
def ParameterSetting():cName = "funcOpt"           # 定义问题名称 YouCans, XUPTnVar = 2                    # 给定自变量数量,y=f(x1,..xn)xMin = [0, 0]               # 给定搜索空间的下限,x1_min,..xn_minxMax = [8, 8]               # 给定搜索空间的上限,x1_max,..xn_maxtInitial = 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 = TruekBetter += 1else:  # 容忍解:如果新解的目标函数比当前解差,则以一定概率接受新解pAccept = math.exp(-deltaE / tNow)  # 计算容忍解的状态迁移概率if pAccept > random.random():accept = True  # 接受劣质解kBadAccept += 1else:accept = False  # 拒绝劣质解kBadRefuse += 1# 保存新解if accept == True:  # 如果接受新解,则将新解保存为当前解xNow[:] = xNew[:]fxNow = fxNewif fxNew < fxBest:  # 如果新解的目标函数好于最优解,则将新解保存为最优解fxBest = fxNewxBest[:] = xNew[:]totalImprove += 1scale = 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 * alfakIter = kIter + 1fxBest = 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!")returnelse: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# = 关注 Youcans,分享原创系列 https://blog.csdn.net/youcans =
# 主程序
def main(): # YouCans, XUPT# 参数设置,优化问题参数定义,模拟退火算法参数设置[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、运行结果

x_Initial:2.000000,7.000000,	f(x_Initial):17.000000
i:0,t(i):100.00, badAccept:0.814286, f(x)_best:-152.000000
i:10,t(i):81.71, badAccept:0.635135, f(x)_best:-98.000000
i:20,t(i):66.76, badAccept:0.782051, f(x)_best:-98.000000
...
i:200,t(i):1.76, badAccept:0.090000, f(x)_best:-98.000000
i:210,t(i):1.44, badAccept:0.120000, f(x)_best:-98.000000
i:220,t(i):1.17, badAccept:0.130000, f(x)_best:-98.000000
improve:7Optimization by simulated annealing algorithm:x[0] = 8.0x[1] = 2.0f(x) = -98.0

参考文献:

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


版权说明:

= 关注 Youcans,分享原创系列 https://blog.csdn.net/youcans =
原创作品
Copyright 2021 YouCans, XUPT
Crated:2021-05-01

关注 Youcans,分享原创系列 https://blog.csdn.net/youcans

Python数模笔记-PuLP库(1)线性规划入门
Python数模笔记-PuLP库(2)线性规划进阶
Python数模笔记-PuLP库(3)线性规划实例
Python数模笔记-StatsModels 统计回归(1)简介
Python数模笔记-StatsModels 统计回归(2)线性回归
Python数模笔记-StatsModels 统计回归(3)模型数据的准备
Python数模笔记-StatsModels 统计回归(4)可视化
Python数模笔记-Sklearn (1)介绍
Python数模笔记-Sklearn (2)聚类分析
Python数模笔记-Sklearn (3)主成分分析
Python数模笔记-Sklearn (4)线性回归
Python数模笔记-Sklearn (5)支持向量机
Python数模笔记-模拟退火算法(1)多变量函数优化
Python数模笔记-模拟退火算法(2)约束条件的处理
Python数模笔记-模拟退火算法(3)整数规划问题
Python数模笔记-模拟退火算法(4)旅行商问题

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/news/566178.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

数据结构之算法特性及分类

数据结构之算法特性及分类 算法的特性 1.通用性。2.有效性。3.确定性4.有穷性。基本算法分类 1.穷举法顺序查找K值2.回溯,搜索八皇后&#xff0c;树和图遍历3.递归分治二分查找K值&#xff0c;快速排序&#xff0c;归并排序。4.贪心法Huffman编码树&#xff0c;最短路Dijkstra…

Python数模笔记-模拟退火算法(4)旅行商问题

1、旅行商问题(Travelling salesman problem, TSP) 旅行商问题是经典的组合优化问题&#xff0c;要求找到遍历所有城市且每个城市只访问一次的最短旅行路线&#xff0c;即对给定的正权完全图求其总权重最小的Hamilton回路&#xff1a;设有 n个城市和距离矩阵 D[dij]&#xff0…

神经网络概述

神经网络概述 以监督学习为例&#xff0c;假设我们有训练样本集 &#xff0c;那么神经网络算法能够提供一种复杂且非线性的假设模型 &#xff0c;它具有参数 &#xff0c;可以以此参数来拟合我们的数据。 为了描述神经网络&#xff0c;我们先从最简单的神经网络讲起&#x…

Python数模笔记-StatsModels 统计回归(1)简介

1、关于 StatsModels statsmodels&#xff08;http://www.statsmodels.org&#xff09;是一个Python库&#xff0c;用于拟合多种统计模型&#xff0c;执行统计测试以及数据探索和可视化。 欢迎关注 Youcans 原创系列&#xff0c;每周更新数模笔记 Python数模笔记-PuLP库 Pyth…

Python数模笔记-StatsModels 统计回归(2)线性回归

1、背景知识 1.1 插值、拟合、回归和预测 插值、拟合、回归和预测&#xff0c;都是数学建模中经常提到的概念&#xff0c;而且经常会被混为一谈。 插值&#xff0c;是在离散数据的基础上补插连续函数&#xff0c;使得这条连续曲线通过全部给定的离散数据点。 插值是离散函数…

Python数模笔记-StatsModels 统计回归(3)模型数据的准备

1、读取数据文件 回归分析问题所用的数据都是保存在数据文件中的&#xff0c;首先就要从数据文件读取数据。 数据文件的格式很多&#xff0c;最常用的是 .csv&#xff0c;.xls 和 .txt 文件&#xff0c;以及 sql 数据库文件的读取 。 欢迎关注 Youcans 原创系列&#xff0c;每…

神经网络反向传导算法

假设我们有一个固定样本集 &#xff0c;它包含 个样例。我们可以用批量梯度下降法来求解神经网络。具体来讲&#xff0c;对于单个样例 &#xff0c;其代价函数为&#xff1a; 这是一个&#xff08;二分之一的&#xff09;方差代价函数。给定一个包含 个样例的数据集&#xff…

Python数模笔记-StatsModels 统计回归(4)可视化

1、如何认识可视化&#xff1f; 图形总是比数据更加醒目、直观。解决统计回归问题&#xff0c;无论在分析问题的过程中&#xff0c;还是在结果的呈现和发表时&#xff0c;都需要可视化工具的帮助和支持。  欢迎关注 Youcans 原创系列&#xff0c;每周更新数模笔记 Python数…

梯度检验与高级优化

众所周知&#xff0c;反向传播算法很难调试得到正确结果&#xff0c;尤其是当实现程序存在很多难于发现的bug时。举例来说&#xff0c;索引的缺位错误&#xff08;off-by-one error&#xff09;会导致只有部分层的权重得到训练&#xff0c;再比如忘记计算偏置项。这些错误会使你…

Python数模笔记-Sklearn (1)介绍

1、SKlearn 是什么 Sklearn&#xff08;全称 SciKit-Learn&#xff09;&#xff0c;是基于 Python 语言的机器学习工具包。 Sklearn 主要用Python编写&#xff0c;建立在 Numpy、Scipy、Pandas 和 Matplotlib 的基础上&#xff0c;也用 Cython编写了一些核心算法来提高性能。…

自编码算法与稀疏性

目前为止&#xff0c;我们已经讨论了神经网络在有监督学习中的应用。在有监督学习中&#xff0c;训练样本是有类别标签的。现在假设我们只有一个没有带类别标签的训练样本集合 &#xff0c;其中 。自编码神经网络是一种无监督学习算法&#xff0c;它使用了反向传播算法&#…

Python数模笔记-Sklearn(2)聚类分析

1、分类的分类 分类的分类&#xff1f;没错&#xff0c;分类也有不同的种类&#xff0c;而且在数学建模、机器学习领域常常被混淆。 首先我们谈谈有监督学习&#xff08;Supervised learning&#xff09;和无监督学习&#xff08;Unsupervised learning&#xff09;&#xff…

可视化自编码器训练结果

训练完&#xff08;稀疏&#xff09;自编码器&#xff0c;我们还想把这自编码器学到的函数可视化出来&#xff0c;好弄明白它到底学到了什么。我们以在1010图像&#xff08;即n100&#xff09;上训练自编码器为例。在该自编码器中&#xff0c;每个隐藏单元i对如下关于输入的函数…

Python数模笔记-Sklearn(3)主成分分析

主成分分析&#xff08;Principal Components Analysis&#xff0c;PCA&#xff09;是一种数据降维技术&#xff0c;通过正交变换将一组相关性高的变量转换为较少的彼此独立、互不相关的变量&#xff0c;从而减少数据的维数。 1、数据降维 1.1 为什么要进行数据降维&#xff1…

稀疏自编码器一览表

下面是我们在推导sparse autoencoder时使用的符号一览表&#xff1a; 符号含义训练样本的输入特征&#xff0c;.输出值/目标值. 这里 可以是向量. 在autoencoder中&#xff0c;.第 个训练样本输入为 时的假设输出&#xff0c;其中包含参数 . 该输出应当与目标值 具有相同的…

Python数模笔记-Sklearn(4)线性回归

1、什么是线性回归&#xff1f; 回归分析&#xff08;Regression analysis)是一种统计分析方法&#xff0c;研究自变量和因变量之间的定量关系。回归分析不仅包括建立数学模型并估计模型参数&#xff0c;检验数学模型的可信度&#xff0c;也包括利用建立的模型和估计的模型参数…

Python数模笔记-Sklearn(5)支持向量机

支持向量机&#xff08;Support vector machine, SVM&#xff09;是一种二分类模型&#xff0c;是按有监督学习方式对数据进行二元分类的广义线性分类器。 支持向量机经常应用于模式识别问题&#xff0c;如人像识别、文本分类、手写识别、生物信息识别等领域。 1、支持向量机&…

矢量化编程

当使用学习算法时&#xff0c;一段更快的代码通常意味着项目进展更快。例如&#xff0c;如果你的学习算法需要花费20分钟运行完成&#xff0c;这意味着你每个小时能“尝试”3个新主意。但是假如你的程序需要20个小时来运行&#xff0c;这意味着你一天只能“尝试”一个新主意&am…

Python数模笔记-NetworkX(1)图的操作

1、NetworkX 图论与网络工具包 NetworkX 是基于 Python 语言的图论与复杂网络工具包&#xff0c;用于创建、操作和研究复杂网络的结构、动力学和功能。 NetworkX 可以以标准和非标准的数据格式描述图与网络&#xff0c;生成图与网络&#xff0c;分析网络结构&#xff0c;构建…