论文复现-基于NSGA-Ⅲ的梯级水电-火电机组多目标优化调度
火电是目前我国电网中提供调节能力的主要能源，而随着碳中和目标的提出，势必要降低火电在电网中的分量。梯级水电则是可以替代火电发挥调节作用的能源。基于电网中存在大量火电的背景，梯级水电如何与火电机组耦合调度呢？本代码通过优于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);
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.');
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
142
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
