calu1.m

可以计算电力系统概率潮流

function pui=CalU1(x1, x2, x3, x4, x5, x6, x7, x8, x9)

% calculate Central Moments
afa=zeros(1,9);

afa(1,1)=x1;
afa(1,2)=x2+x1^2;
afa(1,3)=x3+3*x1*x2+x1^3;
afa(1,4)=x4+3*x2^2+4*x1*x3+6*x1^2*x2+x1^4;
afa(1,5)=x5+5*x1*x4+10*x2*x3+10*x1^2*x3+15*x1*x2^2+10*x1^3*x2+x1^5;
afa(1,6)=x6+10*x3^2+15*x2*x4+6*x1*x5+15*x2^3+15*x1^2*x4+60*x1*x2*x3+20*x1^3*x3+45*x1^2*x2^2+15*x1^4*x2+x1^6;
afa(1,7)=x7+7*x1*x6+21*x2*x5+35*x3*x4+21*x1^2*x5+105*x2^2*x3+70*x3^2*x1+105*x1*x2*x4+35*x1^3*x4+105*x2^3*x1+210*x1^2*x2*x3+35*x1^4*x3+105*x1^3*x2^2+21*x1^5*x2+x1^7;
afa(1,8)=x8+35*x4^2+8*x1*x7+28*x2*x6+56*x3*x5+28*x1^2*x6+210*x2^2*x4+280*x3^2*x2+168*x1*x2*x5+280*x1*x3*x4+105*x2^4+56*x1^3*x5+280*x1^2*x3^2+420*x1^2*x2*x4+840*x1*x2^2*x3+70*x1^4*x4+560*x1^3*x2*x3+420*x1^2*x2^3+56*x1^5*x3+210*x1^4*x2^2+28*x1^6*x2+x1^8;
afa(1,9)=x9+9*x1*x8+36*x2*x7+84*x3*x6+126*x4*x5+280*x3^3+36*x1^2*x7+378*x2^2*x5+315*x4^2*x1+252*x1*x2*x6+504*x1*x3*x5+2520*x2*x3*x4+84*x1^3*x6+1260*x2^3*x3+1260*x1^2*x3*x4+756*x1^2*x2*x5+1890*x2^2*x1*x4+2520*x3^2*x1*x2+126*x1^4*x5+945*x2^4*x1+840*x1^3*x3^2+1260*x1^3*x2*x4+3780*x1^2*x2^2*x3+126*x1^5*x4+1260*x1^4*x2*x3+1260*x1^3*x2^3+84*x1^6*x3+378*x1^5*x2^2+36*x1^7*x2+x1^9;

pui=zeros(1,9);
pui(1,1)=0;
pui(1,2)=afa(1,2)-afa(1,1)^2;
pui(1,3)=afa(1,3)-3*afa(1,1)*afa(1,2)+2*afa(1,1)^3;
pui(1,4)=afa(1,4)-4*afa(1,1)*afa(1,3)+6*afa(1,1)^2*afa(1,2)-3*afa(1,1)^4;
pui(1,5)=afa(1,5)-5*afa(1,1)*afa(1,4)+10*afa(1,1)^2*afa(1,3)-10*afa(1,1)^3*afa(1,2)+4*afa(1,1)^5;
pui(1,6)=afa(1,6)-6*afa(1,1)*afa(1,5)+15*afa(1,1)^2*afa(1,4)-20*afa(1,1)^3*afa(1,3)+15*afa(1,1)^4*afa(1,2)-5*afa(1,1)^6;
pui(1,7)=afa(1,7)-7*afa(1,1)*afa(1,6)+21*afa(1,1)^2*afa(1,5)-35*afa(1,1)^3*afa(1,4)+35*afa(1,1)^4*afa(1,3)-21*afa(1,1)^5*afa(1,2)+6*afa(1,1)^7;
pui(1,8)=afa(1,8)-8*afa(1,1)*afa(1,7)+28*afa(1,1)^2*afa(1,6)-56*afa(1,1)^3*afa(1,5)+70*afa(1,1)^4*afa(1,4)-56*afa(1,1)^5*afa(1,3)+28*afa(1,1)^6*afa(1,2)-7*afa(1,1)^8;
pui(1,9)=afa(1,9)-9*afa(1,1)*afa(1,8)+36*afa(1,1)^2*afa(1,7)-84*afa(1,1)^3*afa(1,6)+126*afa(1,1)^4*afa(1,5)-126*afa(1,1)^5*afa(1,4)+84*afa(1,1)^6*afa(1,3)-36*afa(1,1)^7*afa(1,2)+8*afa(1,1)^9;

return