双层优化入门(3)—基于智能优化算法的求解方法(附matlab代码)

简介: 除了数学规划方法之外,还可采用智能优化算法求解双层优化问题,一般在上层优化中采用智能优化算法,下层优化使用数学规划方法;也可以在上下层优化中都采用智能优化算法,这篇博客将进行详细介绍。算例依旧使用上面两篇博客中的线性双层优化问题,由于这个优化问题比较简单,我们采用最基础的粒子群算法进行求解。​。

         前面两篇博客介绍了双层优化的基本原理和使用KKT条件求解双层优化的方法,以及使用yalmip工具箱求解双层优化的方法。

       除了数学规划方法之外,还可采用智能优化算法求解双层优化问题,一般在上层优化中采用智能优化算法,下层优化使用数学规划方法;也可以在上下层优化中都采用智能优化算法,这篇博客将进行详细介绍。

       算例依旧使用上面两篇博客中的线性双层优化问题,由于这个优化问题比较简单,我们采用最基础的粒子群算法进行求解。

image.gif

1.粒子群算法

       1995年,受到鸟群觅食行为的规律性启发,James Kennedy和Russell Eberhar建立了一个简化算法模型,经过多年改进最终形成了粒子群优化算法(Particle Swarm Optimization, PSO) ,也可称为粒子群算法。

       关于算法的原理及步骤可以参考这篇博客:

粒子群优化算法(Particle Swarm Optimization, PSO)的详细解读 - 知乎 (zhihu.com)

       我们假设目标函数是平方和最小,也就是:

image.gif

       假设种群数为200,问题维度为10,最大迭代次数为200,位置上下限分别为10,-10,速度上下限为5,-5,惯性权重取固定值0.9,自我学习因子以及群体学习因子都取1.5,Matlab代码如下:

%% 清除变量
clc
clear
close all
warning off
%% 设置种群参数
sizepop = 200;                      % 初始种群个数
dim = 10;                           % 空间维数
ger = 200;                          % 最大迭代次数    
x_max = 10*ones(1,dim);             % 位置上限
x_min = -10*ones(1,dim);            % 位置下限
v_max = 5*ones(1,dim);              % 速度上限
v_min = -5*ones(1,dim);             % 速度下限
w=0.9;                              % 惯性权重
c_1 = 1.5;                          % 自我学习因子
c_2 = 1.5;                          % 群体学习因子 
%% 种群初始化
pop=x_min+rand(sizepop,dim).*(x_max-x_min);     % 初始化种群
pop_v=v_min+rand(sizepop,dim).*(v_max-v_min);   % 初始化种群速度        
pop_zbest=pop(1,:);                             % 初始化群体最优位置
pop_gbest=pop;                                  % 初始化个体最优位置
fitness=zeros(1,sizepop);                       % 所有个体的适应度
fitness_zbest=inf;                              % 初始化群体最优适应度
fitness_gbest=inf*ones(1,sizepop);              % 初始化个体最优适应度
% 初始的适应度
for k=1:sizepop
    % 计算适应度值
    fitness(k)=sum(pop(k,:).^2);
    if fitness(k)<fitness_zbest
        fitness_zbest=fitness(k);
        pop_zbest=pop(k,:);
    end
