diff --git a/FormLw.m b/FormLw.m index 42081c5..ab128f2 100644 --- a/FormLw.m +++ b/FormLw.m @@ -1,5 +1,5 @@ function Lw=FormLw(Loadi,Mat_G,Init_U) -upper=0.01*sparse(ones(length(Loadi)*1,1)); +upper=0.03*sparse(ones(length(Loadi)*1,1)); Lw=Mat_G+Init_U-upper; end \ No newline at end of file diff --git a/IPMLoop.m b/IPMLoop.m index 9332c8d..9ca35ee 100644 --- a/IPMLoop.m +++ b/IPMLoop.m @@ -1,4 +1,4 @@ -function [ V1r,V1i,I1r,I1i ] = IPMLoop(Vmeasurement,Imeasurement,BalI1r,BalI1i,busNum,Loadi,fsY1,Balance,Vref ) +function [ V1r,V1i,I1r,I1i ] = IPMLoop(V1measurement,wV1r,wV1i,I1measurement,wI1r,wI1i,BalI1r,BalI1i,busNum,Loadi,fsY1,Balance,Vref ) %把每个序的循环写在这个函数中。其实也就是内点法循环。 V1r=1*ones(busNum,1); V1i=0*ones(busNum,1); @@ -33,7 +33,7 @@ while(abs(Gap)>0.000001) L_1Z=diag(Init_Z./Init_L); U_1W=diag(Init_W./Init_U); %% 形成海森阵 - deltdeltF=func_deltdeltF(busNum,fsY1,Loadi); + deltdeltF=func_deltdeltF(busNum,fsY1,Loadi,wV1r,wV1i,wI1r,wI1i); % deltdeltF=0; %% 形成ddHy % ddh=func_ddh(busNum,Loadi,Init_Z,Init_W); @@ -41,7 +41,7 @@ while(abs(Gap)>0.000001) %% 开始构建ddg ddg=func_ddg(busNum,Loadi,Init_Z,Init_W); %% 开始构建deltF - deltF=func_deltF(Vmeasurement,Imeasurement,busNum,fsY1,Loadi,V1r,V1i,I1r,I1i); + deltF=func_deltF(V1measurement,wV1r,wV1i,I1measurement,wI1r,wI1i,busNum,fsY1,Loadi,V1r,V1i,I1r,I1i); % deltF=0; %% Luu=Init_U.*Init_W+Init_u*ones(RestraintCount,1); @@ -62,7 +62,7 @@ while(abs(Gap)>0.000001) fprintf('Gap %f\n',full(Gap)); KK=KK+1; end -f=sum(([real(Imeasurement);imag(Imeasurement)]-[-I1r;-I1i]).^2)+sum((real(Vmeasurement)-V1r(Loadi)).^2)+sum((imag(Vmeasurement)-V1i(Loadi)).^2); +f=sum(([real(I1measurement);imag(I1measurement)]-[-I1r;-I1i]).^2)+sum((real(V1measurement)-V1r(Loadi)).^2)+sum((imag(V1measurement)-V1i(Loadi)).^2); end diff --git a/func_deltF.m b/func_deltF.m index 1cde7a6..b5ca5ad 100644 --- a/func_deltF.m +++ b/func_deltF.m @@ -1,4 +1,4 @@ -function deltF=func_deltF(Vmeasurement,Imeasurement,busNum,fsY11,Loadi,V1r,V1i,I1r,I1i) +function deltF=func_deltF(V1measurement,wV1r,wV1i,I1measurement,wI1r,wI1i,busNum,fsY11,Loadi,V1r,V1i,I1r,I1i) %t1=PG(setdiff(PVi,Balance)); % t2=Volt'*Volt; % t3=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat); @@ -9,9 +9,11 @@ function deltF=func_deltF(Vmeasurement,Imeasurement,busNum,fsY11,Loadi,V1r,V1i,I %% deltF=[ %zeros(busNum*2,1); - sparse(Loadi,1,-2*(real(Vmeasurement)-V1r(Loadi)),busNum,1); - sparse(Loadi,1,-2*(imag(Vmeasurement)-V1i(Loadi)),busNum,1); - -2*( [real(Imeasurement);imag(Imeasurement)]-[I1r;I1i]); + sparse(Loadi,1,-2*(real(V1measurement)-V1r(Loadi))./wV1r./wV1r,busNum,1); + sparse(Loadi,1,-2*(imag(V1measurement)-V1i(Loadi))./wV1i./wV1i,busNum,1); + %-2*( [real(Imeasurement);imag(Imeasurement)]-[I1r;I1i]); + sparse(1:length(Loadi),1,-2*(real(I1measurement)-V1r(Loadi))./wI1r./wI1r,length(Loadi),1); + sparse(1:length(Loadi),1,-2*(imag(I1measurement)-V1i(Loadi))./wI1i./wI1i,length(Loadi),1); ]; end \ No newline at end of file diff --git a/func_deltdeltF.m b/func_deltdeltF.m index e687865..b034cf9 100644 --- a/func_deltdeltF.m +++ b/func_deltdeltF.m @@ -1,10 +1,10 @@ -function deltdeltF=func_deltdeltF(busNum,fsY11,Loadi) +function deltdeltF=func_deltdeltF(busNum,fsY11,Loadi,wV1r,wV1i,wI1r,wI1i) deltdeltF=[ %zeros(busNum*2,busNum*2+length(Loadi)*2); - sparse(Loadi,Loadi,2,busNum,busNum*2+length(Loadi)*2); - sparse(Loadi,busNum+Loadi,2,busNum,busNum*2+length(Loadi)*2); - zeros(length(Loadi)*2,busNum*2),2*eye(length(Loadi)*2); + sparse(Loadi,Loadi,2./[wV1r.*wV1r],busNum,busNum*2+length(Loadi)*2); + sparse(Loadi,busNum+Loadi,2./[wV1i.*wV1i],busNum,busNum*2+length(Loadi)*2); + zeros(length(Loadi)*2,busNum*2),2*eye(length(Loadi)*2)*diag(1./[wI1r.^2;wI1i.^2]); %sparse(Loadi,busNum*2+Loadi,2,busNum,busNum*2+length(Loadi)*2); ]; end \ No newline at end of file diff --git a/run.m b/run.m index cc9cf4e..2262c58 100644 --- a/run.m +++ b/run.m @@ -221,6 +221,8 @@ If2=conj(f012(3,:)'); %试着算一下正序电流 % fsY11*V1; %形成负荷序电流的测量值 +%% 设定量测误差 +sigma=0.03; mIf0=-If0; mIf1=-If1; % mIf1(3)=-mIf1(2); @@ -236,31 +238,40 @@ fprintf(' BalI1r=real(-sum(mIf1)); BalI1i=imag(-sum(mIf1)); %电压 +%制作量测量 mfsY11=fsY11; mfsY11(:,Balance)=0; mfsY11(Balance,:)=0; mfsY11=mfsY11+sparse(Balance,Balance,ones(length(Balance),1),busNum,busNum); -mV1=inv(mfsY11)*(mIf1)+1; -Vmeasurement=mV1(Loadi); -Imeasurement=-mIf1(Loadi); -[ V1r,V1i,I1r,I1i ]=IPMLoop(Vmeasurement,Imeasurement,BalI1r,BalI1i,busNum,Loadi,fsY11,Balance,1 ); -f=sum(([real(Imeasurement);imag(Imeasurement)]-[-I1r;-I1i]).^2)+sum((real(Vmeasurement)-V1r(Loadi)).^2)+sum((imag(Vmeasurement)-V1i(Loadi)).^2); +rV1=inv(mfsY11)*(mIf1)+1; +sigmaV1=normrnd(0,sigma,length(Loadi),1); +V1measurement=rV1(Loadi).*(1+sigmaV1); +wV1r=abs(real( rV1(Loadi).*sigmaV1 )); +wV1i=abs(imag( rV1(Loadi).*sigmaV1 )); +sigmaI1=normrnd(0,sigma,length(Loadi),1); +I1measurement=-mIf1(Loadi).*(1+sigmaI1); +wI1r=abs( real(-mIf1(Loadi).*sigmaI1) ); +wI1i=abs( imag(-mIf1(Loadi).*sigmaI1) ); + +[ V1r,V1i,I1r,I1i ]=IPMLoop(V1measurement,wV1r,wV1i,I1measurement,wI1r,wI1i,BalI1r,BalI1i,busNum,Loadi,fsY11,Balance,1 ); +% f=sum(([real(Imeasurement);imag(Imeasurement)]-[-I1r;-I1i]).^2)+sum((real(Vmeasurement)-V1r(Loadi)).^2)+sum((imag(Vmeasurement)-V1i(Loadi)).^2); +f=sum(([real(-mIf1(Loadi));imag(-mIf1(Loadi))]-[-I1r;-I1i]).^2)+sum((real(rV1(Loadi))-V1r(Loadi)).^2)+sum((imag(rV1(Loadi))-V1i(Loadi)).^2); fprintf('目标值 %f\n',full(f)); %% 算负序的 fprintf('负序\n'); BalI2r=real(-sum(mIf2)); BalI2i=imag(-sum(mIf2)); -measurement=-mIf2(Loadi); -[ V2r,V2i,I2r,I2i ]=IPMLoop(measurement,BalI2r,BalI2i,busNum,Loadi,fsY22,Balance,0 ); -f=sum([real(measurement);imag(measurement)]-[-I2r;-I2i]); +Imeasurement=-mIf2(Loadi); +[ V2r,V2i,I2r,I2i ]=IPMLoop(Vmeasurement,Imeasurement,BalI1r,BalI1i,busNum,Loadi,fsY1,Balance,0 ); +f=sum(([real(Imeasurement);imag(Imeasurement)]-[-I1r;-I1i]).^2)+sum((real(Vmeasurement)-V2r(Loadi)).^2)+sum((imag(Vmeasurement)-V2i(Loadi)).^2); fprintf('目标值 %f\n',full(f)); %% 算零序 fprintf('零序\n'); BalI0r=real(-sum(mIf0)); BalI0i=imag(-sum(mIf0)); -measurement=-mIf0(Loadi); -[ V0r,V0i,I0r,I0i ]=IPMLoop(measurement,BalI0r,BalI0i,busNum,Loadi,fsY00,Balance,0 ); -f=sum([real(measurement);imag(measurement)]-[-I0r;-I0i]); +Imeasurement=-mIf0(Loadi); +[ V0r,V0i,I0r,I0i ]=IPMLoop(Vmeasurement,Imeasurement,BalI1r,BalI1i,busNum,Loadi,fsY1,Balance,0 ); +f=sum(([real(Imeasurement);imag(Imeasurement)]-[-I1r;-I1i]).^2)+sum((real(Vmeasurement)-V0r(Loadi)).^2)+sum((imag(Vmeasurement)-V0i(Loadi)).^2); fprintf('目标值 %f\n',full(f)); %状态量 % SEVoltpA=sparse(ones(busNum,1));