%function [leastCost,Pt ,U]=UC( coeff, Pmin, Pmax, Pload,lambda)
% leastCost为机组的最小发电成本,是某一时间段内所有机组的发电成本=========>>经济调度
% Pt为机组在各时段的发电功率,行对应一个机组,列对应某一时间段。最终需要的是Pt,而不是P,因为Pt才是经济调度的结果,满足负荷平衡约束
% U为机组在各时段的运行状态,行对应一个机组,列对应某一时间段。若其值为1,则表示机组投运;若值为0,则表示机组停机
% coeff为各机组的成本系数矩阵,设为二次函数形式[ 二次项系数 一次项系数 常数项系数 ]
% Pmin,Pmax为各机组的发电功率下限和上限向量(列向量)
% Pload为各时段的负荷向量(行向量)
% lambda为拉格朗日系数初值
coeff=[0.002 10 500;
0.0025 8 300;
0.005 6 100]; % 每一行对应一个机组
Pmin=[100 ;100 ;50];
Pmax=[600 ;400 ;200];
Pload=[170 520 1100 330];
n=size(coeff,1);% n 为机组数
T=size(Pload,2);% T 为需要计算的总时间段数
lambda=zeros(1,T);
% num=nargin;
% if( num==6)
% lambda=zeros(1,T);
% end
q=0; % 对偶问题
minCost=zeros(n,1); %对应每个机组在所有时间段内的发电成本=========>>动态规划
P=zeros(n,T); %对应每个机组=============>>动态规划
U=zeros(n,T); % 每个机组在所有时间段内的运行状态,若其值为1,则表示机组投运;若值为0,则表示机组停机
lostPower=zeros(1,T);% 各时段内的系统发电功率缺额=========>>动态规划结束后
J=0; % 原问题
leastCost=zeros(1,T);%对应每个时间段内所有机组的发电成本========>>经济调度
Pt=zeros(n,T); %对应每个时间段==========>>经济调度
lackCost=10000; % 当机组容量不能满足经济调度时,强制设定此时的总发电成本为10000
flag=0;
iter=1; % 计算迭代次数
epsilon=0.05; % 设定容许误差
gap=0.05; % 设定对偶间隙初值
while( abs(gap)>=epsilon )
for i=1:n % n为机组数
[ minCost(i),P(i,:),U(i,:) ] = dp( coeff(i,:),Pmin(i),Pmax(i),T,lambda);% 利用动态规划求解每个机组在所有时段的运行状态和最优发电功率
end
%=========================================================================>>
% 本部分用于求对偶问题最优值q
sumPower=zeros(1,T);% 系统各时段的总发电功率
for t=1:T
for i=1:n
sumPower(t)=sumPower(t) + U(i,t)*P(i,t);% 求各时段内的发电功率
end
end
lostPower=Pload-sumPower;% 求得各时段的系统发电功率缺额
q=0;
for i=1:n
q = q + minCost(i); % 叠加各机组的最小发电成本
end
for t=1:T
q = q + lambda(t)*lostPower(t); % 叠加拉格朗日乘子项
end
%=========================================================================>>
%本部分用于求原问题最优值J
for t=1:T
for i=1:n
if U(i,t)==1 % 判断某时段内是否有机组投运
flag=1;
end
end
if flag==1 % 如果所有机组停运则不必求解
[leastCost(t),Pt(:,t)] = ecoDisp( coeff,Pmin,Pmax,U(:,t),Pload(t) );% 对每个时段进行经济调度计算
end
flag=0;% 恢复初始状态
end
J=0;
for t=1:T
if( sum(Pt(:,t))>0) % 判断是否实现了经济调度
J=J + leastCost(t);
else
J=J + lackCost;% 当机组容量不能满足经济调度时,强制设定此时的总发电成本为10000
end
end
%=========================================================================>>
%本部分用于求相对对偶间隙
if q~=0
gap=(J-q)/q;
end
%=========================================================================>>
%本部分输出迭代结果
fprintf('当前的迭代次数是 %4d\n',iter);
fprintf('Hour lambda ');
for i=1:n
fprintf('u%d ',i);
end
for i=1:n
fprintf(' P%d ',i);
end
fprintf('Pload(t)-sum(Pi*Ui) ');
for i=1:n
fprintf(' Pt%d ',i);
end
fprintf('\n');
for t=1:T
fprintf('%d %8.4f ',t,lambda(t));
for i=1:n
fprintf('%d ',U(i,t));
end
for i=1:n
fprintf('%4d ',P(i,t));
end
fprintf(' %8.2f ',lostPower(t));
for i=1:n
fprintf(' %10.4f',Pt(i,t));
end
fprintf('\n');
end
fprintf('\n q*= %f , J*= %f ',q,J);
if q~=0
fprintf('dual_gap = (J*-q*)/q* = %f \n\n\n\n',gap);
else
fprintf('dual_gap = (J*-q*)/q* = undefined \n\n\n\n');
end
%=========================================================================>>
if ( abs(gap)>=epsilon ) % 判断是否有必要为下一步迭代更新数据
leastCost=zeros(1,T);% 恢复初始状态
Pt=zeros(n,T);% 恢复初始状态
%===============================>>
%本部分用于修正拉格朗日乘子lambda
for t=1:T % 计算步长因子
if lostPower(t)>0
alpha=0.01;
elseif lostPower(t)<0
alpha=0.002;
end
lambda(t)=lambda(t)+lostPower(t)*alpha; % 修正拉格朗日乘子lambda
end
iter=iter+1;
end
%=========================================================================>>
end % while循环终止处