分时电价模型，削峰填谷，转移24小时一天中用电率，减少谷峰差

close all; clear; clc; loading = [5110 5325 5200 4910 5435 5660 5960 6435 7175 7340 7525 6700 6825 ... 6925 6880 6860 7290 7350 7780 7750 7290 6810 5815 5535]; %% 其他数据 I = 3; J = 4; Lmin = min(loading); Lmax = max(loading); Pp = 0.42; Afg = 0.15; Afp = 0.15; Apg = 0.15; Kfg = 0.25; Kfp = 0.25; Kpg = 0.25; LambdaFG_MAX = 0.08; LambdaFP_MAX = 0.08; LambdaPG_MAX = 0.08; %j种谷峰电价一共有三种，分为平时电价假设为1，峰时电价假设为平时的1.5倍，谷时电价为平时的0.8倍 % 时间划分 F1f = [7 11;17 21]; F1p = [11 17;21 23]; F1g = [23 7]; F2f = [9 11;17 21]; F2p = [7 9;11 17;21 23]; F2g = [23 9]; F3f = [8 12;18 22]; F3p = [12 18;22 24]; F3g = [0 8]; Ff = {F1f F2f F3f}; Fp = {F1p F2p F3p}; Fg = {F1g F2g F3g}; K = [2.5 1.5;2.8 1.37;3.0 1.6;3.2 1.45]; kX = size(K,1); loading_result = {}; index = 1; for j = 1 : 3 %3F for i = 1 : kX %4K fengK = K(i,2)*Pp; guK = K(i,2)/K(i,1)*Pp; pingK = Pp; deltaFG = fengK - guK; deltaFP = fengK - pingK; deltaPG = pingK - guK; lambdaFG = getLambdaFG(deltaFG,Afg,Kfg,LambdaFG_MAX); lambdaFP = getLambdaFP(deltaFP,Afp,Kfp,LambdaFP_MAX); lambdaPG = getLambdaPG(deltaPG,Apg,Kpg,LambdaPG_MAX); if( j == 3 && i == kX) j i end simData = loading; %峰 currentFTime = cell2mat(Ff(j)); lengthOfCurrentFTime = size(currentFTime,1); for w = 1 : lengthOfCurrentFTime for k = currentFTime(w,1)+1:currentFTime(w,2) simData(k) = (1 - lambdaFG - lambdaFP) * simData(k); end end %平 currentPTime = cell2mat(Fp(j)); lengthOfCurrentPTime = size(currentPTime,1); for w = 1 : lengthOfCurrentPTime for k = currentPTime(w,1)+1:currentPTime(w,2) simData(k) = (1 + lambdaFP - lambdaPG) * simData(k); end end %谷 currentGTime = cell2mat(Fg(j)); lengthOfCurrentGTime = size(currentGTime,1); for w = 1 : lengthOfCurrentGTime if(currentGTime(w,1)+1 < 24) for k = currentGTime(w,1)+1 : 24 simData(k) = (1 + lambdaFG + lambdaPG) * simData(k); end for k = 1 : currentGTime(w,2) simData(k) = (1 + lambdaFG + lambdaPG) * simData(k); end elseif (currentGTime(w,1)+1 == 24) simData(24) = (1 + lambdaFG + lambdaPG) * simData(24); for k = 1:currentGTime(w,2) simData(k) = (1 + lambdaFG + lambdaPG) * simData(k); end else for k = currentGTime(w,1)+1:currentGTime(w,2) simData(k) = (1 + lambdaFG + lambdaPG) * simData(k); end end end loading_result(index) = {simData}; index = index + 1; end end DD = []; for i = 1: 3*kX curLoading = loading_result{i}; cMax = max(curLoading); cMin = min(curLoading); DD = [DD;cMax cMin cMax-cMin]; end % Find P2 = find(DD(:,1) == min(DD(:,1))); P3 = find(DD(:,3) == min(DD(:,3))); OPT = intersect(P2,P3); plot(loading,'-r'); for i = 1:size(OPT) hold on; plot(loading_result{OPT(i)},'--b'); end legend({'Original Loading Curve','Optimized'},'Location','NorthEast'); xlabel({'Time(Hour)'},'FontSize',11); ylabel({'Load(MW)'},'FontSize',11); title('Typical daily load curve','FontWeight','bold','FontSize',14); axis([0 24 0 10000]); set(gca,'XTick',0:4:24); grid on;