$$
\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_{1}, s_{2}, · · · , s_{n})^{T}$为当前状态, 代表优化问题的一个候选解; $\boldsymbol{A}_{k}$ 或$\boldsymbol{B}_{k}$ 为状态转移矩阵, 可以看成是算法中算子所产生的的变换效果; $\boldsymbol{u}_{k}$ 为当前状态及历史状态的函数, 可以看成是控制变量; $f(·)$为目标函数或者评价函数。
2.2 基本状态转移算法的理论实现
状态转移算法中专门设计了局部、全局和启发式搜索算子，其中全局搜索算子是为了保证有一定概率使产生的候选解集形成的邻域中包含全局最优解； 局部搜索算子具有在较小邻域进行精细搜索的功能;；启发式搜索算子用来产生具有一定潜在价值的候选解，避免搜索的盲目性。 另一方面，状态转移算法中设计了特定的选择与更新策略， 用来适应不同优化问题的需要。 更进一步地， 状态转移算法中实施了智能化策略， 比如全局和局部算子的交替调用，可以很大程度上缩短寻找全局最优解的时间和避免陷入局部最优； 采样方式的使用,可以选择有代表性的候选解， 避免穷举邻域内的所有解，将极大地缩短搜索时间。一般地，状态转移算法的基本流程如下:
i ) 从当前最好解出发，通过某种状态变换算子生成具有某种特定性质的邻域；
ii ) 以某种采样机制从该邻域内采集一定数目的样本；
iii ) 根据选择与更新策略，比较当前最好解与采集样本中最好解的优劣， 并以某种更新机制更新当前最好解； iv ) 返回步骤i )， 并交替使用各种状态变换算子， 直到满足某种终止条件。
$$
s_{k+1}=s_{k}+\alpha \frac{1}{n\left\|s_{k}\right\|_{2}} \boldsymbol{R}_{r} s_{k}{(2)}
$$
其中$\alpha$为旋转因子；$\boldsymbol{R}_{r} \in \mathbb{R}^{n \times n}$为随机矩阵，其元素取值于[-1,1]之间均匀分布；$\|\cdot\|_{2}$表示取向量的模。旋转变换能够以当前解$s_{k}$为原点，在给定半径$\alpha$内进行局部搜索。
$$
s_{k+1}=x_{k}+\beta \boldsymbol{R}_{t} \frac{s_{k}-s_{k-1}}{\left\|s_{k}-s_{k-1}\right\|_{2}}
$$
其中$\beta$为平移因子；$\boldsymbol{R}_{t} \in \mathbb{R}^{n \times n}$是一个随机向量，其元素取值于[0,1]之间均匀分布。平移变换沿着当前解$s_{k}$和前一个历史解$s_{k-1}$所形成的的直线上进行启发式搜索。
$$
s_{k+1}=s_{k}+\gamma \boldsymbol{R}_{e} s_{k}
$$
其中$\gamma$为伸缩因子；$\boldsymbol{R}_{e} \in \mathbb{R}^{n \times n}$是一个非零元素取值服从正态分布的随机对角矩阵。从式子中不难理解，伸缩变换能够以正态分布的概率，将当前解$s_{k}$的每个维度上的值伸缩到实轴上的任意位置，具有全局搜索的能力。
$$
s_{k+1}=s_{k}+\delta \boldsymbol{R}_{a} s_{k}
$$
其中$\delta$为坐标因子；$\boldsymbol{R}_{a} \in \mathbb{R}^{n \times n}$是一个非零元素取值服从正态分布的稀疏对角矩阵，在基本状态转移算法中，$\boldsymbol{R}_{a}$仅有一个位置为非零元素。坐标变换能够从当前解$s_{k}$出发，随机挑选一个维度进行搜索，能够增强单一维度的搜索能力。
2.2.2 邻域、采样、选择与更新
邻域：
由当前状态$s_{k}$出发，通过状态变换算子能够产生候选状态$s_{k+1}$。不难看出，由于随机性，某种特定的算子能够反复产生不相等但是具有同质性的一系列候选状态。这些候选解将构成一个集合，形成一个邻域。理论上，邻域中具有无数个候选解，但在实际中，必须采用一定的采样机制选出有限个候选解。
采样：
对于邻域中的具有同质性的候选解，以随机策略采样$SE$个候选解，产生一个大小有限且可控的集合$State$作为实际当中的状态候选集合。$SE$可理解为“搜索强度”、“采样力度”或者“样本个数”。
选择与更新：
在算法运行的阶段中，考虑$State$中的$SE$个候选解和当前最优解，采取“贪婪策略”来进行选择与更新最优解。每次选取具有最优适应度函数值的解来更新当前最优解。
2.2.3 交替轮换
为了实现可控、最优、全局、稳定的搜索，具有不同搜索功能的状态变换算子被交替轮换地调用。这种设计思路，从一开始就清晰地规划了优化中各种可能需要的搜索形式，并将其抽象封装为功能各异的状态变换算子。于是，在解决实际问题的时候，不同搜索功能可以有选择地被强化或者削弱，以适应不同场景的具体需要。进一步改进优化算法的思路也就从漫无目的地构思新的动物算法，变成了如何为已知的状态变换算子设计合适的调用策略。
2.2.4 参数设置
在基本的状态转移算子中，参数设置为：$\alpha_{max}=1$，$\alpha_{min}=1e-4$，$\beta=1$，$\gamma=1$，$\delta=1$，$SE=30$，$fc=2$。其中平移、伸缩和坐标因子的大小都保持恒定，旋转因子以$1/fc$为衰减率，在最大值$\alpha_{max}$和最小值$\alpha_{min}$之间周期变化。
如图为基本状态转移算法的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 源码解析
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] = fitness(funfcn,State);%计算适应度函数，返回当前最优解Best和最优值fBest
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 % can be changed 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 = rand(SE,Dim).*repmat(Pop_Ub-Pop_Lb,SE,1)+repmat(Pop_Lb,SE,1);%生成初始状态集合State，大小为Dim*SE
fitness.m :
function [Best,fBest] = fitness(funfcn,State)
% calculate fitness
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';%转置