连续状态转移算法（STA）的实现（python版）

背景介绍

所有的随机优化算法（GA、ABC、PSO、STA、模拟退火）都是一个套路

通过当前的一个解或者一些产生一堆新的候选解，这一步叫产生候选解
然后通过某个评价函数（评价标准）更新解，就是找一堆里面好的一个或者一些，这叫更新当前解
不断循环，直到满足终止条件

算子

所谓的算子其实就是产生候选解的一种规则，你可以理解成，给算子一定输入，他就给你返回一个或者多个候选解出来

连续STA的算子说明

基本的算子有四个，分别是，旋转变换，伸缩变化，轴向变换、平移变换

旋转变换

利用旋转变换公式产生新解，物理上是对 Best 为中心 alpha 为半径的超球体内进行采样SE个解

伸缩变换

伸缩变换公式产生新解，物理上是对 Best 的每个维度都有可能伸缩到（-∞，+∞），然后进行约束到定义域边界，属于全局搜索

轴向变换

利用轴向变换公式产生新解，物理上是对 Best 的单个维度都有可能伸缩到（-∞，+∞），属于局部搜索，增强的是单维搜索能力

平移变换

利用平移公式产生新解，物理上是在 oldBest —— newBest 这条直线上进行搜索产生新解，我们认为新旧两个最优解的连线上大概率会出现好的解，比如两个解在谷底两侧时，就是上面的图只有一维进行缩放的情况

文件结构

代码

benchmark.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-from functools import reduce
from operator import mul
import numpy as np
import mathdef Sphere(s):square = map(lambda y: y * y ,s)return sum(square)def Rosenbrock(s):square = map(lambda x, y: 100 * (y - x**2)**2 + (x - 1)**2 , s[:len(s)-1],s[1:len(s)])return  sum(square)def Rastrigin(s):square = map(lambda x: x**2 - 10 * math.cos(2*math.pi*x) + 10, s)return  float(sum(square))def Michalewicz(s):n = len(s)t1 = -sum(map(lambda x, y: math.sin(x) * (math.sin( y*x**2/math.pi ))**20, s, range(1,n+1)))return t1def Griewank(s): t1 = sum(map(lambda x: 1/4000 * x**2, s))n = len(s)t2 = map(lambda x, y: math.cos(x/np.sqrt(y)), s, range(1,n+1))t3 = reduce(mul, t2)return t1 - t3 + 1

STA.py

# -*- coding: utf-8 -*-
"""
Created on Sat Nov  6 10:57:13 2021@author: 中南大学自动化学院  智能控制与优化决策课题组
"""import numpy as np
import timeclass STA():# 初始化def __init__(self,funfcn,Range,SE,Maxiter,dict_opt):self.funfcn = funfcnself.Range = Rangeself.SE = SEself.Maxiter = Maxiterself.dict_opt = dict_optself.Dim = self.Range.shape[1]def initialization(self):self.Best = np.array(self.Range[0,:] + self.Range[1,:]-self.Range[1,:]*np.random.uniform(0,1,(1,self.Range.shape[1])))self.fBest = self.funfcn(self.Best[0])self.history = []self.history.append(self.fBest)"""利用旋转变换公式产生新解，物理上是对 Best 为中心 alpha 为半径的超球体内进行采样SE个解"""def op_rotate(self):R1 = np.random.uniform(-1,1,(self.SE,self.Dim))R2 = np.random.uniform(-1,1,(self.SE,1))a = np.tile(self.Best,(self.SE,1))b = np.tile(R2,(1,self.Dim))c = R1/np.tile(np.linalg.norm(R1,axis=1,keepdims = True),(1,self.Dim))State = a + self.dict_opt['alpha']*b*creturn State"""利用伸缩变换公式产生新解，物理上是对 Best 的  每个  维度都有可能伸缩到（-∞，+∞），然后进行约束到定义域边界，属于全局搜索"""def op_expand(self): a = np.tile(self.Best,(self.SE,1))b = np.random.randn(self.SE,self.Dim)State = a + self.dict_opt['gamma'] * b * areturn State"""利用轴向变换公式产生新解，物理上是对 Best 的  单个 维度都有可能伸缩到（-∞，+∞），属于局部搜索，增强的是单维搜索能力"""def op_axes(self):#实现一State = np.tile(self.Best,(self.SE,1))for i in range(self.SE):index = np.random.randint(0,self.Dim)State[i,index] = State[i,index] + self.dict_opt['delta']*np.random.randn()*State[i,index]return State#实现二
#        a = np.tile(self.Best,(self.SE,1))
#        A = np.zeros((self.SE,self.Dim))
#        A[np.arange(self.SE),np.random.randint(0,self.Dim,self.SE)] = np.random.randn(self.SE)
#        State = a + A*a
#        return State"""利用平移公式产生新解，物理上是在 oldBest —— newBest 这条直线上进行搜索产生新解，我们认为新旧两个最优解的连线上大概率会出现好的解，比如两个解在谷底两侧时"""def op_translate(self,oldBest):a = np.tile(self.Best,(self.SE,1)) # SE * nb = np.random.uniform(0,1,(self.SE,1))c = self.dict_opt['beta']*(self.Best - oldBest)/(np.linalg.norm(self.Best - oldBest)) # 1*nState = a + b * creturn Statedef selection(self,State):fState = np.zeros((self.SE,1)) #for i in range(self.SE):fState[i] = self.funfcn(State[i,:])index = np.argmin(fState)return State[index,:],fState[index,:]def bound(self,State):Pop_Lb = np.tile(self.Range[0],(State.shape[0],1))Pop_Ub = np.tile(self.Range[1],(State.shape[0],1))changeRows = State > Pop_UbState[changeRows] = Pop_Ub[changeRows]changeRows = State < Pop_LbState[changeRows] = Pop_Lb[changeRows]return State       def run(self):time_start = time.time()self.initialization()for i in range(self.Maxiter):if self.dict_opt['alpha'] < 1e-4:self.dict_opt['alpha'] = 1.0else:self.dict_opt['alpha'] = 0.5*self.dict_opt['alpha']"""循环使用三个算子产生新解并且更新最优解"""dict_op = {"rotate": self.op_rotate(), "expand": self.op_expand(),"axes": self.op_axes()}for key in dict_op:oldBest = self.BestState = self.bound(dict_op[key])newBest, fnewBest = self.selection(State)if fnewBest < self.fBest:      # 如果算子（三个中任意一个）产生得解确实比历史得好，那再调用translate算子试试能不能找到更好得self.fBest, self.Best = fnewBest, newBestState = self.bound(self.op_translate(oldBest))newBest, fnewBest = self.selection(State)if fnewBest < self.fBest:self.fBest, self.Best = fnewBest, newBestself.history.append(self.fBest[0])print("第{}次迭代的最优值是:{:.5e}".format(i,self.fBest[0]))# """让算法提前停止,如果两次最优值的更新量小于 1e-10"""# if (abs(self.history[i+1] - self.history[i]) < 1e-10) & (self.dict_opt['alpha'] < 1e-4):#     self.history = np.array(self.history[1:])#     breaktime_end = time.time()print("总耗时为:{:.5f}".format(time_end - time_start))

