状态转移算法

1 相关资源：
　　2020年发表在自动化学报的论文《状态转移算法原理及应用》，对状态转移算法作了广泛而深入的介绍：
http://www.aas.net.cn/cn/article/doi/10.16383/j.aas.c190624
　　各种状态转移算法的源代码可以在中南大学周晓君教授的个人主页免费下载：
https://faculty.csu.edu.cn/michael_x_zhou/zh_CN/article/198862/content/2802.htm#article
  或者GitHub开源平台：
https://github.com/tiezhongyu2005/STA

2 基本状态转移算法

　　状态转移算法(state transition algorithm)是一种基于结构主义学习的智能优化算法，它抓住最优化算法的本质、目的和要求，以全局性、最优性、快速性、收敛性、可控性为核心结构要素，需要什么才学习什么。2012年，周晓君教授正式提出这种新颖的智能型随机性全局优化方法,它的基本思想是：将最优化问题的一个解看成是一个状态, 解的迭代更新过程看成是状态转移过程, 利用现代控制理论的状态空间表达式来作为产生候选解的统一框架, 基于此框架来设计状态变换算子。与大多数基于种群的进化算法不同, 基本的状态转移算法是一种基于个体的进化算法, 它基于给定当前解, 通过采样方式, 多次独立运行某种状态变换算子产生候选解集, 并与当前解进行比较, 迭代更新当前解，直到满足某种终止条件。 值得一提的是, 基本状态转移算法中的每种状态变换算子都能够产生具有规则形状、可控大小的几何邻域, 它设计了包括旋转变换、平移变换、伸缩变换、坐标轴搜索等不同的状态变换算子以满足全局搜索、局部搜索以及启发式搜索等功能需要, 并且以交替轮换的方式适时地使用各种不同算子, 使得状态转移算法能够以一定概率很快找到全局最优解。

2.1 状态转移算法的统一框架

　　借鉴于现代控制理论中的状态空间模型表示法，状态转移算法的基本框架可表述如下：
$$ \left\{\begin{array}{l} s_{k+1}=\boldsymbol{A}_{k} s_{k}+\boldsymbol{B}_{k} \boldsymbol{u}_{k} \\ y_{k+1}=f\left(s_{k+1}\right) \end{array}\right. $$

　　其中,s_k = (s₁, s₂, ···, s_n)^T为当前状态, 代表优化问题的一个候选解; 𝑨_k 或 𝑩_k 为状态转移矩阵, 可以看成是算法中算子所产生的的变换效果; 𝒖_k为当前状态及历史状态的函数, 可以看成是控制变量; f(·)为目标函数或者评价函数。

2.2 基本状态转移算法的理论实现

　　状态转移算法中专门设计了局部、全局和启发式搜索算子，其中全局搜索算子是为了保证有一定概率使产生的候选解集形成的邻域中包含全局最优解； 局部搜索算子具有在较小邻域进行精细搜索的功能;；启发式搜索算子用来产生具有一定潜在价值的候选解，避免搜索的盲目性。 另一方面，状态转移算法中设计了特定的选择与更新策略， 用来适应不同优化问题的需要。 更进一步地， 状态转移算法中实施了智能化策略， 比如全局和局部算子的交替调用，可以很大程度上缩短寻找全局最优解的时间和避免陷入局部最优； 采样方式的使用,可以选择有代表性的候选解， 避免穷举邻域内的所有解，将极大地缩短搜索时间。一般地，状态转移算法的基本流程如下:

　　i ) 从当前最好解出发，通过某种状态变换算子生成具有某种特定性质的邻域；

　　ii ) 以某种采样机制从该邻域内采集一定数目的样本；

　　iii ) 根据选择与更新策略，比较当前最好解与采集样本中最好解的优劣， 并以某种更新机制更新当前最好解；
iv ) 返回步骤i )， 并交替使用各种状态变换算子， 直到满足某种终止条件。

