带约束条件的运筹规划问题求解（模拟退火算法实现）

转载作者：我是一只小鸟更新时间：2023-04-19 06:31:16

0. 写在前面

超级简单的模拟退火算法实现ε٩(๑> ₃ <)۶з搭配最简单的线性规划模型进行讲解！但是如果需要的话可以直接修改编程非线性问题哦(´つヮ⊂︎) 。

1. 模型描述及处理

1.1 线性规划模型

\[max\,f(x)=10x_1+9x_2 \]

\(s.t.\) 。

\[6x_1+5x_2\leq{60}\tag{1} \]

\[10x_1+20x_2\leq{150}\tag{2} \]

\[0\leq{x_1}\leq{8}\tag{3} \]

\[0\leq{x_2}\leq{8}\tag{4} \]

1.2 引入惩罚函数处理模型

对约束条件引入惩罚函数:

对约束条件(1)，惩罚函数为： \(p_1=max(0,6x_1+5x_2-60)^2\) 。
对约束条件(2)，惩罚函数为： \(p_2=max(0,10x_1+20x_2-150)^2\) 。

那么，该问题的惩罚函数可以表示为:

\[P(x)=p_1+p_2 \]

由此，可将该问题的约束条件放入目标函数中，此时模型变为:

\[min\,g(x)=-(10x_1+9x_2)+P(x)\quad\forall{x_1,x_2}\in{[0,8]} \]

2. 程序实现

                        
                          # 模拟退火算法 程序：求解线性规划问题（整数规划）
# 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))**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(): # 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()

输出结果:

                        
                          x_Initial:0.000000,4.000000,	f(x_Initial):-36.000000
i:0,t(i):100.00, badAccept:0.925373, f(x)_best:-152.000000
i:10,t(i):81.71, badAccept:0.671053, f(x)_best:-98.000000
i:20,t(i):66.76, badAccept:0.722892, f(x)_best:-98.000000
i:30,t(i):54.55, badAccept:0.704225, f(x)_best:-98.000000
i:40,t(i):44.57, badAccept:0.542169, f(x)_best:-98.000000
i:50,t(i):36.42, badAccept:0.435294, f(x)_best:-98.000000
i:60,t(i):29.76, badAccept:0.359551, f(x)_best:-98.000000
i:70,t(i):24.31, badAccept:0.717647, f(x)_best:-98.000000
i:80,t(i):19.86, badAccept:0.388235, f(x)_best:-98.000000
i:90,t(i):16.23, badAccept:0.555556, f(x)_best:-98.000000
i:100,t(i):13.26, badAccept:0.482353, f(x)_best:-98.000000
i:110,t(i):10.84, badAccept:0.527473, f(x)_best:-98.000000
i:120,t(i):8.85, badAccept:0.164948, f(x)_best:-98.000000
i:130,t(i):7.23, badAccept:0.305263, f(x)_best:-98.000000
i:140,t(i):5.91, badAccept:0.120000, f(x)_best:-98.000000
i:150,t(i):4.83, badAccept:0.422680, f(x)_best:-98.000000
i:160,t(i):3.95, badAccept:0.111111, f(x)_best:-98.000000
i:170,t(i):3.22, badAccept:0.350000, f(x)_best:-98.000000
i:180,t(i):2.63, badAccept:0.280000, f(x)_best:-98.000000
i:190,t(i):2.15, badAccept:0.310000, f(x)_best:-98.000000
i:200,t(i):1.76, badAccept:0.390000, f(x)_best:-98.000000
i:210,t(i):1.44, badAccept:0.390000, f(x)_best:-98.000000
i:220,t(i):1.17, badAccept:0.380000, f(x)_best:-98.000000
improve:10

Optimization by simulated annealing algorithm:
	x[0] = 8.0
	x[1] = 2.0

	f(x) = -98.0

最后此篇关于带约束条件的运筹规划问题求解（模拟退火算法实现）的文章就讲到这里了,如果你想了解更多关于带约束条件的运筹规划问题求解（模拟退火算法实现）的内容请搜索CFSDN的文章或继续浏览相关文章，希望大家以后支持我的博客！。

文章推荐： OI数论中的上界估计与时间复杂度证明

文章推荐：项目讲解之常见安全漏洞

文章推荐： Elasticsearch搜索功能的实现（五）--实战