Test.py

# -*- coding: utf-8 -*-
"""
Created on Sat Nov  6 10:57:13 2021@author: 中南大学自动化学院  智能控制与优化决策课题组
"""from STA import STA
import Benchmark
import numpy as np
import matplotlib.pyplot as plt#参数设置
funfcn = Benchmark.Sphere
Dim = 100
Range = np.repeat([-30,30],Dim).reshape(-1,Dim)
SE = 30
Maxiter = 1000
dict_opt = {'alpha':1,'beta':1,'gamma':1,'delta':1}
# 算法调用
sta = STA(funfcn,Range,SE,Maxiter,dict_opt)
sta.run()
# 画图
#x = np.arange(1,Maxiter+2).reshape(-1,1)
y = sta.history
x = np.arange(len(y)).reshape(-1,1)
plt.semilogy(x,y,'b-o')
plt.xlabel('Ierations')
plt.ylabel('fitness')
plt.show()

结果

这是 Sphere 函数的

后记

可以自行更换函数测试，更多的测试函数可以看无约束 benchmark 函数，通过编写函数，自行测试

连续状态转移算法（STA）的实现（python版）相关推荐

为什么只有状态转移算法才是真正意义上的智能优化算法，其它的都是“假冒伪劣”？
0 旁白长期以来,"论文为王"的观念被推崇至上,没有论文就没有发言权,发表了大量学术论文才能高谈阔论.坐而论道,以至于出现"大家都忙着写论文,没有时间搞科研" ...
python算法与程序设计基础第二版-算法与程序设计基础(Python版) - 吴萍
基本信息书名:21世纪高等学校计算机基础实用规划教材:算法与程序设计基础(Python版) 定价:39.00元作者:吴萍21世纪高出版社:清华大学出版社出版日期:2015_2_1 ISBN:9 ...
算法设计与分析(python版)-作业一
参考教材:算法设计与分析(Python版) 作者:王秋芬 1 . 容易 (4分)2 n=O(100n ^2) 错误 2 . 容易 (3分)10=θ(log10) 正确 3 . 容易 ( ...
算法设计与分析(python版)-作业三
参考教材:算法设计与分析(Python版) 作者:王秋芬 1 . 普通 (5分)以下问题中,哪些问题的分治算法消耗的时间与输入序列无关.() A. 二分查找 B. 合并排序 C. 快速 ...
算法刷题（蓝）【基础篇+算法篇】【Python版】
[前言] 记录自己在刷蓝桥杯题目的一些做题思路,在构思的过程中,会参考一些大佬的代码( 用到了会提供相应的学习链接). 内容有不恰当之处,请各位大佬们批评指正,我会第一时间进行更改. 语言:pytho ...
算法的时间复杂度（python版容易理解）+常用的时间复杂度、python代码--数据结构
https://blog.csdn.net/hanhanwanghaha宝藏女孩欢迎您的关注! 欢迎关注微信公众号:宝藏女孩的成长日记如有转载,请注明出处(如不注明,盗者必究) 目录一.算法时间 ...
啊哈算法-DFS解救小哈python版
DFS 啊哈算法-解救小哈 maze_map = [] n,m = map(int,(input()).split()) maze_map = [input().split() for i in ra ...
算法面试用c还是python_排序算法（C语言+Python版）宝宝再也不怕面试官写排序算法了...
直接插入排序过程: 1. 数据可分看成两个部分,前面的数据是有序的 2. 从后面的数据取出一个元素,插到前面有序数据的合适位置从右端开始查找,到找到比此元素大的时候,则此元素向后移动,以空出多余的 ...
Python版基于递归的冒泡排序算法
应读者要求,写个基于递归的冒泡排序算法代码,之前发过的排序算法代码请参考Python版快速排序算法,Python版选择排序算法,Python版冒泡法排序算法. from random import r ...