模拟退火算法就是在搜索过程中，不仅接受更优解，而且在满足一定概率的情况下接受劣解。接受更优解使解不断收敛，接受劣解有助于跳出某个局部最优解。
模拟退火算法
Part 1 模拟退火算法的自然界原型
Point 1 退火
顾名思义，模拟退火算法就是在模拟退火。
说到“退火”这个名词，机械类和材料类的学生一定记得那些年我们拿砂纸打磨过的小铁块。
退火是指将固体加热到足够高的温度，使分子呈随机排列状态，然后逐步降温使之冷却，最后分子以低能状态排列，固体达到某种稳定状态。
在初始状态，固体内部的粒子存在非均匀状态。
加热到一定温度，增强其粒子的热运动，使粒子偏离其平衡位置，消除原来的非均匀状态。
退火过程中要让温度慢慢降低，在每一个温度下要达到热平衡状态，粒子逐渐均匀有序，最终达到内能最小的状态。
Point 2 Metropolis准则
在退火过程中，某一个温度为T的时刻，固体从一个状态转变到另一个状态遵循Metropolis准则，即以概率接受新状态。Metropolis准则如下：
在某一个温度T时
若新状态j的内能 
 小于当前状态的内能
，即 
 ，则直接接受j为当前状态
否则，若 
 ，则仍然接受j为当前状态，如果不成立，则保留当前状态i。 
也就是说，在某个温度下固体分子从一个状态转移到另一个状态时，如果新状态的内能小，则无条件接受；如果新状态的内能大，则满足判断条件再接受。
其中判断条件与状态j的内能与状态i的内能差距 
 有关， 
 比 
 大的越多，则 
 越小，
 越小，rand小于p的概率越小，则接受的概率越小。
并且判断条件还跟温度T的大小有关系，温度T越小，则 
 越小， 
 越小，rand小于p的概率越小，则接受的概率越小。温度高时更容易接受差的状态，温度低时更不容易接受差的状态。
Point 3 关于温度T
正是由于Metropolis准则中温度T变化的特性，T的选择对于最后的结果有很大影响。降温步骤太快，降低了粒子往其他方向运动的可能性，容易陷入局部最优；降温速度过慢，又会使整个退火时间过长。
Part 2 由自然界原型到模拟退火算法
Point 1 相似点
固体退火与优化问题有着一定的相似性：
根据这种相似性，并依据Metropolis准则取舍进行迭代就形成了模拟退火算法。
Point 2 具体实现
模拟退火算法的实现主要有四个步骤：
在当前温度（状态）下产生一个位于解空间的新解x1
计算新解x1的函数值
判断是否接受新解，判断依据为函数值下降及Metropolis准则。；如果不接受则重新产生一个x1，直到接受x1为止
如果接受新解，则在本轮的基础上进行下一次迭代
基本模拟退火算法流程图如下：
在模拟退火方法的应用过程中，主要是围绕着以下几个方面对其进行调整：
如何产生在解空间里的新解
如何应用Metropolis准则（作出一定的改进）
如何缩减温度
Part 3 模拟退火方法解决TSP问题
在TSP问题中，粒子状态就是城市的排列顺序，内能就是路径的长度。
用交换某两个随机位置的城市顺序来实现在原有解的基础上变异产生新解。
直接使用原有的Metropolis准则
使用 
 来缩减温度，其中a为一个小于1的数
以求10个城市的TSP问题为例，10个城市的坐标如下：
%----------求解TSP问题----------%
clc;
clear all;
close all;
%给定城市坐标City_Coord，每一列代表一个城市
City_Coord = [0.4000, 0.2439, 0.1707, 0.2293, 0.5171, 0.8732, 0.6878, 0.8488, 0.6683, 0.6195;
              0.4439, 0.1463, 0.2293, 0.7610, 0.9414, 0.6536, 0.5219, 0.3609, 0.2536, 0.2634];
%获取城市数量city_quantity
city_quantity = size(City_Coord, 2);     %City_Coord的列数即为城市数量
%计算城市之间的距离City_Distance矩阵    %city_quantity行city_quantity列，第i行第j列代表城市i到城市j的距离
for i = 1:city_quantity
    for j = 1:city_quantity
        if i <= j
            City_Distance(i,j) = norm(City_Coord(:,i)-City_Coord(:,j));
            City_Distance(j,i) = City_Distance(i,j);
        end
    end
end

