（MATLAB代码分享、可运行）基于NSGA-2算法的多目标水电站电力调度优化

问题描述

光伏发电的随机性、波动性和间歇性限制了电网对光伏发电的接纳能力，随着光伏发电的快速发展，提高光伏的消纳能力是一个亟待解决的大问题。利用水电和光伏进行联合发电是一条提高光伏消纳能力新的途径。本文建立了短期的优化调度模型，以调度期内调峰能力最大和互补系统的出力和负荷需求的偏差最小为目标函数，采用模拟退火粒子群算法进行求解，处理多目标问题时采用拥挤距离排序，并采用外部档案的方法进行数据维护，最终得出水光互补可以明显提高系统的调峰能力。

建模分析

1.Pw：水电出力 约束：0Mw≤Pw≤891Mw 单位：Mw

2.Pp：光伏出力 约束：0≤Pp≤277（已知值） 单位：Mw

时间/h	1	2	3	4	5	6	7	8	9	10	11	12
水电出力（Mw）	0	0	0	0	0	0	2	31	97	169	192	258
时间/h	13	14	15	16	17	18	19	20	21	22	23	24
光伏出力（Mw）	274	277	266	226	176	79	38	9	0	0	0	0

3.水电电量约束：个人感觉应该为24个小时，水电所发的电量之和为定值

具体数值：24个小时水电出力Pw相加17803Mwh

∆t=1h，单位一个小时

4. 功率平衡约束：Pp+pw≤Pd

其中：

(1)Pp：光伏出力，即24个小时的光伏 单位：Mw

时间/h	1	2	3	4	5	6	7	8	9	10	11	12
水电出力（Mw）	0	0	0	0	0	0	2	31	97	169	192	258
时间/h	13	14	15	16	17	18	19	20	21	22	23	24
光伏出力（Mw）	274	277	266	226	176	79	38	9	0	0	0	0

（2）Pd ：发电计划，即24个小时所对应的发电计划 单位：Mw

时间/h	1	2	3	4	5	6	7	8	9	10	11	12
发电计划（Mw）	750	780	770	788	756	728	867	899	1098	1256	1340	1250
时间/h	13	14	15	16	17	18	19	20	21	22	23	24
发电计划（Mw）	988	580	572	601	656	782	890	850	760	750	768	729

（3）Pw：水电出力 单位：Mw

由发电的流量决定，而发电流量由水位决定，水头和水位不是同一概念。（仅供参考）我认为决策变量只有水位。在粒子群算法中首先生成水位序列（随机位置），然后生成水位变化序列（进化的速度）。进而由水位所对应的发电流量求出水电的出力，就可以计算出两个目标函数）（仅供参考）

5. Vt：水量 单位亿m³ 这个理论上就是库容约束，可忽略。

6. Vt：库容 约束：0亿m³-----2862.8亿m³

7.Zt：水位 约束：2242m-----2469.15m 单位米

8.Qt：发电流量 约束：200m/s-----7000m/s 单位米/秒

9. Ht：水头高度 约束：99.5m-------150m 单位米（计算时用算术平均值124.75米）

MATLAB程序

代码运行过程示意图

结果示意图

main.m主程序如下：（完整代码可评论私信留邮箱）

clc;
clear;
close all;

%% Problem Definition

CostFunction=@(x) myfun(x);      % Cost Function

nVar=24;             % 时间段

VarSize=[1 nVar];   % Size of Decision Variables Matrix

VarMin= 2451.5.*ones(1,24);        % 水位下界 %这里为了产生初始解做了相关调整
VarMax= 2455.*ones(1,24);          % 水位上界

% Number of Objective Functions
% nObj=numel(CostFunction(unifrnd(VarMin,VarMax,VarSize)));
nObj = 2;

%% NSGA-II Parameters

MaxIt=70;      % Maximum Number of Iterations

nPop=80;        % Population Size