2.2.1 状态变换算子

　　（1）旋转变换
$$ s_{k+1}=s_{k}+\alpha \frac{1}{n\left\|s_{k}\right\|_{2}} \boldsymbol{R}_{r} s_{k} $$

其中 α 为旋转因子；R_r ∈ ℝ^n×n 为随机矩阵，其元素取值于 [-1,1] 之间均匀分布；‖·‖₂ 表示取向量的模。旋转变换能够以当前解 s_k 为原点，在给定半径 α 内进行局部搜索。

　　（2）平移变换
$$ s_{k+1}=x_{k}+\beta \boldsymbol{R}_{t} \frac{s_{k}-s_{k-1}}{\left\|s_{k}-s_{k-1}\right\|_{2}} $$

其中 β 为平移因子；R_t ∈ ℝ 是一个随机变量，其取值于 [0,1] 之间均匀分布。平移变换沿着当前解 s_k 和前一个历史解 s_k-1 所形成的直线上进行启发式搜索。

　　（3）伸缩变换
$$ s_{k+1}=s_{k}+\gamma \boldsymbol{R}_{e} s_{k} $$

其中 γ 为伸缩因子；R_e ∈ ℝ^n×n 是一个非零元素取值服从正态分布的随机对角矩阵。从式子中不难理解，伸缩变换能够以正态分布的概率，将当前解 s_k 的每个维度上的值伸缩到实轴上的任意位置，具有全局搜索的能力。

　　（4）轴向变换
$$ s_{k+1}=s_{k}+\delta \boldsymbol{R}_{a} s_{k} $$

其中 δ 为坐标因子；R_a ∈ ℝ^n×n 是一个非零元素取值服从正态分布的稀疏对角矩阵，在基本状态转移算法中，R_a 仅有一个位置为非零元素。坐标变换能够从当前解 s_k 出发，随机挑选一个维度进行搜索，能够增强单一维度的搜索能力。

2.2.2 邻域、采样、选择与更新

　　邻域：

　　由当前状态s_k出发，通过状态变换算子能够产生候选状态s_k+1。不难看出，由于随机性，某种特定的算子能够反复产生不相等但是具有同质性的一系列候选状态。这些候选解将构成一个集合，形成一个邻域。理论上，邻域中具有无数个候选解，但在实际中，必须采用一定的采样机制选出有限个候选解。　　

　　采样：

　　对于邻域中的具有同质性的候选解，以随机策略采样SE个候选解，产生一个大小有限且可控的集合State作为实际当中的状态候选集合。SE可理解为“搜索强度”、“采样力度”或者“样本个数”。

　　选择与更新：

　　在算法运行的阶段中，考虑State中的SE个候选解和当前最优解，采取“贪婪策略”来进行选择与更新最优解。每次选取具有最优适应度函数值的解来更新当前最优解。

2.2.3 交替轮换

　　为了实现可控、最优、全局、稳定的搜索，具有不同搜索功能的状态变换算子被交替轮换地调用。这种设计思路，从一开始就清晰地规划了优化中各种可能需要的搜索形式，并将其抽象封装为功能各异的状态变换算子。于是，在解决实际问题的时候，不同搜索功能可以有选择地被强化或者削弱，以适应不同场景的具体需要。进一步改进优化算法的思路也就从漫无目的地构思新的动物算法，变成了如何为已知的状态变换算子设计合适的调用策略。

2.2.4 参数设置

在基本的状态转移算子中，参数设置为：αmax = 1，αmin = 1e-4，β = 1，γ = 1，δ = 1，SE = 30，fc = 2。其中平移、伸缩和坐标因子的大小都保持恒定，旋转因子以 1/fc 为衰减率，在最大值 αmax 和最小值 αmin之间周期变化。

3 MATLAB源码描述与解析

如图为基本状态转移算法的MATLAB源码的全部文件。本部分将首先对各个M文件给出概括性的功能描述，然后针对具体代码给出给出解释与分析。