end
history_pso=zeros(1,ger);            % 粒子群历史最优适应度值
%% 迭代求最优解
iter=1;
while iter <= ger
    for k=1:sizepop
        % 更新速度并对速度进行边界处理 
        pop_v(k,:)= w * pop_v(k,:) + c_1*rand*(pop_gbest(k,:)-pop(k,:))+c_2*rand*(pop_zbest-pop(k,:));
        for kk=1:dim
            if  pop_v(k,kk) > v_max(kk)
                pop_v(k,kk) = v_max(kk);
            end
            if  pop_v(k,kk) < v_min(kk)
                pop_v(k,kk) = v_min(kk);
            end
        end
        % 更新位置并对位置进行边界处理
        pop(k,:)=pop(k,:)+pop_v(k,:);
        for kk=1:dim
            if  pop(k,kk) > x_max(kk)
                pop(k,kk) = x_max(kk);
            end
            if  pop(k,kk) < x_min(kk)
                pop(k,kk) = x_min(kk);
            end
        end
        % 更新适应度值
        fitness(k)=sum(pop(k,:).^2);
        if fitness(k)<fitness_zbest
            fitness_zbest=fitness(k);
            pop_zbest=pop(k,:);
        end
        if fitness(k)<fitness_gbest(k)
            fitness_gbest(k)=fitness(k);
            pop_gbest(k,:)=pop(k,:);
        end
    end
    history_pso(iter)=fitness_zbest;
    disp(['PSO第',num2str(iter),'次迭代最优适应度=',num2str(fitness_zbest)])
    iter=iter+1;
end
disp(['最优解:x=',num2str(pop_zbest)])
disp(['最优函数值=',num2str(fitness_zbest)])
plot(history_pso,'linewidth',1)
ylabel('最优适应度值')
xlabel('迭代次数')

image.gif

       运行结果如下:

image.gif

image.gif

       我们知道最优解是当所有变量都取0时,最优函数值为0,上面显示采用粒子群算法求出的显然并不是全局最优解,只是一个局部最优解。这也是智能优化算法无法避免的问题,即使是一个非常简单的目标函数,求出的结果也无法保证是全局最优,那么当目标函数变复杂时,情况将会更糟糕。现在对智能优化算法的研究非常多,各种动植物园算法、各种改进都层出不穷,但还是无法从根本上解决算法无法保证全局收敛的问题。

       所以,只有在数学模型比较复杂,非线性条件很多,而且对结果的误差是可以接受的情况下,才建议使用智能优化算法进行求解。

       下面以粒子群算法为例,介绍采用智能优化算法求解双层优化问题的方法。

2.智能优化算法中对约束条件的处理

       上面提供的代码中,含有的约束条件是x∈[-10,10],直接可以通过粒子位置的上下限进行约束,但其他形式的约束就没办法这样处理了,常用的处理方式有三类,以上面提到的双层优化问题的上层优化为例,分别介绍这几种方法:

image.gif

1)通过加罚函数的方式将有约束的问题转为无约束问题:

       首先定义罚函数:

image.gif

        然后将罚函数加入原来的目标函数中形成新的目标函数:

image.gif

其中,p是惩罚因子。这样处理,可以使得违反约束的粒子目标函数值变得很大,从而促使粒子朝着遵守约束的区域聚集。这种方式在过去非常常用,但我自己不太喜欢这种方式,原因待会我会解释。假设惩罚因子为1,新的目标函数可以用matlab代码表示出来:

function fitness=fitness_fun_method1(pop)
    x=pop(1);
    y=pop(2);
    punish=0;
    p=1;
    if -2*x+3*y > 12
        punish=punish+p*(-2*x+3*y-12);
    end
    if x+y > 14
        punish=punish+p*(x+y-14);
    end
    fitness=-x-2*y+punish;
end

image.gif

       运行结果如下:

image.gif

       很奇怪的是,明明最优解是x=6,y=8,最优函数值是-22,但求出的结果却完全不一样,甚至是违反约束条件的。

       这就是罚函数的缺点,因为你把约束条件当作罚函数加进目标函数里,当惩罚因子选的不合适时,实际上是改变了目标,很有可能导致最优解和最优函数值都发生了变化,和你原来的目标背道而驰了。解决的办法就是要选择一个合适的惩罚因此,比如我改成9,再重新求解:

image.gif

       这时候得到的结果就是正常的了。

       但说实话,我在运行之前也不知道我滴惩罚因子到底多大合适,如果要求解的问题比较复杂,运行时间很长,可能一个错误的惩罚因子就能让你白忙活半天。所以我不太喜欢这种方式。

