(文章复现)基于改进秃鹰算法的微电网群经济优化调度研究

参考文献：

[1]周辉,张玉,肖烈禧,等.基于改进秃鹰算法的微电网群经济优化调度研究[J].太阳能学报,2024,45(02):328-335.

1.基本原理

        微电网群由3个独立的微电网(microgrid , MG)组成,各微电网内部包含光伏(photovoltaic , PV)、风力发电机(windturbine,WT)、电动汽车(electric vehicle,EV)、微型燃气轮机( microturbine, MT)、储能电池(battery , BT)及负荷,通过能量管理中心实现彼此之间电能交互及向配电网买卖电,如图1所示。

1.1 微电网系统模型

1.1.1 微型燃气轮机

1.2.2 储能系统

1.2.3 电动汽车负荷

1.2目标函数

        考虑微电网群系统运行成本及环境成本，构建微电网群经济优化调度模型，其中运行成本包括可控分布式发电单元的发电成本、设备运行维护成本、BT 运行成本、电能交易成本；环境成本为 CO2、SO2、NOx污染物气体的排放惩罚成本，以微电网群系统总运行成本为目标函数进行优化调度。

1.2.1 运行成本

1.2.2 环境成本

1.3约束条件

1.3.1 系统功率平衡约束

1.4.2 MT 爬坡功率约束

1.4.3 BT 约束

1.4.4 MG 之间交互功率约束

1.4.5 MG 与配电网之间交互功率约束

2.求解方法 

2.1基本秃鹰搜索算法

        BES 是一种针对秃鹰狩猎行为提出的自然启发式算法，包含选择搜索空间、搜索猎物及俯冲3个阶段 。

2.1.1 选择搜索空间阶段

        该阶段秃鹰根据猎物数量随机选择搜索区域，飞到当前最优个体附近，计算公式为：

2.1.2 搜索猎物阶段

        该阶段秃鹰在选定的搜索空间内螺旋飞行，加速搜索猎物，计算公式为：

2.1.3 俯冲阶段

        该阶段秃鹰从当前最优位置以螺旋飞行的方式冲向猎物，计算公式为：

2.2改进秃鹰搜索算法

2.2.1 反向学习

        反向学习可根据当前解获得相对中心的反向解，通过对比二者保留更好的解，从而扩大算法寻优范围，计算公式为：

2.2.2 柯西变异

        柯西分布产生的随机数分布范围较广，可有效提高算法的全局搜索能力，跳出局部最优。柯西分布函数为：

2.2.3 融合反向学习和柯西变异策略

        综上所述，在秃鹰算法搜索猎物空间步骤中，通过反向学习扩大秃鹰算法搜索范围，并采用柯西变异跳出局部最优，两者通过选择概率 P 来决定执行，计算公式为：

2.2.4 算法步骤

        综上所述，本文提出的 IBES 流程如图 2 所示。

3.编程思路分析

3.1相关参数和决策变量定义

表1 相关参数的定义

表2 决策变量的定义

3.2编程思路

        根据对文献内容的解读，可以设计下面的编程思路：

步骤1：输入所需数据

        这一步比较简单。算例分析用到的大部分数据可以从原文中找到，其他数据可以自己假设一下。然后将所有需要的数据，按照表1的定义格式输入即可。

步骤2：改进秃鹰搜索算法的实现

        实现基本的秃鹰搜索算法和改进秃鹰搜索算法，并使用文中提到的四个测试函数进行验证。文中所提供基准测试函数中，f3和f4都存在问题，其中测试函数f3的公式表述不清，测试函数f4没有提供参数aij的取值。因此我在标准测试函数中随便找了两个当作测试函数。

步骤3：将改进秃鹰搜索算法用于微网调度问题

        如表2所示，原文中共包含4个决策变量，都是3×24的变量，也就是说决策变量的总数是3×24×4=288，需要将其组合成一个1×288的变量，便于使用改进秃鹰搜索算法求解，也就意味中文中涉及的优化问题维度为288，具体变量设置如下：

        智能优化算法中对于约束条件的处理有很多种形式，我在代码中使用罚函数法进行处理，即将不满足约束的部分作为惩罚添加到目标函数中。

        此外，文献中忽略了MG之间的交互功率总和需要为0这一约束，也就是下面的公式：

        不加上这一项约束的话，MG之间交互功率即使都为正，也是满足文中所提到的约束，但很明显是不符合实际情况的。另外，如果只遵循上下限约束生成初始种群，可能导致生成的种群中大部分甚至全部的个体都是不满足约束的，因此生成种群和更新种群时，还是要满足一定规则，使得生成的种群尽可能都满足约束。

