蒙特卡洛模拟法计算随机潮流

clear all;
mpc.version = ‘2‘;
%%-----  Power Flow Data  -----%%
%% system MVA base
mpc.baseMVA = 10;   
wind_in=3;
wind=10;
wind_out=25;
p_wind=1;
N=3000;
U= [];
P= [];
Q= [];
U_AVER=[];
P_AVER=[];
Q_AVER=[];
mpopt=mpoption（‘VERBOSE‘0‘OUT_ALL‘0）;
for i = 1:1:N
%% bus data
%	bbus_i	type	Pd			Qd		Gs	Bs	 area	Vm	 Va	 baseKV	zone  Vmax	Vmin
V37 = wblrnd（10.8182.842）;
if（V37 > wind_out || V37 < wind_in）
    P37=0;
elseif（wind <=  V37 && V37 <= wind_out） 
    P37=p_wind;
else
    P37=（V37-wind_in）/（wind-wind_in）*p_wind;
end
V38 = wblrnd（10.8182.842）;
if（V38 > wind_out || V38 < wind_in）
    P38=0;
elseif（wind <=  V38 && V38 <= wind_out） 
    P38=p_wind;
else
    P38=（V38-wind_in）/（wind-wind_in）*p_wind;
end
V39 = wblrnd（10.8182.842）;
if（V39 > wind_out || V39 < wind_in）
    P39=0;
elseif（wind <=  V39 && V39 <= wind_out） 
    P39=p_wind;
else
    P39=（V39-wind_in）/（wind-wind_in）*p_wind;
end
P2  = normrnd（0.100 0.04）;
P3  = normrnd（0.090 0.04）;
P4  = normrnd（0.120 0.04）;
P5  = normrnd（0.060 0.04）;
P6  = normrnd（0.060 0.04）;
P7  = normrnd（0.200 0.04）;
P8  = normrnd（0.200 0.04）;
P9  = normrnd（0.060 0.04）;
P10 = normrnd（0.060 0.04）;
P11 = normrnd（0.045 0.04）;
P12 = normrnd（0.060 0.04）;
P13 = normrnd（0.060 0.04）;
P14 = normrnd（0.120 0.04）;
P15 = normrnd（0.060 0.04）;
P16 = normrnd（0.060 0.04）;
P17 = normrnd（0.060 0.04）;
P18 = normrnd（0.090 0.04）;
P19 = normrnd（0.090 0.04）;
P20 = normrnd（0.090 0.04）;
P21 = normrnd（0.090 0.04）;
P22 = normrnd（0.090 0.04）;
P23 = normrnd（0.090 0.04）;
P24 = normrnd（0.420 0.04）;
P25 = normrnd（0.420 0.04）;
P26 = normrnd（0.060 0.04）;
P27 = normrnd（0.060 0.04）;
P28 = normrnd（0.060 0.04）;
P29 = normrnd（0.120 0.04）;
P30 = normrnd（0.200 0.04）;
P31 = normrnd（0.150 0.04）;
P32 = normrnd（0.210 0.04）;
P33 = normrnd（0.060 0.04）;

 Q2   = normrnd（0.060 0.01）;
 Q3   = normrnd（0.040 0.01）;
 Q4   = normrnd（0.080 0.01）;
 Q5   = normrnd（0.030 0.01）;
 Q6   = normrnd（0.020 0.01）;
 Q7   = normrnd（0.100 0.01）;
 Q8   = normrnd（0.100 0.01）;
 Q9   = normrnd（0.020 0.01）;
 Q10  = normrnd（0.020 0.01）;
 Q11  = normrnd（0.030 0.01）;
 Q12  = normrnd（0.035 0.01）;
 Q13  = normrnd（0.035 0.01）;
 Q14  = normrnd（0.080 0.01）;
 Q15  = normrnd（0.010 0.01）;
 Q16  = normrnd（0.020 0.01）;
 Q17  = normrnd（0.020 0.01）;
 Q18  = normrnd（0.040 0.01）;
 Q19  = normrnd（0.040 0.01）;
 Q20  = normrnd（0.040 0.01）;
 Q21  = normrnd（0.040 0.01）;
 Q22  = normrnd（0.040 0.01）;
 Q23  = normrnd（0.050 0.01）;
 Q24  = normrnd（0.200 0.01）;
 Q25  = normrnd（0.200 0.01）;
 Q26  = normrnd（0.025 0.01）;
 Q27  = normrnd（0.025 0.01）;
 Q28  = normrnd（0.020 0.01）;
 Q29  = normrnd（0.070 0.01）;
 Q30  = normrnd（0.600 0.01）;
 Q31  = normrnd（0.070 0.01）;
 Q32  = normrnd（0.100 0.01）;
 Q33  = normrnd（0.040 0.01）;

mpc.bus = [
 	1	     3      0	  	    0	     0  0     1     1.0  0   12.66   1     1.1      0.9;
	2        1      P2  	    Q2     0