带约束的优化算法

问题

O b j e c t F u n c t i o n : m i n X f ( X ) / m a x X f ( X ) s . t . { a ⩽ g ( X ) ⩽ b c ⩽ X ⩽ d . . . \begin{aligned}ObjectFunction:&\boldsymbol{\mathop{min}\limits_{X}f(X)}/\boldsymbol{\mathop{max}\limits_{X}f(X)} \\\boldsymbol{s.t.}&\begin{cases}\boldsymbol{a\leqslant\\g(X)\leqslant\\b} \\ \boldsymbol{c\leqslant\\X\leqslant\\d} \\ ... \end{cases} \end{aligned} ObjectFunctions.t.Xminf(X)/Xmaxf(X) ag(X)bcXd...
其中, X = [ x 1 , x 2 , x 3 , . . , x m ] \boldsymbol{X}=[x_{1},x_{2},x_{3},..,x_{m}] X=[x1,x2,x3,..,xm]

一、遗传算法(Genetic Algorithms,GA)

遗传算法是模拟达尔文生物进化论的自然选择和遗传学机理的生物进化过程的计算模型,是一种通过模拟自然进化过程搜索最优解的方法。遗传算法主要特点是直接对结构对象进行操作,不存在求导和函数连续性的限定;具有内在的隐并行性和更好的全局寻优能力;采用概率化的寻优方法,不需要确定的规则就能自动获取和指导优化的搜索空间,自适应地调整搜索方向。

种群

数学问题的多组可行解

个体

数学问题的某一组可行解

染色体

数学问题某一组可行解中的某个数据

基因

数学问题某一组可行解中的某个数据的二进制数串中的0和1

适应度函数

某些情况下,适应度函数可直接用目标函数来代替,但大多数情况下需要对目标函数做出变换才能作为适应度函数,比如目标函数之间差别很小时,个体被选择的概率也就相差不大,这样算法的选择功能将被大大弱化,这时就需要适应度函数设计和变换。
适应度函数需要根据业务需求来设计,对于所设计的适应度函数,其值相对偏差既不能过小也不能过大,差别过小会使得选择压力小,弱化了选优功能,差别过大会使得选择压力大,容易陷入局部最优,而且适应度函数值不能为负,否则概率无法计算。

适应度函数设计方法:
1、静态线性变换

F = a f + b F=af+b F=af+b,其中 f f f为目标函数, F F F为适应度函数
当目标是求 f f f的最大值时: a = 1 , b = − f m i n + ϵ a=1,b=-f_{min}+\epsilon a=1,b=fmin+ϵ
当目标是求 f f f的最小值时: a = − 1 , b = f m a x + ϵ a=-1,b=f_{max}+\epsilon a=1,b=fmax+ϵ
其中, ϵ \epsilon ϵ是一个较小的数,目的是使种群中最差的个体仍然有繁殖的机会,避免因为某个个体计算出的概率为0而导致永远没有繁殖的机会,增加种群的多样性。

2、动态线性变换

F = a k f + b k F=a^{k}f+b^{k} F=akf+bk,其中 f f f为目标函数, F F F为适应度函数
当目标是求 f f f的最大值时: a k = 1 , b k = − f m i n + ϵ k a^{k}=1,b^{k}=-f_{min}+\epsilon^{k} ak=1,bk=fmin+ϵk
当目标是求 f f f的最小值时: a k = − 1 , b k = f ( m a x ) + ϵ k a^{k}=-1,b^{k}=f(max)+\epsilon^{k} ak=1,bk=f(max)+ϵk
其中, ϵ 0 = m , ϵ k = ϵ k − 1 r \epsilon^{0}=m,\epsilon^{k}=\epsilon^{k-1}r ϵ0=m,ϵk=ϵk1r,可以通过调节 m m m r r r来调节 ϵ k \epsilon^{k} ϵk,加入 ϵ k \epsilon^{k} ϵk既可以使种群中最差的个体仍然有繁殖的机会,避免因为某个个体计算出的概率为0而导致永远没有繁殖的机会,也可以动态的调节选择压力(开始选择压力小,后来选择压力大),扩大搜索范围,增加种群的多样性。

3、幂函数变换

F = f α F=f^{\alpha} F=fα,其中 f f f为目标函数, F F F为适应度函数
α > 1 \alpha>1 α>1时选择压力增大, α < 1 \alpha<1 α<1时选择压力减小。

4、对数函数变换

F = a l g f + b F=algf+b F=algf+b,其中 f f f为目标函数, F F F为适应度函数
对数函数变换可以缩小目标函数值的差别。

5、指数函数变换

F = a e b f + c F=ae^{bf}+c F=aebf+c,其中 f f f为目标函数, F F F为适应度函数
指数函数变换既保证了适应度函数非负性,也可以使得目标函数值越大适应度函数值越小,目标函数值越小适应度函数越大。 α \alpha α取值越小,那些变换前具有较大适应度个体的新适应度与变换前较小适应度个体的新适应度值的差距减小;反之, α \alpha α取值越大,可将变换前适应度接近的个体的新适应度值差距加大。

6、正规化变换

F = f − f m i n + r f m a x − f m i n + r F=\frac{f-f_{min}+r}{f_{max}-f_{min}+r} F=fmaxfmin+rffmin+r,其中 f f f为目标函数, F F F为适应度函数
正规化变换可将 f f f的函数值映射到 ( 0 , 1 ) (0,1) (0,1)区间,抑制超级染色体。

交叉

从上一代的染色体中寻找两条染色体,将这两条染色体的某一个位置切断,各取一部分拼接在一起,从而生成一条新的染色体,其中交叉的两个个体采用轮盘赌法选择

变异

新染色体上随机选择若干个基因,然后随机修改基因的值,从而给现有的染色体引入新的基因

复制

上一代中适应度最高的几条染色体直接复制给下一代

应用举例

f ( x , y ) = x 2 s i n ( x + π ) − y c o s ( x y ) f(x,y)=x^{2}sin(x+\pi)-ycos(xy) f(x,y)=x2sin(x+π)ycos(xy)的最大值,其中 − 1 < x < 1 , − 1 < y < 1 -1<x<1,-1<y<1 1<x<1,1<y<1
流程:
1、初始化种群

2、计算初始种群中每个个体的适应度
每个个体的目标函数值如下:

由于目标函数值存在负数,因此不能简单的采用目标函数值作为适应度函数值,需进行一定的转化,本例采用正规化变换即 F = f − f m i n + r f m a x − f m i n + r F=\frac{f-f_{min}+r}{f_{max}-f_{min}+r} F=fmaxfmin+rffmin+r来作为适应度函数,其中 r r r取0.05。

3、利用轮盘赌选择个体进行交叉
p i = F ( x i , y i ) ∑ i = 1 n F ( x i , y i ) p_{i}=\frac{F(x_{i},y_{i})}{\sum_{i=1}^{n}F(x_{i},y_{i})} pi=i=1nF(xi,yi)F(xi,yi)

设置交叉概率,交叉概率一般较大以保证优胜劣汰,避免寻优过程缓慢,本例中设为0.7。
假设选中的交配个体为 [ − 0.6848 , 0.6250 ] [-0.6848,0.6250] [0.6848,0.6250] [ − 0.4430 , 0.6006 ] [-0.4430,0.6006] [0.4430,0.6006]

4、变异
设置变异概率,变异概率一般较小以避免破坏了优胜劣汰而成为随机搜索,本例中设为0.01。
假设选中的变异个体为KaTeX parse error: Expected 'EOF', got '#' at position 92: …3a932519ba7.bmp#̲pic_center" wid…[-0.8049,-0.0292] 和 和 [0.0938,-0.7162]$适应度更强,可直接选择复制成为下一代。
6、判断
判断是否满足终止条件,满足条件则终止迭代,反之则更新种群,再转至第二步,继续进行迭代寻优。

clear;
% 问题: f(x,y) = x * x * sin(x+pi) - y * cos(x*y) 其中,-5 =< x,y <= 5
% GA算法
lowBound = [-5,-5];
upperBound = [5,5];
maxIter = 1000;
iter = zeros(maxIter,1);
maxObjFuncValue = zeros(maxIter,1);
minObjFuncValue = zeros(maxIter,1);
meanFitness = zeros(maxIter,1);
maxFitness = zeros(maxIter,1);
CrossProbability = 0.75;                  % 交叉概率
MutationProbability = 0.01;               % 变异概率  
ChromosomeNum = 2;                        % 染色体数目
PopulationNum = 150;                      % 种群数目
KnockOutNum = floor(PopulationNum * 0.05);% 淘汰数目
% --------------------种群进化-------------------- %
k = 1;
iter(k) = k;
PopulationMatrix = rand(PopulationNum,ChromosomeNum) .* (upperBound - lowBound) + lowBound;
ObjFuncValue = zeros(PopulationNum,1);
Fitness = zeros(PopulationNum,1);
% ----------计算目标函数值---------- %
for i = 1:PopulationNum
    x = PopulationMatrix(i,1);
    y = PopulationMatrix(i,2);
    ObjFuncValue(i) = x * x * sin(x+pi) - y * cos(x*y);
end
% ----------计算适应度---------- %
maxValue = max(ObjFuncValue);
minValue = min(ObjFuncValue);
maxObjFuncValue(k) = maxValue;
minObjFuncValue(k) = minValue;
r = (maxValue - minValue) * 0.02;
for i = 1:PopulationNum
    Fitness(i) = (ObjFuncValue(i) - minValue + r) / (maxValue - minValue + r);
end
meanFitness(k) = sum(Fitness) / PopulationNum;
% 更新最优个体信息以及最劣个体信息
[minFitness(k),bestIndex] = min(Fitness);
[maxFitness(k),worstIndex] = max(Fitness);
bestIndividuality = PopulationMatrix(bestIndex,:);
worstIndividuality = PopulationMatrix(worstIndex,:);
% ----------优胜劣汰---------- %
while(k <= maxIter)
    iter(k) = k;
    % ----------计算概率---------- %
    p = Fitness / sum(Fitness);
    CumProbability = cumsum(p);
    % ----------交叉---------- %
    if (rand < CrossProbability)
        fatherIndex = min(find(CumProbability > rand,PopulationNum,'first'));
        father = PopulationMatrix(fatherIndex,:);
        motherIndex = min(find(CumProbability > rand,PopulationNum,'first'));
        mother = PopulationMatrix(motherIndex,:);
        child = mother + rand * (father - mother);
        x = child(1,1);
        y = child(1,2);
        ObjFuncValue_cross = x * x * sin(x+pi) - y * cos(x*y);
        Fitness_cross = (ObjFuncValue_cross - minValue + r) / (maxValue - minValue + r);
        if (Fitness_cross > min(Fitness))
             % 用child取代最劣个体
             PopulationMatrix(worstIndex,:) = child;
             ObjFuncValue(worstIndex) = ObjFuncValue_cross;
        end
    end
    % ----------变异---------- %
    NewIndividual = zeros(1,2);
    if (rand < MutationProbability)
        MutatedIndex = floor(PopulationNum * rand + 1);
        MutatedIndividual = PopulationMatrix(MutatedIndex,:);
        for i = 1:ChromosomeNum
            if (rand > 0.5)
                % 增加
                NewIndividual(i) = MutatedIndividual(i) + rand * (upperBound(i) - MutatedIndividual(i));
            else
                % 减小
                NewIndividual(i) = MutatedIndividual(i) + rand * (lowBound(i) - MutatedIndividual(i));
            end
        end
        % 计算目标函数值
        x = NewIndividual(1,1);
        y = NewIndividual(1,2);
        ObjFuncValue_Mutated = x * x * sin(x+pi) - y * cos(x*y);
        % 更新种群以及目标函数值
        PopulationMatrix(MutatedIndex,:) = NewIndividual;
        ObjFuncValue(MutatedIndex) = ObjFuncValue_Mutated;
    end
    % ----------复制---------- %
    % 冒泡排序
    FitnessSorted = Fitness;
    IndexSorted = (1:PopulationNum)';
    for m = 1 : PopulationNum
        for n = 1 : PopulationNum - m
            value =  FitnessSorted(n);
            if (value > FitnessSorted(n+1))
                FitnessSorted(n) = FitnessSorted(n+1);
                FitnessSorted(n+1) = value;
                index = IndexSorted(n);
                IndexSorted(n) = IndexSorted(n+1);
                IndexSorted(n+1) = index;
            end
        end
    end
    % 前KnockOutNum个优质个体替换KnockOutNum个劣质个体
    for i = 1 : KnockOutNum
        PopulationMatrix(IndexSorted(i),:) = PopulationMatrix(IndexSorted(PopulationNum-i+1),:);
        ObjFuncValue(IndexSorted(i)) = ObjFuncValue(IndexSorted(PopulationNum-i+1));
    end
    % ----------计算适应度---------- %
    maxValue = max(ObjFuncValue);
    minValue = min(ObjFuncValue);
    maxObjFuncValue(k) = maxValue;
    minObjFuncValue(k) = minValue;
    r = 0.5;
    for i = 1:PopulationNum
        Fitness(i) = (ObjFuncValue(i) - minValue + r) / (maxValue - minValue + r);
    end
    meanFitness(k) = sum(Fitness) / PopulationNum;
    % 更新最优个体信息以及最劣个体信息
    [minFitness(k),bestIndex] = min(Fitness);
    [maxFitness(k),worstIndex] = max(Fitness);
    bestIndividuality = PopulationMatrix(bestIndex,:);
    worstIndividuality = PopulationMatrix(worstIndex,:);
    % 迭代器自增
    k = k + 1;
end
figure(1);
plot(iter,maxFitness,'-','Linewidth',2.5);
hold on;
plot(iter,meanFitness,'-','Linewidth',2.5);
figure(2);
plot(iter,maxObjFuncValue,'-','Linewidth',2.5);

二、粒子群算法(Particle Swarm Optimization,PSO)

在粒子群算法中,将被优化问题的每个解当作一只鸟,称之为粒子。所有的粒子都有一个由被优化的函数决定的适应度值(解的好坏程度),每个粒子还有一个速度(通俗点说就是调整位置的系数)决定它们飞翔的方向。然后,粒子们就追随当前的最优粒子在解空间中搜索。粒子群算法首先在给定的解空间中随机初始化粒子群,每个粒子有了初始位置与初始速度,然后通过迭代寻优。在每一次迭代中,每个粒子通过跟踪个体极值,和全局极值来更新自己在解空间中的空间位置与飞行速度,最终找到全局最优。

参数

i i i个粒子在 d d d维空间中的位置

X i d = ( x i 1 , x i 2 , x i 3 , . . . , x i d ) X_{id}=(x_{i1},x_{i2},x_{i3},...,x_{id}) Xid=(xi1,xi2,xi3,...,xid)

i i i个粒子在 d d d维空间中的速度

V i d = ( v i 1 , v i 2 , v i 3 , . . . , v i d ) V_{id}=(v_{i1},v_{i2},v_{i3},...,v_{id}) Vid=(vi1,vi2,vi3,...,vid)

i i i个粒子在 d d d维空间搜索到的最优位置

P i d , p B e s t = ( p i 1 , p B e s t , p i 2 , p B e s t , p i 3 , p B e s t , . . . , p i d , p B e s t ) P_{id,pBest}=(p_{i1,pBest},p_{i2,pBest},p_{i3,pBest},...,p_{id,pBest}) Pid,pBest=(pi1,pBest,pi2,pBest,pi3,pBest,...,pid,pBest)

粒子群体在 d d d维空间搜索到的最优位置

P d , g B e s t = ( p 1 , g B e s t , p 2 , g B e s t , p i , g B e s t , . . . , p i , g B e s t ) P_{d,gBest}=(p_{1,gBest},p_{2,gBest},p_{i,gBest},...,p_{i,gBest}) Pd,gBest=(p1,gBest,p2,gBest,pi,gBest,...,pi,gBest)

i i i个粒子在 d d d维空间搜索到的最优位置的适应值

f p f_{p} fp-个体历史最优适应值

粒子群体在 d d d维空间搜索到的最优位置的适应值

f g f_{g} fg-群体历史最优适应值

粒子速度更新公式

v i d k + 1 = ω v i d k + c 1 r 1 ( p i d , p B e s t k − x i d k ) + c 2 r 2 ( p d , g B e s t k − x i d k ) v_{id}^{k+1}=\omega v_{id}^{k}+c_{1}r_{1}(p_{id,pBest}^{k}-x_{id}^{k})+c_{2}r_{2}(p_{d,gBest}^{k}-x_{id}^{k}) vidk+1=ωvidk+c1r1(pid,pBestkxidk)+c2r2(pd,gBestkxidk)

第一项:惯性部分

由惯性权重和粒子自身速度构成,表示粒子对先前自身运动状态的信任。

第二项:认知部分

表示粒子本身的思考,即粒子自己经验的部分,可理解为粒子当前位置与自身历史最优位置之间的距离和方向。

第三项:社会部分

表示粒子之间的信息共享与合作,即来源于群体中其他优秀粒子的经验,可理解为粒子当前位置与群体历史最优位置之间的距离和方向。
其中, N N N为粒子群规模, i i i为粒子序号, i = 1 , 2 , 3 , . . . , N i=1,2,3,...,N i=1,2,3,...,N D D D为粒子维度, d d d为粒子维度序号, d = 1 , 2 , 3 , . . . , D d=1,2,3,...,D d=1,2,3,...,D k k k为迭代次数; ω \omega ω为惯性权重; c 1 c_{1} c1为个体学习因子; c 2 c_{2} c2为群体学习因子; r 1 、 r 2 r_{1}、r_{2} r1r2为区间 [ 0 , 1 ] [0,1] [01]内的随机数,增加搜索的随机性; x i d k x_{id}^{k} xidk为第 i i i个粒子在第 k k k次迭代中第 d d d维度的位置; v i d k v_{id}^{k} vidk为第 i i i个粒子在第 k k k次迭代中第 d d d维度的速度; p i d , p B e s t k p_{id,pBest}^{k} pidpBestk为第 i i i个粒子在第 k k k次迭代中第 d d d维度的历史最优位置,即在第 k k k次迭代后第 i i i个粒子搜索得到的最优解; p d , g B e s t k p_{d,gBest}^{k} pdgBestk为群体在第 k k k次迭代中第 d d d维度的历史最优位置,即在第 k k k次迭代后整个粒子群的最优解;

粒子位置更新公式

x i d k + 1 = x i d k + v i d k + 1 x_{id}^{k+1}=x_{id}^{k}+v_{id}^{k+1} xidk+1=xidk+vidk+1

参数的选择

粒子群规模 N N N

推荐范围: [ 20 , 1000 ] [20,1000] [20,1000]。简单问题一般取2040,较难或特定类别的问题可以取100200。较小的种群规模容易陷入局部最优;较大的种群规模可以提高收敛性,更快找到全局最优解,但是相应地每次迭代的计算量也会增大;当种群规模增大至一定水平时,再增大将不再有显著的作用。

粒子维度 D D D

粒子搜索的空间维数即为自变量的个数。

迭代次数 k k k

推荐范围: [ 50 , 100 ] [50,100] [50,100]。典型取值:60、70、100;这需要在优化的过程中根据实际情况进行调整,迭代次数太小的话解不稳定,太大的话非常耗时,没有必要。

惯性权重 ω \omega ω

在解决实际优化问题时,往往希望先采用全局搜索,使搜索空间快速收敛于某一区域,然后采用局部精细搜索以获得高精度的解。因此提出了自适应调整的策略,即随着迭代的进行,线性地减小 ω \omega ω的值。这里提供一个简单常用的方法——线性变化策略:随着迭代次数的增加,惯性权重 ω \omega ω不断减小,从而使得粒子群算法在初期具有较强的全局收敛能力,在后期具有较强的局部收敛能力。
ω = i t e r i t e r m a x ω m i n + ( 1 − i t e r i t e r m a x ) ω m a x \omega=\frac{iter}{iter_{max}}\omega_{min}+(1-\frac{iter}{iter_{max}})\omega_{max} ω=itermaxiterωmin+(1itermaxiter)ωmax
其中, ω m i n 、 ω m a x \omega_{min}、\omega_{max} ωminωmax分别为最小惯性权重和最大惯性权重; i t e r m i n 、 i t e r m a x iter_{min}、iter_{max} iterminitermax分别为当前迭代次数和最大迭代次数;

学习因子 c 1 、 c 2 c_{1}、c_{2} c1c2

学习因子也称为加速系数或加速因子,推荐范围: [ 0 , 4 ] [0,4] [0,4],典型取值有 c 1 = c 2 = 2 c_{1}=c_{2}=2 c1=c2=2 c 1 = 1.6 、 c 2 = 1.8 c_{1}=1.6、c_{2}=1.8 c1=1.6c2=1.8以及 c 1 = 1.6 、 c 2 = 2 c_{1}=1.6、c_{2}=2 c1=1.6c2=2
c 1 c_{1} c1表示粒子下一步动作来源于自身经验部分所占的权重,将粒子推向个体最优位置 P i d , p B e s t P_{id,pBest} Pid,pBest的加速权重;
c 2 c_{2} c2表示粒子下一步动作来源于其它粒子经验部分所占的权重,将粒子推向群体最优位置 P d , g B e s t P_{d,gBest} Pd,gBest的加速权重;

速度限制 v m a x v_{max} vmax

为了平衡算法的探索能力与开发能力,需要设定一个合理的速度范围,限制粒子的最大速度 v m a x v_{max} vmax,即粒子下一步迭代可以移动的最大距离。 v m a x v_{max} vmax较大时,粒子飞行速度快,探索能力强,但粒子容易飞过最优解; v m a x v_{max} vmax较小时,粒子飞行速度快,探索能力强,但粒子容易飞过最优解;因此, v m a x v_{max} vmax一般设为粒子变化范围的10%~20%,可根据实际情况调试,但不能大于粒子(解)的变化范围。(北理工韩宝玲教授的"基于粒子群算法的四足机器人静步态优化方法[4]"论文中采用了拉丁方抽样方法来解决粒子初始化问题,可以尝试采用这种方法初始化种群)

初始种群

粒子群算法优化的结果受很多因素的影响,其中受粒子初始值的影响比较大,而且较难调控。如果粒子初始值是随机初始化的,在不改变任何参数的情况下,多次优化的结果不一定都收敛到一个全局或局部最优解,也可能会得到一个无效解。所以粒子初始化是一个十分重要的步骤,它关系到整个优化过程中优化收敛的速度与方向。如果粒子的初始化范围选择得好的话可以大大缩短优化的收敛时间,也不易于陷入局部最优解。我们需要根据具体的问题进行分析,如果根据我们的经验判断出最优解一定在某个范围内,则就在这个范围内初始化粒子。如果无法确定,则以粒子的取值边界作为初始化范围。

应用举例

f ( x , y ) = x 2 s i n ( x + π ) − y c o s ( x y ) f(x,y)=x^{2}sin(x+\pi)-ycos(xy) f(x,y)=x2sin(x+π)ycos(xy)的最大值,其中 − 1 < x < 1 , − 1 < y < 1 -1<x<1,-1<y<1 1<x<1,1<y<1
流程:
1、初始化种群(群体规模、每个粒子的位置以及速度)
2、计算每个粒子的适应度
3、记录每个粒子的最优位置
4、寻找种群的最优位置
5、迭代更新每个粒子的位置和速度,并对边界条件进行处理
6、判断极值是否满足条件,若满足条件则结束迭代;反之则返回步骤二继续迭代寻优

clear;
% 问题: f(x,y) = x * x * sin(x+pi) - y * cos(x*y) 其中,-5 =< x,y <= 5
% PSO算法
maxIter = 1000;
iter = zeros(maxIter,1);
ParticleNum = 150;                              % 粒子群群数目
ParticleDim = 2;                                % 粒子群维度
c1 = 1.5;                                       % 学习因子1
c2 = 1.5;                                       % 学习因子2
w = 0.75;                                       % 惯性权重
minPos = [-5,-5];                               % 最小位置
maxPos = [5,5];                                 % 最大位置
minVel = [-1,-1];                               % 最小速度
maxVel = [1,1];                                 % 最大速度
gOptimalObjFuncValue = zeros(maxIter,1);        % 全局最优值
gOptimalPos = zeros(maxIter,ParticleDim);       % 全局最优位置
% --------------------种群进化-------------------- %
k = 1;
iter(k) = k;
ParticlePosition = rand(ParticleNum,ParticleDim) .* (maxPos - minPos) + minPos;
ParticleVelocity = rand(ParticleNum,ParticleDim) .* (maxVel - minVel) + minVel;
ObjFuncValue = zeros(ParticleNum,1);
% ----------计算目标函数值---------- %
for i = 1:ParticleNum
    x = ParticlePosition(i,1);
    y = ParticlePosition(i,2);
    ObjFuncValue(i) = x * x * sin(x+pi) - y * cos(x*y);
end
% ----------初始化最优值和最优位置---------- %
pOptimalPos = ParticlePosition;% 每个个体最优位置
pOptimalObjFuncValue = ObjFuncValue;% 每个个体最优函数值
[gOptimalObjFuncValue(k),Index] = min(pOptimalObjFuncValue);
gOptimalPos(k,:) = pOptimalPos(Index,:);% 全局最优位置
% ----------迭代寻优---------- %
while(k <= maxIter)
    % -----更新位置和速度----- %
    for i = 1:ParticleNum
        ParticleVelocity(i,:) = w * ParticleVelocity(i,:)...
            + c1 * rand * (pOptimalPos(i,:) - ParticlePosition(i,:))...
            + c2 * rand * (gOptimalPos(k,:) - ParticlePosition(i,:));
        ParticlePosition(i,:) = ParticlePosition(i,:) + ParticleVelocity(i,:);
    end
    % -----边界条件处理----- %
    for i = 1:ParticleNum
        for j = 1:ParticleDim
            if (ParticlePosition(i,j) < minPos(1,j) || ParticlePosition(i,j) > maxPos(1,j))
                ParticlePosition(i,j) = rand * (maxPos(1,j) - minPos(1,j)) + minPos(1,j);
            end
            if (ParticleVelocity(i,j) < minVel(1,j) || ParticleVelocity(i,j) > maxVel(1,j))
                ParticleVelocity(i,j) = rand * (maxVel(1,j) - minVel(1,j)) + minVel(1,j);
            end
        end
    end
    % 迭代器自增
    k = k + 1;
    iter(k) = k;
    % ----------更新最优值和最优位置---------- %
    for i = 1:ParticleNum
        x = ParticlePosition(i,1);
        y = ParticlePosition(i,2);
        ObjFuncValue(i) = x * x * sin(x+pi) - y * cos(x*y);
        if (pOptimalObjFuncValue(i) > ObjFuncValue(i))
            pOptimalPos(i,:) = ParticlePosition(i,:);
            pOptimalObjFuncValue(i) = ObjFuncValue(i);
        end
    end
    [gOptimalObjFuncValue(k),Index] = min(pOptimalObjFuncValue);
    gOptimalPos(k,:) = pOptimalPos(Index,:);% 全局最优位置
end
figure(1);
plot(iter,gOptimalObjFuncValue,'-','Linewidth',2.5);
figure(2);
plot(gOptimalPos(:,1),gOptimalPos(:,2),'o','Linewidth',2.5);

三、拉格朗日乘子乘子(Lagrange Multiplier,LM)

拉格朗日乘子法是一种寻找多元函数在一组约束下的极值的方法。通过引入拉格朗日乘子,可将有 m m m个变量与 n n n个约束条件的最优化问题转化为具有 m + n m+n m+n个变量的无约束优化问题求解。

情况一:最优解 X ∗ \boldsymbol{X^{*}} X在不等式约束边界上时,在最优解 X ∗ \boldsymbol{X^{*}} X处,梯度 ▽ f ( X ) \boldsymbol{\bigtriangledown f(X)} f(X)与梯度 ▽ g ( X ) \boldsymbol{\bigtriangledown g(X)} g(X)共线

对于约束 a ⩽ g ( X ) ⩽ b \boldsymbol{a\leqslant\\g(X)\leqslant\\b} ag(X)b的边界上任意一点 X \boldsymbol{X} X,该点 X \boldsymbol{X} X的梯度 ▽ g ( X ) \boldsymbol{\bigtriangledown g(X)} g(X)正交与约束 g ( X ) \boldsymbol{g(X)} g(X);同理,对于目标函数 f ( X ) 上 \boldsymbol{f(X)}上 f(X)任意一点 X \boldsymbol{X} X,该点处的梯度 ▽ f ( X ) \boldsymbol{\bigtriangledown f(X)} f(X)正交与目标函数 f ( X ) \boldsymbol{f(X)} f(X);因此,在最优解 X ∗ \boldsymbol{X^{*}} X处,梯度 ▽ f ( X ) \boldsymbol{\bigtriangledown f(X)} f(X)与梯度 ▽ g ( X ) \boldsymbol{\bigtriangledown g(X)} g(X)共线,即存在 λ ≠ 0 \lambda≠0 λ=0使得 ▽ f ( X ) + λ ▽ g ( X ) = 0 \boldsymbol{\bigtriangledown f(X)}+\lambda\boldsymbol{\bigtriangledown g(X)}=0 f(X)+λg(X)=0

情况二:最优解 X ∗ \boldsymbol{X^{*}} X在不等式约束内部时,约束不起作用

当约束不起作用时,可直接利用 ▽ f ( X ) = 0 \boldsymbol{\bigtriangledown f(X)}=0 f(X)=0来求解,即 ▽ f ( X ) + λ ▽ g ( X ) = 0 \boldsymbol{\bigtriangledown f(X)}+\lambda\boldsymbol{\bigtriangledown g(X)}=0 f(X)+λg(X)=0,其中 λ = 0 \lambda=0 λ=0
总结:对于上述不等约束下的优化问题可以转化为新约束下的拉格朗日函数优化问题,即: L ( X , λ ) = f ( X ) + λ g ( X ) \boldsymbol{L(X},\lambda\boldsymbol{)}=\boldsymbol{f(X)}+\lambda\boldsymbol{g(X)} L(X,λ)=f(X)+λg(X) a ⩽ g ( X ) ⩽ b \boldsymbol{a\leqslant\\g(X)\leqslant\\b} ag(X)b以及 λ g ( X ) = 0 \lambda\boldsymbol{g(X)}=0 λg(X)=0下的优化问题

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值