%----------设置参数----------%
T = 1000;         %初始温度
T_min = 1;     %最低温度
a = 0.99;         %降温系数
n = 0;           %迭代计数
K = 0.01;
%----------初始化----------%
X = X_init(city_quantity);                      %初始化
X_Function = GetFunction(city_quantity, City_Distance, X);    %计算路径长度
%----------主循环----------%
while T > T_min
    while 1
        %变异得到新解
        X1 = mutate(city_quantity, X);
        X1_Function = GetFunction(city_quantity, City_Distance, X1);

        %计算函数值差值
        DX = X1_Function - X_Function;

        if DX < 0
            break
        elseif rand < exp(-DX / (K * T))
            break
        end
    end
    
    %更新解
    X = X1;
    X_Function = X1_Function;
    
    %将当前路径存储进历史记录
    n = n + 1;
    History_X(n, :) = X;              %存储当前路径
    History_FX(n, :) = X_Function;    %存储当前路径的长度
    %降温
    T = T * a;
end
%**********************打印结果**********************%
fprintf('模拟退火算法：n')
fprintf('****************迭代结果****************n')
fprintf('最优路径为： %d %d %d %d %d %d %d %d %d %dn', X)
fprintf('最优路径的长度为： %fn', X_Function)
fprintf('****************************************n')
fprintf('endn')                                                                       
%**********************************************************************%
%****************************种群X初始化函数****************************%
%**********************************************************************%
function X = X_init(city_quantity)                   
        X = randperm(city_quantity);     %随机生成城市顺序
end
%**********************************************************************%
%*********************计算函数值（路径总长度）函数**********************%
%**********************************************************************%
function X_Function = GetFunction(city_quantity, City_Distance, X)
    X_Function = 0;
    for i = 1:(city_quantity - 1)
        X_Function = X_Function + City_Distance(X(i), X(i+1));
    end
    X_Function = X_Function + City_Distance(X(city_quantity), X(1));
end
%**********************************************************************%
%****************************变异函数****************************%
%**********************************************************************%
 function X = mutate(city_quantity, X)
    first = round(rand * (city_quantity - 1)) + 1;              %变异点1
    last =  round(rand * (city_quantity - first)) + first;      %变异点2
    temp = X(first);
    X(first) = X(last);
    X(last) = temp;
 end
Part 4 模拟退火算法解决函数极值问题
在极值问题中，粒子状态就是点的位置，内能就是函数值。
在上下界中随机产生一个数作为新解。
直接使用原有的Metropolis准则
使用 
 来缩减温度，其中a为一个小于1的数
以函数 
 求极值为例：
目标函数定义为：
function f = func(x)
    f = x^(1/2) + 7 * sin(5 * x) + 9 * cos(3 * x) + 2 * sin(x);
模拟退火算法：
clc;
clear all;
close all;
%----------设置参数----------%
T = 1000;         %初始温度
T_min = 1;     %最低温度
a = 0.99;         %降温系数
n = 0;           %迭代计数
K = 0.01;
lb = 0;
ub = 8;

X = lb + rand * (ub - lb);          %在上界、下界之间随机产生初始点

while T > T_min
    while 1
        %随机取新解
        X1 = lb + rand * (ub - lb);

        %计算函数值差值
        DX = func(X1) - func(X);

        if DX < 0
            break
        elseif rand < exp(-DX / (K * T))
            break
        end
    end
    
    %更新解
    X = X1;
    
    %将X存储进历史记录
    n = n + 1;
    History_X(n, :) = X;              %存储当前路径
    History_FX(n, :) = func(X);    %存储当前路径的长度
    %降温
    T = T * a;
end

%**********************打印结果**********************%
fprintf('模拟退火算法：n')
fprintf('****************迭代结果****************n')
fprintf('最小值点为： %fn', X)
fprintf('最小函数值为： %fn', func(X))
fprintf('****************************************n')
fprintf('endn')

模拟退火算法简介
下面以一个实例来讲解模拟退火算法求解TSP问题
第一步就是加载城市坐标，然后处理
第二步计算各个城市之间的距离，并创建距离矩阵，我采用的各个城市之间距离的计算公式为
距离矩阵是一个对称矩阵，对称线上为0
第三步就是随机生成一条线路，取得路线的路径相对较小的一条线路
最后一步就是模拟退火过程
降温系数设置为at=0.999，即new——T=at X old——T，终止温度选为e=10E-30，T设置为1。
在相对路径距离较小的基础上随机选择两个城市号进行交换，组成新的路径
交换i和j两个城市的位置，同时i和j中间的城市倒换顺序交换后的路径为
得到两条路径的距离差为
接受准则就是
即当df<0时，也就是新路径距离比老路径距离小，否则就是依概率接受新路径。
然后就是降温了，也就是判断是否满足迭代条件，如果满足就停止迭代了
最后的运行结果就是
最后附上完整代码

	

