背景介绍
所有的随机优化算法（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 math
def 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 t1
def 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 time
class STA():
# 初始化
def __init__(self,funfcn,Range,SE,Maxiter,dict_opt):
self.funfcn = funfcn
self.Range = Range
self.SE = SE
self.Maxiter = Maxiter
self.dict_opt = dict_opt
self.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*c
return 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 * a
return 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 * n
b = np.random.uniform(0,1,(self.SE,1))
c = self.dict_opt['beta']*(self.Best - oldBest)/(np.linalg.norm(self.Best - oldBest)) # 1*n
State = a + b * c
return State
def 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_Ub
State[changeRows] = Pop_Ub[changeRows]
changeRows = State < Pop_Lb
State[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.0
else:
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.Best
State = self.bound(dict_op[key])
newBest, fnewBest = self.selection(State)
if fnewBest < self.fBest: # 如果算子（三个中任意一个）产生得解确实比历史得好，那再调用translate算子试试能不能找到更好得
self.fBest, self.Best = fnewBest, newBest
State = self.bound(self.op_translate(oldBest))
newBest, fnewBest = self.selection(State)
if fnewBest < self.fBest:
self.fBest, self.Best = fnewBest, newBest
self.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:])
# break
time_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 函数 ，通过编写函数，自行测试