pCrossover=0.7;                         % Crossover Percentage
nCrossover=2*round(pCrossover*nPop/2);  % Number of Parnets (Offsprings)

pMutation=0.4;                          % Mutation Percentage
nMutation=round(pMutation*nPop);        % Number of Mutants

mu=0.02;                    % Mutation Rate

sigma=0.1*(VarMax-VarMin);  % Mutation Step Size


%% Initialization

empty_individual.Position=[];
empty_individual.Cost=[];
empty_individual.Rank=[];
empty_individual.DominationSet=[];
empty_individual.DominatedCount=[];
empty_individual.CrowdingDistance=[];

pop=repmat(empty_individual,nPop,1);

disp('产生初始可行解...')
for i=1:nPop
    
    flag=0;
    while flag==0
        tmp=[];
        for j=1:1:nVar
            tmp = [tmp unifrnd(VarMin(j),VarMax(j),1)];
        end
        flag = test(tmp); % 检查约束 约束不满足就重新生成解
    end
    pop(i).Position=tmp;
    
    pop(i).Cost=CostFunction(pop(i).Position);
    
end

% pause

% Non-Dominated Sorting
[pop, F]=NonDominatedSorting(pop);

% Calculate Crowding Distance
pop=CalcCrowdingDistance(pop,F);

% Sort Population
[pop, F]=SortPopulation(pop);


%% NSGA-II Main Loop

for it=1:MaxIt
    
    % 交叉
    popc=repmat(empty_individual,nCrossover/2,2);
    for k=1:nCrossover/2
        
        i1=randi([1 nPop]);
        p1=pop(i1);
        
        i2=randi([1 nPop]);
        p2=pop(i2);
        
        [popc(k,1).Position, popc(k,2).Position]=Crossover(p1.Position,p2.Position,VarMin,VarMax);
        if test(popc(k,1).Position)+test(popc(k,2).Position)==2
            popc(k,1).Cost=CostFunction(popc(k,1).Position);
            popc(k,2).Cost=CostFunction(popc(k,2).Position);
        else
            popc(k,1)=p1;
            popc(k,2)=p2;
        end
        
    end
    popc=popc(:);
    
    % 变异
    popm=repmat(empty_individual,nMutation,1);
    for k=1:nMutation
        
        i=randi([1 nPop]);
        p=pop(i);
        
        popm(k).Position=Mutate(p.Position,mu,sigma,VarMin,VarMax);
        
        %% 越界处理
%         chrom=[popm(k).Position];
%         chrom(find(chrom<VarMin(1)))=VarMin(1);
%         chrom(find(chrom>VarMax(1)))=VarMax(1);
%         popm(k).Position=chrom;
        
        if test(popm(k).Position)
            
            popm(k).Cost=CostFunction(popm(k).Position);
        else
            popm(k)=p;
        end
        
    end
    
    % 合并父代和子代
    pop=[pop
        popc
        popm]; %#ok
    
    % Non-Dominated Sorting
    [pop, F]=NonDominatedSorting(pop);
    
    % Calculate Crowding Distance
    pop=CalcCrowdingDistance(pop,F);
    
    % Sort Population
    pop=SortPopulation(pop);
    
    % 截取前n个解
    pop=pop(1:nPop);
    
    % Non-Dominated Sorting
    [pop, F]=NonDominatedSorting(pop);
    
    % Calculate Crowding Distance
    pop=CalcCrowdingDistance(pop,F);
    
    % Sort Population
    [pop, F]=SortPopulation(pop);
    
    % Store F1
    F1=pop(F{1});
    
    % Show Iteration Information
    disp(['Iteration ' num2str(it) ': Number of F1 Members = ' num2str(numel(F1))]);
    
    % Plot F1 Costs
    figure(1);
    PlotCosts(F1);
    pause(0.01);
    
end

%% Results
result

最后，博主专注于论文的复现工作，有兴趣的同学可以私信共同探讨。相关代码已经上传到资源共享，点击我的空间查看分享代码。