模拟退火算法（惩罚函数法求约束优化问题）

构造惩罚函数：

p1 = (max(0, 4642 - X[0] - X[1] - X[2] - X[3] - X[4] - X[5] - X[6] - X[7] - X[8] - X[9])) ** 2

p2 = (max(0, 23833 - X[10] - X[11] - X[12] - X[13] - X[14] - X[15] - X[16] - X[17])) ** 2

p3 = (max(0, 8155 - X[18] - X[19] - X[20] - X[21] - X[22] - X[23] - X[24] - X[25] - X[26] - X[27] - X[28] - X[29])) ** 2

若存在等式约束，例如x1+2*x2=m，也可以转化为惩罚函数：
p4=(x1+2*x2-m)**2

P(x)=p1+p2+p3+......

构造增广目标函数
L(x,m(k))=min(fx)+m(k)*P(x)

m(k)：惩罚因子，随外循环迭代次数k逐渐增大，但在内循环中保持不变

#  -*- 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, 4642 - X[0] - X[1] - X[2] - X[3] - X[4] - X[5] - X[6] - X[7] - X[8] - X[9])) ** 2p2 = (max(0, 23833 - X[10] - X[11] - X[12] - X[13] - X[14] - X[15] - X[16] - X[17])) ** 2p3 = (max(0, 8155 - X[18] - X[19] - X[20] - X[21] - X[22] - X[23] - X[24] - X[25] - X[26] - X[27] - X[28] - X[29])) ** 2fx = (X[0] + X[1] + X[2] + X[3] + X[4] + X[5] + X[6] + X[7] + X[8] + X[9]) * 1.2 + (X[10] + X[11] + X[12] + X[13] + X[14] + X[15] + X[16] + X[17]) * 1.1 + \(X[18] + X[19] + X[20] + X[21] + X[22] + X[23] + X[24] + X[25] + X[26] + X[27] + X[28] + X[29])return fx + mk * (p1 + p2 + p3)# 子程序：模拟退火算法的参数设置
def ParameterSetting():cName = "funcOpt"  # 定义问题名称 YouCans, XUPTnVar = 30  # 给定自变量数量，y=f(x1,..xn)# xMin = [0, 0]  # 给定搜索空间的下限，x1_min,..xn_min# xMax = [8, 7.5]  # 给定搜索空间的上限，x1_max,..xn_maxxMin = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]xMax = [3147, 30977, 1724, 966, 971, 7661, 9385, 2521, 699, 36972, 7885, 10207, 1181, 9768, 8181, 1014, 21293,2081, 2816, 21267, 1788, 736, 922, 595, 15114, 23695, 5398, 342, 2005, 381]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):# random.uniform(min,max) 在 [min,max] 范围内随机生成一个实数xInitial[v] = random.uniform(xMin[v], xMax[v])# 调用子函数 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] = xNow[v] + scale * (xMax[v] - xMin[v]) * random.normalvariate(0, 1)# random.normalvariate(0, 1)：产生服从均值为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[{}] = {:.6f}'.format(i, xBest[i]))print('\n\tf(x):{:.6f}'.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)# = 关注 Youcans，分享原创系列 https://blog.csdn.net/youcans =
if __name__ == '__main__':main()

代码改动自：

(1条消息) Python数模笔记-模拟退火算法（2）约束条件的处理_YouCans的博客-CSDN博客_有约束条件的模拟退火算法

模拟退火算法（惩罚函数法求约束优化问题）相关推荐

基于遗传算法和模拟退火算法改进的混合模拟退火算法（解决求函数极值问题，MATLAB代码已实现）
基本思想: 混合模拟退火算法时遗传算法和模拟退火算法的结合,在混合模拟退火算法中使用了大量的样本作为问题的可能解决方案而不是将单个样本作为一个问题的可能解决方案.对遗传算法中适应的概念进行相应改进. ...
Python数模笔记-模拟退火算法（3）整数规划问题
1.整数规划问题整数规划问题在工业.经济.国防.医疗等各行各业应用十分广泛,是指规划中的变量(全部或部分)限制为整数,属于离散优化问题(Discrete Optimization). 线性规划问题的 ...
Python数模笔记-模拟退火算法（2）约束条件的处理
1.最优化与线性规划最优化问题的三要素是决策变量.目标函数和约束条件. 线性规划(Linear programming),是研究线性约束条件下线性目标函数的极值问题的优化方法,常用于解决利用现有的资 ...
模拟退火算法简单理解
模拟退火算法算法流程图 1. 引言模拟退火算法(Simulate Anneal,SA)是一种通用概率演算法,用来在一个大的搜寻空间内找寻命题的最优解 .模拟退火是由S.Kirkpatrick, ...
matlab 约束函数,【优化求解】MATLAB约束优化之惩罚函数法
一.算法原理之前我们了解过的算法大部分都是无约束优化问题,其算法有:黄金分割法,牛顿法,拟牛顿法,共轭梯度法,单纯性法等.但在实际工程问题中,大多数优化问题都属于有约束优化问题.惩罚函数法就可以将约 ...
用模拟退火算法（simulated annealing / SA）求函数最小值
#用模拟退火算法(simulated annealing / SA)求函数最小值已知函数 f(x) = (x-1)^2 + 3,是求其在[ -2,2 ]的最小值时刻的解下面为运用模拟退火算法求解上 ...
java模拟退火算法求函数_模拟退火算法从原理到实战【基础篇】
模拟退火算法来源于固体退火原理,将固体加温至充分高,再让其徐徐冷却,加温时,固体内部粒子随温升变为无序状,内能增大,而徐徐冷却时粒子渐趋有序,在每个温度都达到平衡态,最后在常温时达到基态,内能减为最小 ...
模拟退火算法求函数最大、小值——python实现
模拟退火算法(Simulate Anneal,SA)是一种通用概率演算法,用来在一个大的搜寻空间内找寻命题的最优解.模拟退火是由S.Kirkpatrick, C.D.Gelatt和M.P.Vecchi ...
模拟退火算法求函数极值（含MATLAB代码实现）
二.模拟退火算法 1. 简介模拟退火算法的思想借鉴于固体的退火过程,当固体的温度很高时,内能比较大,固体内的粒子处于快速无序运动状态,当温度慢慢降低,固体的内能减小,粒子逐渐趋于有序,最终 ...

模拟退火算法（惩罚函数法求约束优化问题）

模拟退火算法（惩罚函数法求约束优化问题）相关推荐

最新文章

热门文章