转载自：https://www.omegaxyz.com/2018/01/26/aco/
(一)蚁群算法的由来
蚁群算法(ant colony optimization)最早是由Marco Dorigo等人在1991年提出，他们在研究新型算法的过程中，发现蚁群在寻找食物时，通过分泌一种称为信息素的生物激素交流觅食信息从而能快速的找到目标，据此提出了基于信息正反馈原理的蚁群算法。
蚁群算法的基本思想来源于自然界蚂蚁觅食的最短路径原理，根据昆虫科学家的观察，发现自然界的蚂蚁虽然视觉不发达，但它们可以在没有任何提示的情况下找到从食物源到巢穴的最短路径，并在周围环境发生变化后，自适应地搜索新的最佳路径。
蚂蚁在寻找食物源的时候，能在其走过的路径上释放一种叫信息素的激素，使一定范围内的其他蚂蚁能够察觉到。当一些路径上通过的蚂蚁越来越多时，信息素也就越来越多，蚂蚁们选择这条路径的概率也就越高，结果导致这条路径上的信息素又增多，蚂蚁走这条路的概率又增加，生生不息。这种选择过程被称为蚂蚁的自催化行为。对于单个蚂蚁来说，它并没有要寻找最短路径，只是根据概率选择；对于整个蚁群系统来说，它们却达到了寻找到最优路径的客观上的效果。这就是群体智能。
(二)蚁群算法能做什么
蚁群算法根据模拟蚂蚁寻找食物的最短路径行为来设计的仿生算法，因此一般而言，蚁群算法用来解决最短路径问题，并真的在旅行商问题(TSP，一个寻找最短路径的问题)上取得了比较好的成效。目前，也已渐渐应用到其他领域中去，在图着色问题、车辆调度问题、集成电路设计、通讯网络、数据聚类分析等方面都有所应用。
(三)蚁群算法实现
优化的 函数为F(x,y)= -(x.^2+3*y.^4-0.2*cos(3*pi*x)-0.4*cos(4*pi*y)+0.6)
clear
clc
Ant= 300;%蚂蚁数量
Times= 80;%移动次数
Rou= 0.9;%荷尔蒙发挥系数
P0= 0.2;%转移概率
Lower_1= -1;%搜索范围
Upper_1= 1;
Lower_2= -1;
Upper_2= 1;
fori=1:Ant
    X(i,1)=(Lower_1+(Upper_1-Lower_1)*rand);
    X(i,2)=(Lower_1+(Upper_2-Lower_2)*rand);
    Tau(i)=F(X(i,1),X(i,2));
end
step=0.05;
f='-(x.^2+3*y.^4-0.2*cos(3*pi*x)-0.4*cos(4*pi*y)+0.6)';
[x,y]=meshgrid(Lower_1:step:Upper_1,Lower_2:step:Upper_2);
z=eval(f);
figure(1);
subplot(1,2,1);
mesh(x,y,z);
holdon;
plot3(X(:,1),X(:,2),Tau,'k*')
holdon;
text(0.1,0.8,-0.1,'蚂蚁的初始位置分布');
xlabel('x');ylabel('y');zlabel('f(x,y)');
forT=1:Times
    lamda=1/T;
    [Tau_Best(T),BestIndex]=max(Tau);
    fori=1:Ant
        P(T,i)=(Tau(BestIndex)-Tau(i))/Tau(BestIndex);%计算转移状态概率
    end
    fori=1:Ant
        ifP(T,i)<P0%局部搜索
            temp1=X(i,1)+(2*rand-1)*lamda;
            temp2=X(i,2)+(2*rand-1)*lamda;
        else%全局搜索
            temp1=X(i,1)+(Upper_1-Lower_1)*(rand-0.5);
            temp2=X(i,2)+(Upper_2-Lower_2)*(rand-0.5);
        end
        iftemp1<Lower_1%越界处理
            temp1=Lower_1;
        end
        iftemp1>Upper_1
            temp1=Upper_1;
        end
        iftemp2<Lower_2
            temp2=Lower_2;
        end
        iftemp2>Upper_2
            temp2=Upper_2;
        end
        ifF(temp1,temp2)>F(X(i,1),X(i,2))%更新位置
            X(i,1)=temp1;
            X(i,2)=temp2;
        end
    end
    fori=1:Ant
        Tau(i)=(1-Rou)*Tau(i)+F(X(i,1),X(i,2));%更新荷尔蒙
    end
