热门标签
热门文章
- 1【科学计算包NumPy】NumPy数组的基本操作_numpy库对数组的操作与运算
- 2云原生最佳实践系列 5：基于函数计算 FC 实现阿里云 Kafka 消息内容控制 MongoDB DML 操作_函数计算 fc,创建函数,这里能够配置我云服务器中自己搭建的mongodb吗
- 3码农的自我修养 - ARM编译器的区分_宏定义判断是不是mingw
- 4iptables使用limit控制新建连接速率_iptables limit
- 5turbo码_信息论-Turbo码学习
- 6【已解决】git 设置不参与合并的文件_git不合并某个文件
- 79款最新文生图模型汇总！含华为、谷歌、Stability AI等大厂创新模型（附论文和代码）_文生图大模型有哪些
- 8Java获取当前时间的几种方式_java获取当前日期
- 9二、buildroot-2015.02编译根文件系统支持Qt_buildroot 2015
- 10（附源码）SSM学生宿舍管理系统JAVA计算机毕业设计项目_宿舍管理系统怎么在其他电脑运行
当前位置: article > 正文
论文复现-基于NSGA-Ⅲ的梯级水电-火电机组多目标优化调度（全代码-约3000行）_nsga3能改的数据
作者：羊村懒王 | 2024-05-24 11:47:44
赞
踩
nsga3能改的数据
火电是目前我国电网中提供调节能力的主要能源，而随着碳中和目标的提出，势必要降低火电在电网中的分量。梯级水电则是可以替代火电发挥调节作用的能源。基于电网中存在大量火电的背景，梯级水电如何与火电机组耦合调度呢？本代码通过优于NSGA-Ⅱ的NSGA-Ⅲ优化算法实现了梯级水电和火电机组的联合多目标调度。
参考文献：
【1】An interactive fuzzy satisfying method based on evolutionary programming technique for multi-objective short-term hydrothermal scheduling. Electric Power Systems Research；
【2】An extended NSGA-III for solution multi-objective hydro-thermal-wind scheduling considering wind power cost. Energy Conversion and Management.
订阅专栏即可获取全代码，代码较为复杂，调试期间有疑问可私信咨询。
1、优化目标
（1）成本最低
（2）污染最小
2、约束条件
功率平衡、水电站和火电机组约束，详见参考文献
3、代码
3.1 主函数 Main_hydro_thermal
clear clc global nHydro nTherm T HydroPowerGenrationCo Load HydroCons Inflow nVar % 数据来源：An interactive fuzzy satisfying method based on evolutionary % programming technique for multiobjective short-term hydrothermal scheduling HydroPowerGenrationCo = [ -0.0042, -0.42, 0.030, 0.90, 10.0, -50 ; -0.0040, -0.30, 0.015, 1.14, 9.5, -70 ; -0.0016, -0.30, 0.014, 0.55, 5.5, -40 ; -0.0030, -0.31, 0.027, 1.44, 14.0, -90 ]; % 水力发电系数 Load =[780,819,735,683,703,840,998,1060,1145,1134,1155,1200,1166,1081,1060,1113,1100,1170,1123.50000000000,1100,950,900,880,840]; % 库容最小值、库容最大值、初始库容、末时段库容、最小下泄、最大下泄、最小出力、最大出力 HydroCons = [ 80, 150, 100, 120, 5, 15, 0, 500 ; 60, 120, 80, 70, 6, 15, 0, 500 ; 100, 240, 170, 170, 10, 30, 0, 500 ; 70, 160, 120, 140, 6, 20, 0, 500 ]; Inflow = [10, 9, 8, 7, 6, 7, 8, 9, 10, 11, 12, 10, 11, 12, 11, 10, 9, 8, 7, 6, 7, 8, 9, 10; 8, 8, 9, 9, 8, 7, 6, 7, 8, 9, 9, 8, 8, 9, 9, 8, 7, 6, 7, 8, 9, 9, 8, 8; 8.1, 8.2, 4, 2, 3, 4, 3, 2, 1, 1, 1, 2, 4, 3, 3, 2, 2, 2, 1, 1, 2, 2, 1, 0 ; 2.8, 2.4, 1.6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ]; lb = [HydroCons(1,5)*ones(1,24),HydroCons(2,5)*ones(1,24),HydroCons(3,5)*ones(1,24), ... HydroCons(4,5)*ones(1,24),zeros(1,24),zeros(1,24),zeros(1,24),zeros(1,24),zeros(1,24),zeros(1,24),zeros(1,24),zeros(1,24),zeros(1,24),zeros(1,24)]; ub = [HydroCons(1,6)*ones(1,24),HydroCons(2,6)*ones(1,24),HydroCons(3,6)*ones(1,24), ... HydroCons(4,6)*ones(1,24),ones(1,24),ones(1,24),ones(1,24),ones(1,24),ones(1,24),ones(1,24),ones(1,24),ones(1,24),ones(1,24),ones(1,24)]; % 参数设置 T = 24; % 优化时长-1天24h nHydro = 4; % 梯级水电数量 nTherm = 10; % 火电机组数量 MaxIt = 40; % Maximum Number of Iterations nPop = 20; % Population Size pCrossover = 0.5; % Crossover Percentage nCrossover=30; pMutation = 0.01; % Mutation Percentage nMutation=20; mu = 0.01; % Mutation Rate sigma = 0.15*(ub-lb); % Mutation Step Size nDivision = 10; % 决策变量个数 nObj=2; % 目标函数个数 Zr = GenerateReferencePoints(nObj, nDivision); % 返回一个分区列表 params.nPop = nPop; params.Zr = Zr; params.nZr = size(Zr,2); params.zmin = []; params.zmax = []; params.smin = []; plotFlag = 1; nVar = (nHydro+nTherm)*T; % 决策变量数 [F]=NSGA_3(nPop,nVar,MaxIt,nObj,nCrossover,nMutation,mu,sigma,params);
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
- 38
- 39
- 40
- 41
- 42
- 43
- 44
- 45
- 46
- 47
- 48
3.2 函数-NSGA_3
function [F1]=NSGA_3(nPop,nVar,MaxIt,nObj,nCrossover,nMutation,mu,sigma,params) nDivision = 10; Zr = GenerateReferencePoints(nObj, nDivision); %% Initialization [pop,empty_individual] = Initial(nPop,nVar); [pop, F, params] = SortAndSelectPopulation(pop, params); %% NSGA-III Main Loop for it = 1:MaxIt it it % Crossover 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), popc(k,2)] = Crossover(p1, p2); end popc = popc(:); % Mutation popm = repmat(empty_individual, nMutation, 1); for k = 1:nMutation i = randi([1 nPop]); p = pop(i); popm(k,1) = Mutate(p, mu, sigma); end % Merge pop = [pop popc popm]; for k=1:size(pop,1) pop(k) = Reprocessing(pop(k)); pop(k) = Penalty(pop(k)); end [pop, F, params] = SortAndSelectPopulation(pop, params); F1 = pop(F{1}); disp(['Iteration ' num2str(it) ': Number of F1 Members = ' num2str(numel(F1))]); % 帕累托前沿可视化 figure(1); PlotCosts(F1); pause(0.01); end %% Results disp(['Final Iteration: Number of F1 Members = ' num2str(numel(F1))]); disp('Optimization Terminated.');
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
- 38
- 39
- 40
- 41
- 42
- 43
- 44
- 45
- 46
- 47
3.3 其他函数
function [pop, d, rho] = AssociateToReferencePoint(pop, params) Zr = params.Zr; nZr = params.nZr; rho = zeros(1,nZr); d = zeros(numel(pop), nZr); for i = 1:numel(pop) for j= 1:nZr w = Zr(:,j)/norm(Zr(:,j)); z = pop(i).NormalizedCost; d(i,j) = norm(z - w'*z*w); end [dmin, jmin] = min(d(i,:)); pop(i).AssociatedRef = jmin; pop(i).DistanceToAssociatedRef = dmin; rho(jmin) = rho(jmin) + 1; end end function [p1, p2]=Crossover(p1,p2) x1=p1.Position; x2=p2.Position; N=size(x1,2); alpha=rand(size(x1)); y1=alpha.*x1+(1-alpha).*x2; y2=alpha.*x2+(1-alpha).*x1; for i=97:N if y1(i)>0.5 y1(i)=1; else y1(i)=0; end if y2(i)>0.5 y2(i)=1; else y2(i)=0; end end p1.Position=y1; p2.Position=y2; end function b=Dominates(x,y) if isstruct(x) x=x.Cost; end if isstruct(y) y=y.Cost; end b=all(x<=y) && any(x<y); end function [fitness]=Fitness1(Pt,U,Toff) global nTherm T %目标函数 Tc=[5 5 4 4 4 2 2 0 0 0]; HST=[4500 5000 550 560 900 170 260 30 30 30]; CST=[9000 10000 1100 1120 1800 340 520 60 60 60]; Toffmin=[8 8 5 5 6 3 3 1 1 1]; InitialTOFF=[0,0,5,5,6,3,3,1,1,1]; coeff=[0.00048 16.19 1000; 0.00031 17.26 970; 0.00200 16.60 700; 0.00211 16.50 680; 0.00398 19.70 450; 0.00712 22.26 370; 0.00079 27.74 480; 0.00413 25.92 660; 0.00222 27.27 665 ; 0.00173 27.79 670]; fitness=0; C=zeros(nTherm,T); for k=1:nTherm for t=1:T if t==1 if InitialTOFF(k)<=(Tc(k)+Toffmin(k)) && InitialTOFF(k)>=Toffmin(k) C(k,t)=HST(k); elseif InitialTOFF(k)>(Tc(k)+Toffmin(k)) C(k,t)=CST(k); end elseif t>1 if Toff(k,t-1)<=(Tc(k)+Toffmin(k)) && Toff(k,t-1)>=Toffmin(k) C(k,t)=HST(k); elseif Toff(k,t-1)>(Tc(k)+Toffmin(k)) C(k,t)=CST(k); end end end end U0=[1 1 0 0 0 0 0 0 0 0]; for t=1:T p(t)=sum(Pt(:,t)); for k=1:10 if t==1 F=(coeff(k,:)*[Pt(k,t)*Pt(k,t) Pt(k,t) 1]'+C(k,t)*(1-U0(k)))*U(k,t); else F=(coeff(k,:)*[Pt(k,t)*Pt(k,t) Pt(k,t) 1]'+C(k,t)*(1-U(k,t-1)))*U(k,t); end fitness=fitness+F; end end function Fitness2=Fitness2(Pt,U) global nTherm T coeff=[130 -2.86 0.022; 132 -2.72 0.02; 137.7 -2.94 0.044; 130 -2.35 0.058; 125 -2.36 0.065; 110 -2.28 0.08; 135 -2.36 0.075; 157 -1.29 0.082; 160 -1.14 0.09; 137.7 -2.14 0.084]; Fitness2=0; for t=1:T for k=1:nTherm F=(coeff(k,:)*[1 Pt(k,t) Pt(k,t)*Pt(k,t)]')*U(k,t); Fitness2=Fitness2+F; end end function Zr = GenerateReferencePoints(M, p) Zr = GetFixedRowSumIntegerMatrix(M, p)' / p; end function A = GetFixedRowSumIntegerMatrix(M, RowSum) if M < 1 error('M cannot be less than 1.'); end if floor(M) ~= M % 向下取整 error('M must be an integer.'); end if M == 1 A = RowSum; return; end A = []; for i = 0:RowSum B = GetFixedRowSumIntegerMatrix(M - 1, RowSum - i); A = [A; i*ones(size(B,1),1) B]; end end function [pop] =Hydro(pop, plantID) global HydroCons T nVar deltaV = pop.HV(plantID,25) - HydroCons(plantID,4); counter = 0; while abs(deltaV)>=10 && (counter<10) avgV = deltaV/T; for i = 1:T pop.HQ(plantID,i) = pop.HQ(plantID,i) + avgV; if pop.HQ(plantID,i) < HydroCons(plantID,5) pop.HQ(plantID,i) = HydroCons(plantID,5); end if pop.HQ(plantID,i) > HydroCons(plantID,6) pop.HQ(plantID,i) = HydroCons(plantID,6); end end pop = V_Q(pop,plantID); deltaV = pop.HV(plantID,25) - HydroCons(plantID,4); counter = counter + 1; end pop = N_Q(pop, plantID); % 转换 pop.Position = reshape([pop.HQ' pop.TN'],1,nVar); function [Popp,Pop] = Initial(SearchAgents_no,dim) global HydroCons nHydro nTherm T %% NSGA-Ⅲ Initialization % 预留存储空间 Pop.Position = []; Pop.fitness = []; Pop.Best.Position=[]; Pop.Best.fitness=[]; Pop.Velocity = []; Pop.Dominated=false; Pop.IsDominated = []; Pop.GridIndex = []; Pop.GridSubIndex = []; Pop.HN = []; Pop.TN = []; Pop.HQ = []; Pop.HV = []; Pop.HVViolation = []; Pop.HPViolation = []; Pop.PowerViolation = []; Pop.TotalViolation = []; Pop.HEndVVio = []; Pop.constraint = []; Pop.Toff = []; Pop.Pt = []; Pop.Position = []; Pop.Cost = []; Pop.Rank = []; Pop.DominationSet = []; Pop.DominatedCount = []; Pop.NormalizedCost = []; Pop.AssociatedRef = []; Pop.DistanceToAssociatedRef = []; Pop.HNViolation=[]; Pop.HQViolation=[]; Popp = repmat(Pop,SearchAgents_no,1); % 初始化 for i = 1:SearchAgents_no Popp(i).Position = zeros(1,dim); Popp(i).Velocity = zeros(1,dim); Popp(i).HN = zeros(4,24); for k = 1:T Popp(i).HQ(:,k) = HydroCons(:,5) + rand(4,1).*(HydroCons(:,6) - HydroCons(:,5)); Popp(i).TN(:,k) = ones(10,1); end Popp(i).HV(:,1) = HydroCons(:,3); for j = 1:nHydro Popp(i) = V_Q(Popp(i),j); Popp(i) = N_Q(Popp(i),j); Popp(i) = Hydro(Popp(i),j); end Popp(i) = Thermal(Popp(i)); Popp(i).Position = reshape([Popp(i).HQ' Popp(i).TN'],1,dim); Popp(i) = Penalty(Popp(i)); end function [Pt]=Loaddispatch_thermal(U,Pload) global nTherm T h = [1,2,3,4,5,6,7,8,9,10]; % 火电机组出力顺序 coeff=[0.00048 16.19 1000; 0.00031 17.26 970; 0.00200 16.60 700; 0.00211 16.50 680; 0.00398 19.70 450; 0.00712 22.26 370; 0.00079 27.74 480; 0.00413 25.92 660; 0.00222 27.27 665 ; 0.00173 27.79 670]; n=size(coeff,1); c2=coeff(:,1); c1=coeff(:,2); % P_ramp=[130 130 60 60 90 40 40 40 40 40]; Pt=zeros(nTherm,T); lambda=10*ones(1,T); lambda_max=200*ones(1,T); lambda_min=0*zeros(1,T); epsilon=0.1; lostpower=epsilon*ones(1,T); Pmax=[455 455 130 130 162 80 85 55 55 55]; Pmin=[150 150 20 20 25 20 25 10 10 10]; P_max=zeros(1,nTherm); P_min=zeros(1,nTherm); for t=1:T maxpower=0; for k =1:n maxpower=maxpower+U(k,t)*Pmax(k); end if Pload(t)<=maxpower iter=1; while abs(lostpower(t))>=epsilon for k=1:n if U(h(k),t)==1 Pt(h(k),t)=(lambda(t)-c1(h(k)))/(3*c2(h(k))); if Pt(h(k),t)>Pmax(h(k)) Pt(h(k),t)=Pmax(h(k)); elseif Pt(h(k),t)<Pmin(h(k)) Pt(h(k),t)=Pmin(h(k)); end end end sumpower=0; for k=1:n sumpower=sumpower+U(k,t)*Pt(k,t); end lostpower(t)=Pload(t)-sumpower; if lostpower(t)>0 lambda_max(t)=lambda_max(t); lambda_min(t)=lambda(t); lambda(t)=(lambda_max(t)+lambda_min(t))/2; else lambda_max(t)=lambda(t); lambda_min(t)=lambda_min(t); lambda(t)=(lambda_max(t)+lambda_min(t))/2; end iter=iter+1; if iter==10 temPt=0*zeros(1,nTherm)'; Pt=repmat(temPt,1,T); return end end end end function z=MOP2(x) n=numel(x); z1=1-exp(-sum((x-1/sqrt(n)).^2)); z2=1-exp(-sum((x+1/sqrt(n)).^2)); z=[z1 z2]'; end function p=Mutate(p,mu,sigma) x=p.Position; nVar=numel(x); nMu=ceil(mu*nVar); j=randsample(nVar,nMu); %产生nMu个不同的数，数字的范围为[1，nVar] y=x; if j<97 y(j)=x(j)+sigma(j)*randn(size(j)); else y(j)=~y(j); end p.Position=y; end function pop = N_Q(pop,id) % 计算水电出力 global HydroPowerGenrationCo T pop.HN(id,:) = HydroPowerGenrationCo(id,1)*pop.HV(id,1:T).^2 + HydroPowerGenrationCo(id,2)*pop.HQ(id,:).^2 + ... HydroPowerGenrationCo(id,3)*pop.HV(id,1:T).*pop.HQ(id,:) + HydroPowerGenrationCo(id,4)*pop.HV(id,1:T) + ... HydroPowerGenrationCo(id,5)*pop.HQ(id,:) + HydroPowerGenrationCo(id,6); function [pop, F]=NonDominatedSorting(pop) nPop=numel(pop); for i=1:nPop pop(i).DominationSet=[]; pop(i).DominatedCount=0; end F{1}=[]; for i=1:nPop for j=i+1:nPop p=pop(i); q=pop(j); if Dominates(p,q) p.DominationSet=[p.DominationSet j]; q.DominatedCount=q.DominatedCount+1; end if Dominates(q.Cost,p.Cost) q.DominationSet=[q.DominationSet i]; p.DominatedCount=p.DominatedCount+1; end pop(i)=p; pop(j)=q; end if pop(i).DominatedCount==0 F{1}=[F{1} i]; pop(i).Rank=1; end end k=1; while true Q=[]; for i=F{k} p=pop(i); for j=p.DominationSet q=pop(j); q.DominatedCount=q.DominatedCount-1; if q.DominatedCount==0 Q=[Q j]; q.Rank=k+1; end pop(j)=q; end end if isempty(Q) break; end F{k+1}=Q; k=k+1; end end function [pop, params] = NormalizePopulation(pop, params) params.zmin = UpdateIdealPoint(pop, params.zmin); fp = [pop.Cost] - repmat(params.zmin, 1, numel(pop)); params = PerformScalarizing(fp, params); a = FindHyperplaneIntercepts(params.zmax); for i = 1:numel(pop) pop(i).NormalizedCost = fp(:,i)./a; end end function a = FindHyperplaneIntercepts(zmax) w = ones(1, size(zmax,2))/zmax; a = (1./w)'; end function [pop] = Penalty(pop) % 惩罚函数及适应度值——超约束时可调参 global Load HydroCons nHydro T %1 负荷平衡 pop.PowerViolation = sum( abs( sum(pop.HN) + sum(pop.Pt) - Load ) ); %2 水位约束 tempHVvio = 0; for i = 1:T temp = 0; for j = 1:nHydro if pop.HV(j,i) < HydroCons(j,1) temp = temp + abs(pop.HV(j,i) - HydroCons(j,1)); elseif pop.HV(j,i) > HydroCons(j,2) temp = temp + abs(pop.HV(j,i) - HydroCons(j,2)); end end tempHVvio = tempHVvio + temp; end pop.HVViolation = tempHVvio; %3 末水位约束 pop.HEndVVio = sum( abs( pop.HV(:,25) - HydroCons(:,4) ) ); %4 下泄约束 tempHQvio = 0; for i = 1:T temp = 0; for j = 1:nHydro if pop.HQ(j,i) < HydroCons(j,5) temp = temp + abs(pop.HQ(j,i) - HydroCons(j,5)); elseif pop.HQ(j,i) > HydroCons(j,6) temp = temp + abs(pop.HQ(j,i) - HydroCons(j,6)); end end tempHQvio = tempHQvio + temp; end pop.HQViolation = tempHQvio; %4 出力约束 tempHNvio = 0; for i = 1:nHydro for j = 1:T if pop.HN(i,j) < HydroCons(i,7) tempHNvio = tempHNvio + abs(pop.HN(i,j) - HydroCons(i,7)); elseif pop.HN(i,j) > HydroCons(i,8) tempHNvio = tempHNvio + abs(pop.HN(i,j) - HydroCons(i,8)); end end end pop.HNViolation = tempHNvio; % 总惩罚 pop.TotalViolation = 1000*pop.HVViolation +1000*pop.HQViolation + 1000*pop.HEndVVio + 1000*pop.PowerViolation + 1000*pop.HNViolation; %% 计算目标函数 pop.Cost(1,1) = Fitness1(pop.Pt, pop.TN, pop.Toff)+ pop.TotalViolation; pop.Cost(2,1) = Fitness2(pop.Pt, pop.TN)+ pop.TotalViolation; function params = PerformScalarizing(z, params) nObj = size(z,1); nPop = size(z,2); if ~isempty(params.smin) zmax = params.zmax; smin = params.smin; else zmax = zeros(nObj, nObj); smin = inf(1,nObj); end for j = 1:nObj w = GetScalarizingVector(nObj, j); s = zeros(1,nPop); for i = 1:nPop s(i) = max(z(:,i)./w); end [sminj, ind] = min(s); if sminj < smin(j) zmax(:, j) = z(:, ind); smin(j) = sminj; end end params.zmax = zmax; params.smin = smin; end function w = GetScalarizingVector(nObj, j) epsilon = 1e-10; w = epsilon*ones(nObj, 1); w(j) = 1; end function PlotCosts(pop) Costs=[pop.Cost]; plot(Costs(1,:),Costs(2,:),'r*','MarkerSize',8); xlabel('1st Objective'); ylabel('2nd Objective'); grid on; end function pop = Reprocessing(pop) global HydroCons nVar nHydro nTherm T for i = 1:nHydro pop.HQ(i,:) = pop.Position(T*(i-1) + 1 : T*i); end for i = 1:nTherm pop.TN(i,:) = pop.Position(T*4 + 1 + T*(i - 1) : T*4 + T*i); end pop.HV(:,1) = HydroCons(:,3); for i = 1:nHydro pop = V_Q(pop,i); pop = N_Q(pop,i); pop = Hydro(pop,i); end pop = Thermal(pop); pop.Position = reshape([pop.HQ' pop.TN'],1,nVar); function [pop, F, params] = SortAndSelectPopulation(pop, params) [pop, params] = NormalizePopulation(pop, params); [pop, F] = NonDominatedSorting(pop); nPop = params.nPop; if numel(pop) == nPop return; end [pop, d, rho] = AssociateToReferencePoint(pop, params); newpop = []; for l=1:numel(F) if numel(newpop) + numel(F{l}) > nPop LastFront = F{l}; break; end newpop = [newpop; pop(F{l})]; %#ok end while true [~, j] = min(rho); AssocitedFromLastFront = []; for i = LastFront if pop(i).AssociatedRef == j AssocitedFromLastFront = [AssocitedFromLastFront i]; %#ok end end if isempty(AssocitedFromLastFront) rho(j) = inf; continue; end if rho(j) == 0 ddj = d(AssocitedFromLastFront, j); [~, new_member_ind] = min(ddj); else new_member_ind = randi(numel(AssocitedFromLastFront)); end MemberToAdd = AssocitedFromLastFront(new_member_ind); LastFront(LastFront == MemberToAdd) = []; newpop = [newpop; pop(MemberToAdd)]; %#ok rho(j) = rho(j) + 1; if numel(newpop) >= nPop break; end end [pop, F] = NonDominatedSorting(newpop); end function pop = Thermal(pop) global Load nVar ThermalLoad = Load - sum(pop.HN) ; % 火电机组所承担负荷 [pop.TN, pop.Toff] = Unitcommitment_thermal(pop.TN, ThermalLoad); [pop.Pt]=Loaddispatch_thermal(pop.TN,ThermalLoad); pop.Position = reshape([pop.HQ' pop.TN'],1,nVar); function [U, TOFF]=Unitcommitment_thermal(U, PD) global nTherm T % 决定火电机组启停 G=nTherm; TON=zeros(G,T); TOFF=zeros(G,T); SR=0.05*PD; PMAX=[455,455,130,130,162,80,85,55,55,55]; PMIN=[150 150 20 20 25 20 25 10 10 10]; h = [1,2,3,4,5,6,7,8,9,10]; % 火电机组出力顺序 InitialTON=[8,8,0,0,0,0,0,0,0,0]; InitialTOFF=[0,0,5,5,6,3,3,1,1,1]; MDT=[8,8,5,5,6,3,3,1,1,1]; MUT=[8,8,5,5,6,3,3,1,1,1]; for t=1:T if t==1 TON(:,t)=(InitialTON.*U(G*(t-1)+1:G*t)+U(G*(t-1)+1:G*t))'; TOFF(:,t)=InitialTOFF'.*~TON(:,t)+~TON(:,t); else flag=(TOFF(:,t-1)==0&TON(:,t-1)<MUT')|(TON(:,t-1)==0&TOFF(:,t-1)<MDT'); U(G*(t-1)+1:G*t)=U(G*(t-2)+1:G*(t-1)).*flag'+U(G*(t-1)+1:G*t).*~flag'; TON(:,t)=TON(:,t-1).*U(G*(t-1)+1:G*t)'+U(G*(t-1)+1:G*t)'; TOFF(:,t)=TOFF(:,t-1).*~TON(:,t)+~TON(:,t); end P=PMAX*U(G*(t-1)+1:G*t)'; if P<PD(t)+SR(t) for j=1:G g=h(j); if U(g+G*(t-1))==1 continue; else U(g+G*(t-1))=1; if TOFF(g,t)>MDT(g) if t==1 TON(g,t)=InitialTON(g)+1; else TON(g,t)=TON(g,t-1)+1; end TOFF(g,t)=0; else l=t-TOFF(g,t)+1; if l<=0 l=1; end U(g+G*(l-1:t-1))=1; if l==1 TON(g,l)=InitialTON(g)+1; TON(g,l+1:t)=TON(g,l)+[1:t-l]; else TON(g,l:t)=TON(g,l-1)+[1:t-l+1]; end TOFF(g,l:t)=0; end P=PMAX*U(G*(t-1)+1:G*t)'; if P>=PD(t)+SR(t) break; end end end end for i=1:G g=h(G+1-i); P1=PMAX*U(G*(t-1)+1:G*t)'; P2=PMIN*U(G*(t-1)+1:G*t)'; if U(g+G*(t-1))==1 if P1-PMAX(g)>=PD(t)+SR(t) if TON(g,t)>MUT(g) || TON(g,t)==1 U(g+G*(t-1))=0; TON(g,t)=0; if t==1 TOFF(g,t)=InitialTOFF(g)+1; else TOFF(g,t)=TOFF(g,t-1)+1; end else continue; end else break; end end end end function zmin = UpdateIdealPoint(pop, prev_zmin) if ~exist('prev_zmin', 'var') || isempty(prev_zmin) prev_zmin = inf(size(pop(1).Cost)); end zmin = prev_zmin; N=numel(pop); for i = 1:numel(pop) zmin = min(zmin, pop(i).Cost); end end function pop = V_Q(pop,id) % 计算库容 global Inflow HydroCons if id == 1 || id == 2 for i = 1:24 preV = pop.HV(id,i); nextV = preV + Inflow(id,i) - pop.HQ(id,i); if nextV < HydroCons(id,1) tempQ = preV + Inflow(id,i) - HydroCons(id,1); if tempQ >= HydroCons(id,5) pop.HQ(id,i) = tempQ; nextV = HydroCons(id,1); end elseif nextV > HydroCons(id,2) tempQ = preV + Inflow(id,i) - HydroCons(id,2); if tempQ <= HydroCons(id,6) pop.HQ(id,i) = tempQ; nextV = HydroCons(id,2); end end pop.HV(id,i+1) = nextV; end elseif id == 3 || id == 4 for i = 1:24 preV = pop.HV(id,i); if id == 3 if i < 3 nextV = preV + Inflow(id,i) - pop.HQ(id,i); elseif i < 4 nextV = preV + Inflow(id,i) + pop.HQ(1,i - 2) - pop.HQ(id,i); else nextV = preV + Inflow(id,i) + pop.HQ(1,i - 2) + pop.HQ(2,i - 3) - pop.HQ(id,i); end else if i < 5 nextV = preV + Inflow(id,i) - pop.HQ(id,i); else nextV = preV + Inflow(id,i) + pop.HQ(3,i - 4) - pop.HQ(id,i); end end if nextV < HydroCons(id,1) tempQ = pop.HQ(id,i) - (HydroCons(id,1) - nextV); if tempQ >= HydroCons(id,5) pop.HQ(id,i) = tempQ; nextV = HydroCons(id,1); end elseif nextV > HydroCons(id,2) tempQ = pop.HQ(id,i) + (nextV - HydroCons(id,2)); if tempQ <= HydroCons(id,6) pop.HQ(id,i) = tempQ; nextV = HydroCons(id,2); end end pop.HV(id,i + 1) = nextV; end end
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
- 38
- 39
- 40
- 41
- 42
- 43
- 44
- 45
- 46
- 47
- 48
- 49
- 50
- 51
- 52
- 53
- 54
- 55
- 56
- 57
- 58
- 59
- 60
- 61
- 62
- 63
- 64
- 65
- 66
- 67
- 68
- 69
- 70
- 71
- 72
- 73
- 74
- 75
- 76
- 77
- 78
- 79
- 80
- 81
- 82
- 83
- 84
- 85
- 86
- 87
- 88
- 89
- 90
- 91
- 92
- 93
- 94
- 95
- 96
- 97
- 98
- 99
- 100
- 101
- 102
- 103
- 104
- 105
- 106
- 107
- 108
- 109
- 110
- 111
- 112
- 113
- 114
- 115
- 116
- 117
- 118
- 119
- 120
- 121
- 122
- 123
- 124
- 125
- 126
- 127
- 128
- 129
- 130
- 131
- 132
- 133
- 134
- 135
- 136
- 137
- 138
- 139
- 140
- 141
- 142
- 143
- 144
- 145
- 146
- 147
- 148
- 149
- 150
- 151
- 152
- 153
- 154
- 155
- 156
- 157
- 158
- 159
- 160
- 161
- 162
- 163
- 164
- 165
- 166
- 167
- 168
- 169
- 170
- 171
- 172
- 173
- 174
- 175
- 176
- 177
- 178
- 179
- 180
- 181
- 182
- 183
- 184
- 185
- 186
- 187
- 188
- 189
- 190
- 191
- 192
- 193
- 194
- 195
- 196
- 197
- 198
- 199
- 200
- 201
- 202
- 203
- 204
- 205
- 206
- 207
- 208
- 209
- 210
- 211
- 212
- 213
- 214
- 215
- 216
- 217
- 218
- 219
- 220
- 221
- 222
- 223
- 224
- 225
- 226
- 227
- 228
- 229
- 230
- 231
- 232
- 233
- 234
- 235
- 236
- 237
- 238
- 239
- 240
- 241
- 242
- 243
- 244
- 245
- 246
- 247
- 248
- 249
- 250
- 251
- 252
- 253
- 254
- 255
- 256
- 257
- 258
- 259
- 260
- 261
- 262
- 263
- 264
- 265
- 266
- 267
- 268
- 269
- 270
- 271
- 272
- 273
- 274
- 275
- 276
- 277
- 278
- 279
- 280
- 281
- 282
- 283
- 284
- 285
- 286
- 287
- 288
- 289
- 290
- 291
- 292
- 293
- 294
- 295
- 296
- 297
- 298
- 299
- 300
- 301
- 302
- 303
- 304
- 305
- 306
- 307
- 308
- 309
- 310
- 311
- 312
- 313
- 314
- 315
- 316
- 317
- 318
- 319
- 320
- 321
- 322
- 323
- 324
- 325
- 326
- 327
- 328
- 329
- 330
- 331
- 332
- 333
- 334
- 335
- 336
- 337
- 338
- 339
- 340
- 341
- 342
- 343
- 344
- 345
- 346
- 347
- 348
- 349
- 350
- 351
- 352
- 353
- 354
- 355
- 356
- 357
- 358
- 359
- 360
- 361
- 362
- 363
- 364
- 365
- 366
- 367
- 368
- 369
- 370
- 371
- 372
- 373
- 374
- 375
- 376
- 377
- 378
- 379
- 380
- 381
- 382
- 383
- 384
- 385
- 386
- 387
- 388
- 389
- 390
- 391
- 392
- 393
- 394
- 395
- 396
- 397
- 398
- 399
- 400
- 401
- 402
- 403
- 404
- 405
- 406
- 407
- 408
- 409
- 410
- 411
- 412
- 413
- 414
- 415
- 416
- 417
- 418
- 419
- 420
- 421
- 422
- 423
- 424
- 425
- 426
- 427
- 428
- 429
- 430
- 431
- 432
- 433
- 434
- 435
- 436
- 437
- 438
- 439
- 440
- 441
- 442
- 443
- 444
- 445
- 446
- 447
- 448
- 449
- 450
- 451
- 452
- 453
- 454
- 455
- 456
- 457
- 458
- 459
- 460
- 461
- 462
- 463
- 464
- 465
- 466
- 467
- 468
- 469
- 470
- 471
- 472
- 473
- 474
- 475
- 476
- 477
- 478
- 479
- 480
- 481
- 482
- 483
- 484
- 485
- 486
- 487
- 488
- 489
- 490
- 491
- 492
- 493
- 494
- 495
- 496
- 497
- 498
- 499
- 500
- 501
- 502
- 503
- 504
- 505
- 506
- 507
- 508
- 509
- 510
- 511
- 512
- 513
- 514
- 515
- 516
- 517
- 518
- 519
- 520
- 521
- 522
- 523
- 524
- 525
- 526
- 527
- 528
- 529
- 530
- 531
- 532
- 533
- 534
- 535
- 536
- 537
- 538
- 539
- 540
- 541
- 542
- 543
- 544
- 545
- 546
- 547
- 548
- 549
- 550
- 551
- 552
- 553
- 554
- 555
- 556
- 557
- 558
- 559
- 560
- 561
- 562
- 563
- 564
- 565
- 566
- 567
- 568
- 569
- 570
- 571
- 572
- 573
- 574
- 575
- 576
- 577
- 578
- 579
- 580
- 581
- 582
- 583
- 584
- 585
- 586
- 587
- 588
- 589
- 590
- 591
- 592
- 593
- 594
- 595
- 596
- 597
- 598
- 599
- 600
- 601
- 602
- 603
- 604
- 605
- 606
- 607
- 608
- 609
- 610
- 611
- 612
- 613
- 614
- 615
- 616
- 617
- 618
- 619
- 620
- 621
- 622
- 623
- 624
- 625
- 626
- 627
- 628
- 629
- 630
- 631
- 632
- 633
- 634
- 635
- 636
- 637
- 638
- 639
- 640
- 641
- 642
- 643
- 644
- 645
- 646
- 647
- 648
- 649
- 650
- 651
- 652
- 653
- 654
- 655
- 656
- 657
- 658
- 659
- 660
- 661
- 662
- 663
- 664
- 665
- 666
- 667
- 668
- 669
- 670
- 671
- 672
- 673
- 674
- 675
- 676
- 677
- 678
- 679
- 680
- 681
- 682
- 683
- 684
- 685
- 686
- 687
- 688
- 689
- 690
- 691
- 692
- 693
- 694
- 695
- 696
- 697
- 698
- 699
- 700
- 701
- 702
- 703
- 704
- 705
- 706
- 707
- 708
- 709
- 710
- 711
- 712
- 713
- 714
- 715
- 716
- 717
- 718
- 719
- 720
- 721
- 722
- 723
- 724
- 725
- 726
- 727
- 728
- 729
- 730
- 731
- 732
- 733
- 734
- 735
- 736
- 737
- 738
- 739
- 740
- 741
- 742
- 743
- 744
- 745
- 746
- 747
- 748
- 749
- 750
- 751
- 752
- 753
- 754
- 755
- 756
- 757
- 758
- 759
- 760
- 761
- 762
- 763
- 764
- 765
- 766
- 767
- 768
- 769
- 770
- 771
- 772
- 773
- 774
- 775
- 776
- 777
- 778
- 779
- 780
- 781
- 782
参考
本文内容由网友自发贡献，转载请注明出处：【wpsshop博客】
推荐阅读
相关标签