2)当约束被违反时,给这个粒子赋一个很差的适应度值;

       第二种方式说起来就很简单,就是让违反约束的粒子适应度很差。如果是最小化问题就给他赋一个很大的值,如果是最大化问题就赋一个很小的值。这时候目标函数Matlab代码可以这样写:

function fitness=fitness_fun_method2(pop)
    x=pop(1);
    y=pop(2);
    if -2*x+3*y > 12
        fitness=999;
    end
    if x+y > 14
        fitness=999;
    end
    if -2*x+3*y <= 12 && x+y <= 14
        fitness=-x-2*y;
    end
end

image.gif

运行时可以求出最优解:

image.gif

3)初始化粒子时,仅产生满足约束条件的粒子,当约束被违反时,重新生成一个满足约束的粒子。

       这种方式和其他两种方式的角度不同,是从粒子的初始化和更新的角度进行考虑,而不是从适应度函数的角度出发。

       首先,在初始化种群时,保证生成的粒子都是满足约束条件的,可以采用下面的代码:

while true
    pop0=x_min+rand(sizepop*10,dim).*(x_max-x_min);
    x = pop0(:,1);
    y = pop0(:,2);
    % 使用逐元素运算检查约束条件
    mask = (-2*x+3*y<=12) & (x+y<=14) ;
    % 如果点的数量已经足够了,则退出循环
    if nnz(mask) >= sizepop
        pop=[x(mask), y(mask)];
        break;
    end
end
pop=pop(1:sizepop,:);                           % 初始化种群

image.gif

       目标函数可以不考虑任何约束条件:

function fitness=fitness_fun_method3(pop)
    x=pop(1);
    y=pop(2);
    fitness=-x-2*y;
end

image.gif

       但是当更新粒子之后如果不满足约束条件,需要将该粒子修复为满足约束的粒子:

% 更新位置并对位置进行边界处理
pop(k,:)=pop(k,:)+pop_v(k,:);
x=pop(k,1);
y=pop(k,2);
% 修复不满足约束的粒子
if (-2*x+3*y > 12 )|| (x+y > 14)
    while true
        pop0=x_min+rand(10,dim).*(x_max-x_min);
        x = pop0(:,1);
        y = pop0(:,2);
        % 使用逐元素运算检查约束条件
        mask = (-2*x+3*y<=12) & (x+y<=14) ;
        % 如果数量已经足够了,则退出循环
        if nnz(mask) >= 1
            pop_test=[x(mask), y(mask)];
            pop(k,:)=pop_test(1,:);
            break;
        end
    end
end

image.gif

       优化结果如下:

image.gif

       我们可以看到,在初始化种群以及更新粒子时都考虑约束条件时,也可以求出较优的解。

       此外根据所求优化问题的实际情况不同,也可以将上面几种对约束条件的不同处理方式相结合。例如,在初始化以及更新种群时考虑一部分约束,另一部分约束采用罚函数的方式添加到目标函数中,此处不再赘述。

3.上层优化采用粒子群算法,下层优化使用yalmip或粒子群算法进行求解

       对于双层优化的问题求解,第一种方式是上层优化采用智能优化算法,下层问题使用yalmip等工具箱直接求出最优解。第二种方式是上下层优化都采用智能优化算法。还是以上面那个线性双层优化问题为例,采用第二种方法处理约束条件,说明一下采用智能优化算法求解双层优化问题的具体步骤:

1)设置算法参数

       采用智能优化算法求解上层优化问题时,由于上下层需同时优化,并进行迭代,为避免求解时间过长,种群规模和迭代次数不宜过大,在这里设置为30和50。另外,变量y的取值需要在下层优化中进行决策,上层优化只需要决策变量x,因此上层优化的粒子群算法问题维度设置为1。