步骤4：输出运行结果

        参考文中的图表的格式，输出结果即可。

4.Matlab代码

%% 设置种群参数
sizepop = 40;                       % 初始种群个数
dim = 10;                           % 空间维数
ger = 500;                          % 最大迭代次数   
% 选择适应度函数
[xmax0, xmin0, fun] = Select_fitnessfun(fitnesfun);
x_max = xmax0*ones(1,dim);          % 位置上限
x_min = xmin0*ones(1,dim);          % 位置下限
v_max = x_max*0.5;                  % 速度上限
v_min = -x_max*0.5;                 % 速度下限
a = 2;                              % 位置变化参数
a1 = 10;                            % 搜索点之间角的参数
R = 1.5;                            % 搜索周期数
c1 = 1.8;                           % 增加秃鹰移动强度的随机数
c2 = 1.8;                           % 增加秃鹰移动强度的随机数
[x,y] = polr(a,R,sizepop);          % 搜索猎物阶段的参数    
[x1,y1]=swoo_p(a,R,sizepop);        % 俯冲阶段的参数

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

%% 初始的适应度
for k = 1:sizepop
    % 计算适应度值
    fitness(k) = fun(pop(k,:));
    if fitness(k) < fitness_best
        fitness_best = fitness(k);
        pop_best = pop(k,:);
    end
end
history_BES = zeros(1,ger);                         % BES历史最优适应度值
%% 迭代求最优解
iter = 1;
while iter <= ger
    for k = 1:sizepop
        % 1.选择搜索空间阶段
        pop_new = pop_best + 2*rand(1,dim).*(mean(pop) - pop(k,:));
        fitness_new = fun(pop_new);
        if fitness_new < fitness(k)
            fitness(k) = fitness_new;
            pop(k,:) = pop_new;
        end
        
        % 2.搜索猎物阶段
        if k <= sizepop - 1
            pop_new = pop(k,:) + y(k)*(pop(k,:) - pop(k + 1,:)) + x(k)*(pop(k,:) - mean(pop));
            fitness_new = fun(pop_new);
            if fitness_new < fitness(k)
                fitness(k) = fitness_new;
                pop(k,:) = pop_new;
            end
        end
        
        % 3.俯冲阶段
        pop_new = rand(1,dim).*pop_best + x(k)*(pop(k,:) - c1*mean(pop)) + y(k)*(pop(k,:) - c2*pop_best);
        fitness_new = fun(pop_new);
        if fitness_new < fitness(k)
            fitness(k) = fitness_new;
            pop(k,:) = pop_new;
        end
        
        % 更新位置并对位置进行边界处理
        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
        
        % 更新种群最优
        if fitness(k) < fitness_best
            fitness_best = fitness(k);
            pop_best = pop(k,:);
        end
    end
    history_BES(iter) = fitness_best;
    iter = iter+1;
end
% time0 = toc;
% disp(['运行时间为：',num2str(time0) , '秒'])
% disp(['最优解：x=',num2str(pop_best)])
% disp(['最优函数值=',num2str(fitness_best)])
% plot(history_BES,'linewidth',1)
% ylabel('最优适应度值')
% xlabel('迭代次数')


function [xR,yR] = swoo_p(a,R,N)
th = a*pi*exp(rand(N,1));
xR = R.*sinh(th);
yR = R.*cosh(th);
xR = xR/max(abs(xR));
yR = yR/max(abs(yR));
end

function [xR,yR] = polr(a,R,N)
th = a*pi*rand(N,1);
r = th+R*rand(N,1);
xR = r.*sin(th);
yR = r.*cos(th);
xR = xR/max(abs(xR));
yR = yR/max(abs(yR));
end

        以上仅为基本秃鹰搜索算法的matlab代码主函数部分，完整论文复现的matlab代码可以从这个链接获取：

(文章复现)基于改进秃鹰算法的微电网群经济优化调度研究的matlab代码资源-CSDN文库

5.运行结果