未发表！基于主从博弈的多主体微网/综合能源系统利益分摊与能量管理程序代码！

适用平台：Matlab+Yalmip+Cplex

程序首先建立了综合能源运营商、新能源用户、电动汽车充电商的多主体园区，分析了三种市场交易主体的属性和市场交易机制。其次建立基于主从博弈的三方市场主体能量管理与收益分摊方案，利用改进粒子群算法进行求解。程序中算例丰富、注释清晰、干货满满，创新性很高！下面对文章和程序做简要介绍！

程序创新点：

（1）建立基于主从博弈的多方市场主体能量管理策略，解决多主体之间的多目标优化问题；

（2）电动汽车充电商通过管理车辆的充电行为，以降低自身充电充本；

（3）将新能源与用户划分为一体，可向电动汽车外售电能获取收益，也可购电满足自身负荷需求，具有能源生产/消费者双重属性；

主要工作：

作为耦合电、热、冷、气、交通等多能源网络的底层多能互补集成耦合系统，园区综合能源系统的能源利用效率相比传统的电能园区系统较优，在可再生能源就地消纳、提高需求侧调度灵活性、实现多能互补协同利用等方面具有重要的现实意义。多元市场主体参与市场化交易必然需要制定合适的交易机制，以协调不同市场主体的利益冲突。而博弈论能够解决复杂主体之间的多目标优化问题，实现系统内的最优综合能量管理。程序首先建立了含系统能源运营商、含分布式光伏的用户、电动汽车充电代理商的园区综合能源系统，分析了三种市场交易主体的属性和市场交易机制。其次建立了三方市场主体各自综合能量管理优化运行方案，并运用改进的粒子群算法进行求解。

文章框架：

交易主体及交易机制分析：

园区综合能源系统内的市场交易主体类型主要包括：能源生产者、能源消费者。两者均以自身利益最大化为目标来制定能量产消策略，实现自身能量优化管理。本系统内一共有能源运营商、含分布式光伏的用户、电动汽车充电代理商等3个市场主体，分别在系统内扮演不同的角色。三个市场交易主体构成了以系统能源运营商为主导的三方主从博弈模型。在该非合作交易中，三方均理性追求自身收益最大化，并制定相应的售能和用能策略，三方互动和交流最终达到一种均衡的最佳交易和运行状态。

综合能源运营商：

系统能源运营商是其他市场交易主体的主要供能方。系统能源运营商能量优化分配过程中，以自身净收益最大作为目标，对于含分布式光伏的用户而言，当自身供能不足时，能源运营商能够向其供能以满足额外用能需求。对电动汽车充电代理商而言，能源运营商既可以在较低电价时段出售足够的电能，以满足其较低成本的充电需求，也可以在用户低报价时段出售不足的光伏电能情况下，提供额外电能作为补充。

新能源用户能量管理策略：

含分布式光伏的用户是系统内的一类新兴市场交易主体，广泛参与系统综合能源管理，以自身净收益最大作为目标，具体有两个途径扩大收益：一方面优化自身能源转换设备出力降低运行成本，用户的分布式光伏还享受政府的全电量补贴政策，根据实际的发电量可获取一定的补贴收益；另一方面当用户的光伏出力满足自身需求且有盈余时，可通过制定合理且具有竞争力的报价策略，获得电动汽车充电代理商的售电权，获得售电收益。

电动汽车充电商：

电动汽车充电代理商是系统内的纯电能负荷，通过控制内部车辆的充电电量和充电时间参与系统综合能量管理，以自身运行成本最低作为目标，其运行成本的降低可以通过两个途径来减少：一方面选择比低报价市场主体购得电能，另一方面管理车辆充放电行为减少充电成本。

程序结果：

部分程序：

%% 产生初始粒子和速度
%   k
for j=1:sizepop
    %随机产生一个种群
    for i=1:N     %初始化粒子的维数
        %%%%%%%%%%%%使得角度在（0,2π上均匀分布）%%%%%%%%%%%%%%
        pop(((j-1)*N+i))=rand;    %初始种群位置,即节点的角度
        V(((j-1)*N+i),1)=alfa*randn(1);     %初始化速度,即旋转的角度,每次旋转一度
    end
end
scale=0;
%% 个体极值和群体极值
zbest=pop(1:N,:);   %全局最佳
gbest=pop;    %个体最佳
fitnessgbest=scale(1)*ones(1,sizepop);   %个体最佳适应度值
fitzbest=scale(1);   %全局最佳适应度值
s(1)=fitzbest;
%%%%%%%%%%%%%%%%%%%%%%%%%%
%%需要改，主感知方向作为粒子的位置，用速度表示主感知方向变换的角度，最后统一将调整方向，质心的受力主要用来防止粒子过于密集
%%%%%%%%%%%%%%%%%%%%%%%%%%
for K=2:maxgen     %迭代更新的次数
    w=1; 
    for j=1:sizepop   % 20个粒子
        for i=1:N  %粒子的维数
            %速度更新
            V(((j-1)*N+i),:) = w*V(((j-1)*N+i),:) + c1*rand*(gbest(((j-1)*N+i),:) - pop(((j-1)*N+i),:)) + c2*rand*(zbest(i,:) - pop(((j-1)*N+i),:));
            V(((j-1)*N+i),V(((j-1)*N+i),:)>Vmax)=Vmax;
            V(((j-1)*N+i),V(((j-1)*N+i),:)<Vmin)=Vmin;  %V为质心绕节点顺时针旋转的角度
            
            %种群更新
            pop(((j-1)*N+i),:)=pop(((j-1)*N+i),:)+V(((j-1)*N+i),:);
            pop(((j-1)*N+i),:) = min(max((pop(((j-1)*N+i),:)),0),1);
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        end
        scale(j)=Obj(pop(((j-1)*N+1):((j-1)*N+N),:));
        fitness(j)=scale(j);   
        
        %个体最优更新
        if fitness(j) > fitnessgbest(j)
            gbest(((j-1)*N+1):((j-1)*N+N),:) = pop(((j-1)*N+1):((j-1)*N+N),:);
            fitnessgbest(j) = fitness(j); 
        end
        
        %群体最优更新
        if fitness(j) > fitzbest
            zbest = pop(((j-1)*N+1):((j-1)*N+N),:);
            fitzbest = fitness(j);
        end
    end
    s(K)=fitzbest; 
%     K
    fitzbest
end

欢迎感兴趣的小伙伴关注，小编会不定期更新高质量的学习资料、文章和程序代码，为您的科研加油助力！