随机生产模拟matlab程序

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;