电力负荷分配part3

%function [leastCostPt U]=UC（ coeff Pmin Pmax Ploadlambda）

% leastCost为机组的最小发电成本，是某一时间段内所有机组的发电成本=========>>经济调度
% Pt为机组在各时段的发电功率，行对应一个机组，列对应某一时间段。最终需要的是Pt而不是P因为Pt才是经济调度的结果，满足负荷平衡约束
% U为机组在各时段的运行状态，行对应一个机组，列对应某一时间段。若其值为1，则表示机组投运；若值为0，则表示机组停机

% coeff为各机组的成本系数矩阵，设为二次函数形式[  二次项系数   一次项系数    常数项系数    ]
% PminPmax为各机组的发电功率下限和上限向量（列向量）
% Pload为各时段的负荷向量（行向量）
% lambda为拉格朗日系数初值

coeff=[0.00128  6.48  459;
       0.00194  7.85  310;
       0.00482  7.97  78]; % 每一行对应一个机组
Pmin=[150 ;100 ;50];
Pmax=[600 ;400 ;200];
Pload=[850];
n=size（coeff1）;% n 为机组数
T=size（Pload2）;% T 为需要计算的总时间段数
lambda=zeros（1T）;
      
%        num=nargin;
%        if（ num==6）
%              lambda=zeros（1T）;
%        end
       
       
       q=0;                 % 对偶问题
       minCost=zeros（n1）;  %对应每个机组在所有时间段内的发电成本=========>>动态规划
       P=zeros（nT）;        %对应每个机组=============>>动态规划
       U=zeros（nT）;        % 每个机组在所有时间段内的运行状态，若其值为1，则表示机组投运；若值为0，则表示机组停机
       lostPower=zeros（1T）;% 各时段内的系统发电功率缺额=========>>动态规划结束后
       
       J=0;                 % 原问题
       leastCost=zeros（1T）;%对应每个时间段内所有机组的发电成本========>>经济调度
       Pt=zeros（nT）;       %对应每个时间段==========>>经济调度
       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）Tlambda）;% 利用动态规划求解每个机组在所有时段的运行状态和最优发电功率
          end
          %=========================================================================>>
          % 本部分用于求对偶问题最优值q
          sumPower=zeros（1T）;% 系统各时段的总发电功率
          for t=1:T
              for i=1:n
                  sumPower（t）=sumPower（t） + U（it）*P（it）;% 求各时段内的发电功率
               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（it）==1   %  判断某时段内是否有机组投运
                             flag=1;
                         end    
                    end  
                    if flag==1   % 如果所有机组停运则不必求解
                         [leastCost（t）Pt（:t）] = ecoDisp