参考文献：

[1]高瑜,黄森,陈刘鑫等.基于改进灰狼算法的并网交流微电网经济优化调度[J].科学技术与工程, 2020,20(28):11605-11611.

[2]邓长征,冯朕,邱立等.基于混沌灰狼算法的交直流混合微网经济调度[J].电测与仪表, 2020, 57(04):99-107.

这两篇文章不管是从模型、原理、求解方法还是算例分析上看，都非常相似，鉴于文章两篇文章的作者单位没有重合，肯定是两位不同的作者在不同的时间灵感爆发，提出了相同的idea(甚至文献1中目标函数中成本只有4项，但算例分析部分和文献2一样有5项成本)，证明了学术论文的可重复性。分析表明，这两篇文章绝对没有谁照搬谁的情况，也不可能是两个人抄了同一篇英文文献，更不可能是找了同一家论文代写机构然后被坑了。

综上所示，这篇博客就是对上面两篇不同文章中提到的一种相同的方法进行复现，以实现本人水一篇博客的目的。

1.基本原理

交直流混合微电网运行方式分为孤岛运行方式和并网运行方式，在本次研究中主要考虑并网情况下交直流混合微电网的经济优化调度。交直流混合微电网分为交流侧和直流侧，交流侧与外部大电网互联。交流侧包含交流母线、交流微电源和交流负荷; 直流侧包含直流母线、直流微电源和直流负荷。具体结构如图 1所示。

交直流混合微电网的主要运行方式有4种:

( 1) 满足交流侧和直流侧功率平衡，交流侧和直流侧的微电源相应满足交流负荷和直流负荷的供电需求;

( 2) 通过AC /DC 双向换流器使交流侧和直流侧功率交互，交流侧和直流侧互为补充，满足微网内交流负荷和直流负荷总的供电需求;

( 3) 直流侧储能系统参与微电网削峰填谷和平抑功率，当微网中微电源总出力大于负荷所需供电要求时; 储能系统充电，反之则储能系统放电;

( 4) 当微网中微电源和储能系统在某一时刻的出力不足以满足微网中负荷的需求时，向交流大电网购电; 反之则向大电网售电。储能系统和换流系统在以上运行方式中有重要的作用，因此在建立交直流混合微电网经济优化调度的经济模型中应考虑储能损耗和换流成本。

1.1目标函数

交流微电网内器件维修保养费用 F2 为

1.2约束条件

2.改进灰狼算法求解

2.1基本灰狼算法

来自澳大利亚的研究者 Mirjalili 于2014年提出一种新型群体智能算法——灰狼算法( grey wolf optimizer，GWO) ［12］。灰狼算法通过模拟灰狼狼群在捕食过程中的特征及狼群等级制度，通过灰狼对猎物不停地搜索追寻，不断更新猎物所处的位置，最终灰狼成功捕捉猎物。在 GWO 算法中，将狼群按照等级由高到低分为 α 狼、β 狼、δ 狼和其余狼群个体 ω，灰狼的等级越高，其自身适应度越好，猎物的位置则对应全局最优解。α 狼为头狼，在狼群中起领导及决策作用; β 狼为头狼候选，在狼群中反馈信息给 α 狼; δ 狼负责对 ω 狼群的统领; ω 狼对猎物实施搜寻与追捕。狼群捕食过程首先要包围猎物，灰狼需要知道自身与猎物的距离大小，然后根据自身与猎物的距离来调整自己的位置。灰狼包围猎物的数学公式为

在狼群捕食过程中，α 狼、β 狼、δ 狼与猎物之间的距离可由式( 25) ～式( 27) 计算。由于 α 狼、β狼、δ 狼在狼群中等级制度，这意味着它们的适应度更好，从而离猎物的距离更近，ω 狼根据 α 狼、β 狼、 δ 狼的位置向猎物靠近，ω 狼朝向 α 狼、β 狼、δ 狼前进的方向与前进步长分别可由式( 28) ～式( 30) 计算。

灰狼算法的流程图如图 2 所示。

2.2改进灰狼算法

均衡 GWO 的全局搜索能力和局部搜索能力，是 GWO 能否取得高寻优性能的关键因素。当| A | ＞ 1 时，灰狼群体将扩大包围圈，以寻找更让好的猎物，此时对应于全局搜索; 当 | A | ＜ 1 时，灰狼群体将收缩包围圈，从而对猎物完成最后的攻击行为，此时对应于局部精确搜索［13］。因此 GWO 的全局搜索能力和局部搜索能力与 | A | 的取值之间有着极大关联，另外根据式( 23) 可知 | A| 的取值也与收敛因子 a 的取值紧密相关。由于原始 GWO 的收敛因子是线性的，且收敛因子 a 从 2 线性递减至 0，这种线性递减策略不能完全体现出实际的优化搜索过程。因此，引进一种基于非线性变化的收敛因子更新公式:

采取非线性变化的收敛因子在寻优初始阶段，a 衰减程度低，这时候，算法的全局寻优能力较强，能够搜索到更好的全局最优解;寻优末尾阶段，a 衰减程度高，这时候，算法的局部寻优能力较强，能够搜索到更好的局部最优解。这种基于非线性变化的收敛因子可以更够更好地均衡 GWO 的全局寻优能力和局部寻优能力，从而获得精度更好的解。

3.编程思路分析

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

表1 相关参数的定义

表2 决策变量的定义

3.2编程思路

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

步骤1：输入所需数据

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

步骤2：灰狼算法的实现