%% 设置种群参数
sizepop = 30;                       % 初始种群个数
dim = 1;                            % 空间维数
ger = 50;                           % 最大迭代次数    
x_max = 14*ones(1,dim);             % 位置上限
x_min = zeros(1,dim);               % 位置下限
v_max = 7*ones(1,dim);              % 速度上限
v_min = -7*ones(1,dim);             % 速度下限
w=0.9;                              % 惯性权重
c_1 = 1.5;                          % 自我学习因子
c_2 = 1.5;                          % 群体学习因子

image.gif

2)初始化种群

%% 种群初始化
pop=x_min+rand(sizepop,dim).*(x_max-x_min);     % 初始化种群
pop_v=v_min+rand(sizepop,dim).*(v_max-v_min);   % 初始化种群速度        
pop_zbest=pop(1,:);                             % 初始化群体最优位置
pop_gbest=pop;                                  % 初始化个体最优位置
fitness=zeros(1,sizepop);                       % 所有个体的适应度
fitness_zbest=inf;                              % 初始化群体最优适应度
fitness_gbest=inf*ones(1,sizepop);              % 初始化个体最优适应度

image.gif

3)计算种群初始适应度

for k=1:sizepop
    % 计算适应度值
    fitness(k)=up_fitness_fun(pop(k,:));
    if fitness(k)<fitness_zbest
        fitness_zbest=fitness(k);
        pop_zbest=pop(k,:);
    end
end
history_pso=zeros(1,ger);            % 粒子群历史最优适应度值

image.gif

       粒子的适应度函数(也就是优化问题的目标函数取值)取决于变量x和y的值,而y的取值又由下层优化决定,因此求适应度函数的过程其实就是一个迭代过程,上层优化向下层优化传递变量x,下层优化根据x的取值和自身的目标函数以及约束条件,得到y的取值,然后向上层优化传递,最终根据x和y的取值综合得到目标函数,我们这里选择yalmip求解下层优化问题,代码如下:

%% up_fitness_fun.m
function fitness=up_fitness_fun(pop)
    x=pop;
    y=main_down(x);
    punish=0;
    if -2*x+3*y > 12
        punish=punish+999;
    end
    if x+y > 14
        punish=punish+999;
    end
    fitness=-x-2*y+punish;
end
%% main_down.m
function best_y=main_down(x)
    y=sdpvar(1);
    Constraints=[-3*x+y <= -3 , 3*x+y <= 30 ];
    objective=-y;
    ops=sdpsettings('verbose', 0 , 'solver', 'cplex');
    optimize(Constraints,objective,ops);
    best_y=value(y);
end

image.gif

4)迭代求最优解

       最终优化结果如下:

image.gif

image.gif

       求解结果和前两篇博客都一样(非常接近全局最优解x=8,y=6),但是因为每次迭代时,每个上层粒子都需要将自身的参数传递给下层粒子求解下层优化问题,意味着完成计算过程需要sizepop×ger=1500次下层优化,所以需要的时间很长。如果问题规模更大,模型更复杂,采用智能优化算法求解双层优化问题所需的时间将会更多,这时候就可以考虑将非线性双层优化问题模型简化为线性双层优化,然后再使用KKT条件转为单层优化问题进行求解。这样求解时间会大大缩短,可能也会提高计算精度(智能优化算法的局部收敛问题)