算法~利用zset实现滑动窗口限流
滑动窗口限流滑动窗口限流是一种常用的限流算法，通过维护一个固定大小的窗口，在单位时间内允许通过的请求次数不超过设定的阈值。具体来说，滑动窗口限流算法通常包括以下几个步骤：初始化：设置窗口
【算法】表达式求值
表达式求值：一个只有+,-,*,/的表达式，没有括号一种神奇的做法：使用数组存储数字和运算符，先把优先级别高的乘法和除法计算出来，再计算加法和减法 int GetVal(string s){
【算法】前缀和
【算法】前缀和题目先来看一道题目：（前缀和模板题）已知一个数组A[]，现在想要求出其中一些数字的和。输入格式：先是整数N,M，表示一共有N个数字，有M组询问接下来有N个数，表示A[1]..
【算法】二叉树的各种遍历方式
1.前序遍历根-左-右的顺序遍历，可以使用递归 void preOrder(Node *u){ if(u==NULL)return; printf("%d ",u->val);
【算法】01背包
先看题目物品不能分隔，必须全部取走或者留下，因此称为01背包（只有不取和取两种状态）看第一个样例我们需要把4个物品装入一个容量为10的背包我们可以简化问题，从小到大入手分析 weightva
算法 - 矩阵中被另一种颜色包围的颜色
我最近在一次采访中遇到了这个问题: 给出以下矩阵: [[ R R R R R R], [ R B B B R R], [ B R R R B B], [ R B R R R R]] 找出是否有任
使用Outlook发送电子邮件的C++算法
我正在尝试通过 C++ 算法从我的 outlook 帐户发送一封电子邮件，该帐户已经打开并记录，但真的不知道从哪里开始(对于 outlook-c++ 集成)，谷歌也没有帮我这么多。任何提示将不胜感激。
容器上滑动窗口的C++算法
我发现自己像这样编写了一个手工制作的 while 循环: std::list foo; // In my case, map, but list is simpler auto currentPoin
检测正方形后运行命令的c++算法
我有用于检测正方形的 opencv 代码。现在我想在检测正方形后，代码运行另一个命令。代码如下: #include "cv.h" #include "cxcore.h" #include "high
二值图像的泛洪填充C++算法
我正在尝试模拟一个 matlab 函数“imfill”来填充二进制图像(1 和 0 的二维矩阵)。我想在矩阵中指定一个起点，并像 imfill 的 4 连接版本那样进行洪水填充。这是否已经存在于
算法递归公式
我正在阅读 Robert Sedgewick 的《C++ 算法》。 Basic recurrences section it was mentioned as 这种循环出现在循环输入以消除一个项目的递
算法 - 如何生成日期结构？
我正在思考如何在我的日历中生成代表任务的数据结构(仅供我个人使用)。我有来自 DBMS 的按日期排序的任务记录，如下所示: 买牛奶(18.1.2013) 任务日期 (2013-01-15) 任务标签(
算法:查找恰好出现两次的元素
输入一个未排序的整数数组A[1..n]只有 O(d) :(d int) 计算每个元素在单次迭代中出现在列表中的次数。 map 是balanced Binary Search Tree基于确保 O(nl
算法——基于寻找最大匹配数
我遇到了一个问题，但我仍然不知道如何解决。我想出了如何用蛮力的方式来做到这一点，但是当有成千上万的元素时它就不起作用了。 Problem: Say you are given the followin
算法 - 用于计算成对相互出现的次数
我有一个列表列表。 L1= [[...][...][.......].......]如果我在展平列表后获取所有元素并从中提取唯一值，那么我会得到一个列表 L2。我有另一个列表 L3，它是 L2 的某个
算法 - 在矩阵中求和
我们得到二维矩阵数组(假设长度为 i 和宽度为 j)和整数 k我们必须找到包含这个或更大总和的最小矩形的大小F.e k=7 4 1 1 1 1 1 4 4 Anwser是2，因为4+4=8 >= 7，
算法:根据周数获取下一年日期工作类次类型
我实行 3 类倒制，每周换类。顺序为早类 (m)、晚类 (n) 和下午类 (a)。我固定的订单，即它永远不会改变，即使那个星期不工作也是如此。我创建了一个函数来获取 ISO 周数。当我给它一个日期时
算法 - 找到满足输入元素任意组合的所有集合
假设我们有一个输入，它是一个元素列表: {a, b, c, d, e, f} 还有不同的集合，可能包含这些元素的任意组合，也可能包含不在输入列表中的其他元素: A:{e,f} B:{d,f,a} C:
算法:添加新元素时如何找到集合的子集？
我有一个子集算法，可以找到给定集合的所有子集。原始集合的问题在于它是一个不断增长的集合，如果向其中添加元素，我需要再次重新计算它的子集。有没有一种方法可以优化子集算法，该算法可以从最后一个计算点重新
算法:按预期频率将符号压缩成位串？
我有一个包含 100 万个符号及其预期频率的表格。我想通过为每个符号分配一个唯一(且前缀唯一)的可变长度位串来压缩这些符号的序列，然后将它们连接在一起以表示序列。我想分配这些位串，以使编码序列的预

我是一只小鸟

个人简介

我是一名优秀的程序员,十分优秀！

作者热门文章

滴滴打车优惠券免费领取

全站热门文章

首页

博学

6Ren·AI