实现基本的灰狼优化算法和改进灰狼优化算法，并使用文中提到的四个测试函数进行验证，两篇文章中对灰狼算法改进各不相同，一个是改进了收敛因子，一个是改进了初始化策略，这篇博客中的代码是将两种改进策略结合起来。以Sphere函数为例，基本灰狼优化算法和改进灰狼优化算法的对比如下(优化问题维度取10，种群数和最大迭代次数分别为200，100)：

基本灰狼算法的运行结果：

最优解：x=1.6921e-13 9.8087e-15 -1.527e-13 -3.2711e-15 -3.2045e-14 -7.4947e-14 1.5045e-14 1.6444e-14 -2.0382e-14 2.4686e-14

最优函数值=6.0224e-26

改进灰狼算法的运行结果：

最优解：x=1.5078e-17 -1.7022e-16 1.3285e-18 3.5724e-17 2.0807e-16 3.5095e-17 -2.8155e-16 2.14e-16 1.5832e-16 -9.8643e-17

最优函数值=2.3486e-31

仅针对Sphere函数而言，改进灰狼算法的效果要更好一些。

步骤3：将灰狼算法用于微网调度问题

如表2所示，原文中共包含3个决策变量，都是1×24的变量，需要将其组合成一个1×72的变量，便于使用灰狼算法求解，也就意味中文中涉及的优化问题维度为72，具体变量设置如下：

智能优化算法中对于约束条件的处理有很多种形式，我在代码中参考这篇文章给的罚函数法进行处理：粒子群算法求解带约束优化问题 - 知乎。

另外，如果只遵循上下限约束生成初始种群，可能导致生成的种群中大部分甚至全部的个体都是不满足约束的，因此生成种群和更新种群时，还是要满足一定规则，使得生成的种群尽可能都满足约束。

步骤4：输出运行结果

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

4.Matlab代码

%% 清除变量
clc
clear
close all
warning off

%% 设置种群参数
sizepop = 200;                      % 初始种群个数
dim = 10;                           % 空间维数
ger = 100;                          % 最大迭代次数
x_max = 100*ones(1,dim);            % 位置上限
x_min = -100*ones(1,dim);           % 位置下限

%% 种群初始化
pop = x_min + rand(sizepop,dim).*(x_max-x_min);     % 初始化种群
fitness = zeros(1,sizepop);                         % 所有个体的适应度
for k = 1:sizepop
    fitness(k) = Sphere(pop(k,:));
end

% 初始化Alpha、Beta、Delta狼群
[SortFitness,Index] = sort(fitness);
A_pop = pop(Index(1),:);                            % 初始化α狼群
A_fitness = SortFitness(1);                         % 初始化α狼群适应度
B_pop = pop(Index(2),:);                            % 初始化β狼群
B_fitness = SortFitness(2);                         % 初始化β狼群适应度
D_pop = pop(Index(3),:);                            % 初始化δ狼群
D_fitness = SortFitness(3);                         % 初始化δ狼群适应度
history = zeros(1,ger);                             % 历史最优适应度值
%% 迭代求最优解
iter = 1;
while iter <= ger
    a = 2 - iter*(2/ger);                           % 线性调整a的值
    for k = 1:sizepop
        % 1.根据Alpha狼群更新位置X1
        r1 = rand;
        r2 = rand;
        A1 = 2*a*r1 - a;
        C1 = 2*r2;
        D_A = abs(C1*A_pop - pop(k,:));
        X1 = A_pop - A1*D_A;
        
        % 2.根据Beta狼群更新位置X2
        r1 = rand;
        r2 = rand;
        A2 = 2*a*r1 - a;
        C2 = 2*r2;
        D_B = abs(C2*B_pop - pop(k,:));
        X2 = B_pop - A2*D_B;
        
        % 3.根据Delta狼群更新位置X3
        r1 = rand;
        r2 = rand;
        A3 = 2*a*r1 - a;
        C3 = 2*r2;
        D_D = abs(C3*D_pop - pop(k,:));
        X3 = D_pop - A3*D_D;
        
        % 4.更新后的狼群位置
        pop(k,:) = (X1+X2+X3)/3;
        
        % 修复越限的狼群位置
        pop(k,pop(k,:) > x_max) = x_max(pop(k,:) > x_max);
        pop(k,pop(k,:) < x_min) = x_max(pop(k,:) < x_min);
        
        % 更新狼群的适应度
        fitness(k) = Sphere(pop(k,:));
    end
    % 更新Alpha、Beta、Delta狼群位置
    [SortFitness,Index] = sort(fitness);
    A_pop = pop(Index(1),:);                            % 初始化α狼群
    A_fitness = SortFitness(1);                         % 初始化α狼群适应度
    B_pop = pop(Index(2),:);                            % 初始化β狼群
    B_fitness = SortFitness(2);                         % 初始化β狼群适应度
    D_pop = pop(Index(3),:);                            % 初始化δ狼群
    D_fitness = SortFitness(3);                         % 初始化δ狼群适应度
    history(iter) = A_fitness;
    iter = iter+1;
end

%% 画图
disp(['最优解：x=',num2str(A_pop)])
disp(['最优函数值=',num2str(A_fitness)])
plot(history,'linewidth',1)
ylabel('最优适应度值')
xlabel('迭代次数')
%% 适应度函数
function fitness = Sphere(pop)
fitness = sum(pop.^2);
end

运行结果如下：

最优解：x=-2.7836e-14 1.1833e-13 2.0262e-13 1.0821e-14 -2.534e-14 6.7367e-14 -7.4143e-14 -1.6845e-13 -1.5746e-13 3.9423e-14
最优函数值=1.2135e-25