相关文章
|
4天前
|
存储 缓存 算法
如何通过优化算法和代码结构来提升易语言程序的执行效率?
如何通过优化算法和代码结构来提升易语言程序的执行效率?
|
4天前
|
机器学习/深度学习 人工智能 算法
[大语言模型-算法优化] 微调技术-LoRA算法原理及优化应用详解
[大语言模型-算法优化] 微调技术-LoRA算法原理及优化应用详解
22 0
[大语言模型-算法优化] 微调技术-LoRA算法原理及优化应用详解
|
5天前
|
搜索推荐
插入排序算法的讲解和代码
【10月更文挑战第12天】插入排序是一种基础的排序算法,理解和掌握它对于学习其他排序算法以及数据结构都具有重要意义。你可以通过实际操作和分析,进一步深入了解插入排序的特点和应用场景,以便在实际编程中更好地运用它。
|
7天前
|
机器学习/深度学习 算法 数据挖掘
基于GWO灰狼优化的GroupCNN分组卷积网络时间序列预测算法matlab仿真
本项目展示了基于分组卷积神经网络(GroupCNN)和灰狼优化(GWO)的时间序列回归预测算法。算法运行效果良好,无水印展示。使用Matlab2022a开发,提供完整代码及详细中文注释。GroupCNN通过分组卷积减少计算成本,GWO则优化超参数,提高预测性能。项目包含操作步骤视频,方便用户快速上手。
|
8天前
|
机器学习/深度学习 算法 数据安全/隐私保护
基于WOA鲸鱼优化的GroupCNN分组卷积网络时间序列预测算法matlab仿真
本项目展示了一种基于WOA优化的GroupCNN分组卷积网络时间序列预测算法。使用Matlab2022a开发,提供无水印运行效果预览及核心代码(含中文注释)。算法通过WOA优化网络结构与超参数,结合分组卷积技术,有效提升预测精度与效率。分组卷积减少了计算成本,而WOA则模拟鲸鱼捕食行为进行优化,适用于多种连续优化问题。
|
1天前
|
机器学习/深度学习 算法 数据安全/隐私保护
基于贝叶斯优化卷积神经网络(Bayes-CNN)的多因子数据分类识别算法matlab仿真
本项目展示了贝叶斯优化在CNN中的应用,包括优化过程、训练与识别效果对比,以及标准CNN的识别结果。使用Matlab2022a开发,提供完整代码及视频教程。贝叶斯优化通过构建代理模型指导超参数优化,显著提升模型性能,适用于复杂数据分类任务。
|
7天前
|
算法 决策智能
基于GA-PSO遗传粒子群混合优化算法的TSP问题求解matlab仿真
本文介绍了基于GA-PSO遗传粒子群混合优化算法解决旅行商问题(TSP)的方法。TSP旨在寻找访问一系列城市并返回起点的最短路径,属于NP难问题。文中详细阐述了遗传算法(GA)和粒子群优化算法(PSO)的基本原理及其在TSP中的应用,展示了如何通过编码、选择、交叉、变异及速度和位置更新等操作优化路径。算法在MATLAB2022a上实现,实验结果表明该方法能有效提高求解效率和解的质量。
|
17天前
|
机器学习/深度学习 算法 数据安全/隐私保护
基于MSER和HOG特征提取的SVM交通标志检测和识别算法matlab仿真
### 算法简介 1. **算法运行效果图预览**:展示算法效果,完整程序运行后无水印。 2. **算法运行软件版本**:Matlab 2017b。 3. **部分核心程序**:完整版代码包含中文注释及操作步骤视频。 4. **算法理论概述**: - **MSER**:用于检测显著区域,提取图像中稳定区域,适用于光照变化下的交通标志检测。 - **HOG特征提取**:通过计算图像小区域的梯度直方图捕捉局部纹理信息,用于物体检测。 - **SVM**:寻找最大化间隔的超平面以分类样本。 整个算法流程图见下图。
|
1天前
|
存储
基于遗传算法的智能天线最佳阵列因子计算matlab仿真
本课题探讨基于遗传算法优化智能天线阵列因子,以提升无线通信系统性能,包括信号质量、干扰抑制及定位精度。通过MATLAB2022a实现的核心程序,展示了遗传算法在寻找最优阵列因子上的应用,显著改善了天线接收功率。
|
4天前
|
监控 算法 数据安全/隐私保护
基于三帧差算法的运动目标检测系统FPGA实现,包含testbench和MATLAB辅助验证程序
本项目展示了基于FPGA与MATLAB实现的三帧差算法运动目标检测。使用Vivado 2019.2和MATLAB 2022a开发环境,通过对比连续三帧图像的像素值变化,有效识别运动区域。项目包括完整无水印的运行效果预览、详细中文注释的代码及操作步骤视频,适合学习和研究。

热门文章

最新文章