3.1 概括性的功能描述

　　Test_sta.m为算法的调用示例；STA.m为封装好的基本状态转移算法；initialization.m为解的初始化过程；fitness.m为适应度计算与比较；Get_Test_Functions.m为各类常见优化问题目标函数的汇总封装；rotate.m，expand.m，axesion.m为对应状态变换算子的外层逻辑；op_rotate.m，op_expand.m，op_axesion.m，op_translate.m分别是四种状态变换算子的矩阵运算实现。

3.2 源码解析

在这个部分，将给出基本状态转移算法的主要源码与解析。其中：rotate.m，expand.m，axesion.m的写法与功能类似，故仅以rotate.m为例说明；op_rotate.m，op_expand.m，op_axesion.m，op_translate.m的性质相同，故仅以op_rotate.m为例说明。

Test_sta.m ：

clear all
clc
currentFolder = pwd;
addpath(genpath(currentFolder))
% parameter setting
warning('off')
SE = 30;%搜索强度
Dim = 30;%问题维度定义
Range = repmat([-5.12;5.12],1,Dim);%搜索范围定义
Iterations = 1e3;%最大迭代次数定义
tic
[Best,fBest,history] = STA(@Rastrigin,SE,Dim,Range,Iterations);
%调用sta.m，以SE的搜索强度，迭代最多Iterations次，在搜索范围Range内，
%解决变量数为Dim个的典型优化问题Rastrigin函数。并返回最优解Best，
%最优值fBest以及搜索历史history。
toc
history
semilogy(history)

STA.m ：

%参数初始化
alpha_max = 1;
alpha_min = 1e-4;
alpha = alpha_max;
beta = 1;
gamma = 1;
delta = 1;
fc = 2;
% 初始化状态集合State
State = initialization(SE,Dim,Range);
%计算适应度函数，返回当前最优解Best和最优值fBest
[Best,fBest] = fitness(funfcn,State);
counter = 0;
%迭代过程
for iter = 1:Iterations
    if alpha < alpha_min
        alpha = alpha_max;
    end
    oldfBest = fBest;
    [Best,fBest] = expand(funfcn,Best,SE,Range,beta,gamma);    
    [Best,fBest] = rotate(funfcn,Best,SE,Range,alpha,beta);   
    [Best,fBest] = axesion(funfcn,Best,SE,Range,beta,delta);    
    %终止条件判断
    if norm(oldfBest-fBest) < 1e-8 % 计算矩阵范数表示两代最优值之间的误差
        counter = counter + 1;
        if counter > 10
            break;
        end
    else
        counter = 0;
    end
    fprintf('iter=%d      ObjVal=%g\n',iter,fBest);
    history(iter) = fBest;
    alpha = alpha/fc;%更新参数alpha的值
end

initialization.m ：

function  State=initialization(SE,Dim,Range)
Pop_Lb = Range(1,:);%Range的第一行，全是-5.12
Pop_Ub = Range(2,:);%Range的第二行，全是5.12
%生成初始状态集合State，大小为Dim*SE
State  = rand(SE,Dim).*repmat(Pop_Ub-Pop_Lb,SE,1)+repmat(Pop_Lb,SE,1);

fitness.m :

function [Best,fBest] = fitness(funfcn,State)
%计算适应度
SE = size(State,1);
fState = zeros(SE,1);
for i = 1:SE
    fState(i) = feval(funfcn,State(i,:));
end
[fGBest, g] = min(fState);%返回最小值fGBest和其索引g
fBest = fGBest;
Best = State(g,:);%最值对应那一行的解

op_rotate.m :

function y=op_rotate(Best,SE,alpha)
n = length(Best);
y = repmat(Best',1,SE) + alpha*(1/n/(norm(Best)+eps))*reshape(unifrnd(-1,1,SE*n,n)*Best',n,SE);%旋转算子的矩阵运算实现
y = y';%转置