end
subplot(1,2,2);
mesh(x,y,z);
holdon;
x=X(:,1);
y=X(:,2);
plot3(x,y,eval(f),'k*');
holdon;
text(0.1,0.8,-0.1,'蚂蚁的最终位置分布');
xlabel('x');ylabel('y');zlabel('f(x,y)');  
[max_value,max_index]=max(Tau);
maxX=X(max_index,1);
maxY=X(max_index,2);
maxValue=F(X(max_index,1),X(max_index,2));
优化函数：
functionf= F(x,y)
f= -(x.^2+3*y.^4-0.2*cos(3*pi*x)-0.4*cos(4*pi*y)+0.6);
end
效果：
作者：OmegaXYZ
https://www.omegaxyz.com/2018/01/26/aco/
	

以下为纯小白对卡尔曼滤波算法的MATLAB实现进行简单总结，以便以后复习。
一、 背景介绍
假设一小车作匀加速运动，初速度为0，加速度为5米每二次方秒，小车上装有速度传感器，采样频率为10Hz，传感器测量噪声为高斯白噪声，需要充分利用这些信息来估计车辆的速度状态，并验证卡尔曼滤波算法的实验原理与过程。
二、卡尔曼滤波原理
早在近百年前，就有人开始采用状态变量模型研究随机过程，随后为了解决对非平稳、多输入输出随机序列的估计问题，卡尔曼提出了递推最优估计理论。它采用状态空间法描述系统，由状态方程和测量方程组成，即知道前一个状态的估计值和最近一个观测数据，采用递推算法估计当前的状态值。其具有以下特点：

（1）	算法是一个递推过程，且状态空间法采用在时域内设计滤波器的方法，因而适用于多维随机过程的估计；离散型卡尔曼滤波算法适用于计算机处理。

（2）	采用递推算法计算，不用知道所有过去的值，用状态方程描述状态变量的动态变化规律，卡尔曼滤波同样适用于非平稳过程。

如图所示，为卡尔曼滤波算法的实现流程图，其基本思路是若有一组强而合理的假设，给出系统的历史测量值，可以建立最大化这些早前测量值的后验概率的系统状态模型。其基本假设是被建系统是线性的，影响测量的噪声属于符合高斯分布的白噪声。


三、卡尔曼滤波的状态方程和测量方程
整个实现过程包含预测更新和测量更新两大部分。

（1）	预测更新：

预测状态量：



预测误差协方差矩阵：



（2）	测量更新：

最优估计状态量：



计算误差增益：



误差协方差矩阵：



设k时刻小车速度为 ，则系统状态方程为：



测量方程为：



结合卡尔曼滤波算法的预测和测量更新流程，可有：


四、结果分析
根据以上理论分析，完成相应matlab程序（代码见后），并画出小车速度的观测值、真实值、滤波值的对比曲线图如下：


五、代码
(1) kalman.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %日 期： 2020.3.10
 %程序功能：使用卡尔曼滤波器估计小车匀加速运动时的速度状态
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
clc
clear all
A = [1];
B = [1];
U = 0.5;      %A,B,U初始值来自于状态方程
H = [1];       %H来自于测量方程
step = 500;
v = normrnd(0,10,1,step); %均值为0，方差为10噪声
X_0 = zeros(1,step);         %先验估计初始值
X0 = zeros(1,step);           %后验估计初始值
Z = zeros(1,step);             %观测值初始值
true = 0.5*(1:step);           %小车真实速度
temp = true + v(1,:);
Z(1,:) = temp;    %加入噪声的速度，模拟观测值
Q = [10];           %状态方程噪声协方差矩阵
R = [500];         %测量方程协方差矩阵
P0 = [1];           %初始后验估计协方差
P_0 = [1];
for i = 2:step
    [X_0(i),P_0] = kalmanPredictor(A,B,U,P0,Q,X0(i-1));
    [X0(i),P0] = kalmanUpdater(H,R,P_0,Z(i),X_0(i));
end
figure;
plot(Z,'g.');
hold on;
plot(true,'r.');
hold on;
plot(X0,'b');
grid on;
legend('观测值','真实值','kalman滤波');

(2)kalmanPredictor.m
function [Xp,Pp] = kalmanPredictor(A,B,U,P0,Q,X0)
Xp = A*X0 + B*U;
Pp = A*P0*A' + Q;
end

(3)kalmanUpdater.m
function [Xup,Pup] = kalmanUpdater(H,R,Pp,Z,Xp)
I = eye(size(Xp,1));
K = (Pp*H')/(H*Pp*H' + R);
Xup = Xp + K*(Z-H*Xp);
Pup = (I-K*H)*Pp;
end