tic;
unit=[1515;%電容量
????25.525.5;%熱容量??%%未考慮機組投運順序問題
????0.00220.0022];%3*m數組，各個機組電容量熱容量與故障率
pzzxs=1;%園區電負荷增長系數
hzzxs=1;%園區熱負荷增長系數
load（‘chpdata.mat‘）;%讀入歷史電負荷數據
load（‘chhdata.mat‘）;%讀入歷史熱負荷數據
load（‘pvdata.mat‘）;%讀入歷史光伏出力數據
load（‘ps.mat‘）;?%讀入從電網購電分時電價
load（‘priceuser.mat‘）;?%讀入運營商向用戶售電分時電價
load（‘day.mat‘）;
%load（‘hql.mat‘）;
pv=pvdata‘;
%根據歷史電負荷，考慮一定的負荷增長速度，得到未來一段時間內電負荷
p=pzzxs*chpdata‘;
%熱負荷根據歷史熱負荷，考慮一定的負荷增長速度，得到未來一段時間內熱負荷
h=hzzxs*chhdata‘;
%將光伏出力從電負荷曲線中分離，得到凈電負荷
for?s=1:8760?
????if?p（1s）>=pv（1s）
????????p（1s）=p（1s）-pv（1s）;
????else
????????p（1s）=0;
????end
end
prl=unit（1:）;%求各個機組容量數組
hrl=unit（2:）;%求各個機組容量數組
gz=unit（3:）;%求各個機組故障數組
jizu=length（prl）;%求機組數
T=length（p）;%求模擬時間
delp=1;%步長
delh=1;%步長
NP=ceil（max（p）/delp）;???%求電模擬的離散量
NH=ceil（max（h）/delh）;???%求熱模擬的離散量
x1=zeros（1NH）;
x2=zeros（1NH）;
y1=zeros（1NP）;
y2=zeros（1NP）;
W=zeros（3*NP+12*NH+1）;?%第一列和第一行用來存儲熱電負荷坐標
W（NP+1NH+1）=0;
cu=134.5;?%園區與電網聯絡線傳輸功率上限值?各主變容量之和
r=0.0001;%從電網購電也存在一定概率的宕機
xishu=1.7;?%熱電比

for?i=1:NH
????x1（1i）=i*delh;
????x2（1i）=（i-1）*delh;
????x3（1i）=-i*delh;
????x4=flip（x3）;
end
for?i=1:NP
????y1（1i）=i*delp;
????y2（1i）=（i-1）*delp;
????y3（1i）=-i*delp;
????y4=flip（y3）;
????y5=y4-NP*delp;
????y6=[y5y4];
end
W（11:NH）=x4;
W（1:2*NP1）=y6;
W（2*NP+2:3*NP+11）=y1;
W（1NH+2:2*NH+1）=x1;

%形成電熱負荷聯合密度函數
for?i=1:T
????wp=0;
????for?m=1:NH
????????if?h（i）>x2（1m）&&x1（1m）>=h（i）
????????????for?n=1:NP
????????????????if?p（i）>y2（1n）&&y1（1n）>=p（i）?%采取取上限的計數
????????????????????wp=wp+1;
????????????????end
????????????????if?wp==1
????????????????????W（n+1+2*NPm+1+NH）=?W（n+1+2*NPm+1+NH）+wp;
????????????????????wp=0;
????????????????end
????????????end
????????end
????end
end
W（2:end2:end）=W（2:end2:end）/T;

x=-NH:delh:NH;?
y=-2*NP:delp:NP;
?[XY]=meshgrid（xy）;
?for?j1=2:delh:2*NH+1
?????for?j2=2:delp:3*NP+1
?????????WW（j2j1）=W（j2j1）;
?????end
?end
?z=WW;
?mesh（xyz）
?????????
%?對聯合密度函數函數進行修正
E=cell（1jizu+2）;??%用來存放每次修正后的函數
E{11}=W;
sj=0;????????????%實際投入機組數目
for?i=1:jizu????
????W=E{1i};
????if?ceil（prl（i）/delp）==prl（i）/delp
????????k1=ceil（prl（i）/delp）;
????????k2=ceil（hrl（i）/delh）;
????else
????????k1=ceil（prl（i）/delp）-1;
????????k2=ceil（hrl（i）/delh）-1;
????end
????for?t=NH+2:2*NH+1
????????for?j=2:3*NP+1
????????????if?W（1t）????????????????if?ceil（W（1t）/xishu）==W（1t）/xishu
????????????????????qq=ceil（W（1t）/xishu）;
????????????????else
????????????????????qq=ceil（W（1t）/xishu）-1;
????????????????end
????????????????%?????????????????if?W（j1）+qq<=W（2*NP+11）&&2*W（1t）<=W（12*NH+1）
????????????????%?????????????????????W（jt）=gz（i）*W（jt）+（1-gz（i））*W（j+qqt+W（1t））;???%以熱定電，確定出力
????????????????%?????????????????else
????????????????%?????????????????????W（jt）=gz（i）*W（jt）;
????????????????%?????????????????end????
????????????????if?W（jt）~=0;
????????????????W（j-qqNH+1）=?W（j-qqNH+1）+W（jt）;??
????????????????W（jt）=0;
???