diff --git a/run.m b/run.m
index fc6b79c..a90264e 100644
--- a/run.m
+++ b/run.m
@@ -127,33 +127,46 @@ fprintf('
 %% 开始进入状态估计
 % clear PD QD PG QG;
 %准备量测量
+sigma=0.03;
 iterPhaseASpotLoadP=phaseASpotLoadP;
 iterPhaseBSpotLoadP=phaseBSpotLoadP;
 iterPhaseCSpotLoadP=phaseCSpotLoadP;
 iterPhaseASpotLoadQ=phaseASpotLoadQ;
 iterPhaseBSpotLoadQ=phaseBSpotLoadQ;
 iterPhaseCSpotLoadQ=phaseCSpotLoadQ;
+%加误差
+mphaseASpotLoadP=phaseASpotLoadP.*(1+normrnd(0,sigma,length(phaseASpotLoadP),1));
+mphaseBSpotLoadP=phaseBSpotLoadP.*(1+normrnd(0,sigma,length(phaseBSpotLoadP),1));
+mphaseCSpotLoadP=phaseCSpotLoadP.*(1+normrnd(0,sigma,length(phaseCSpotLoadP),1));
+mphaseASpotLoadQ=phaseASpotLoadQ.*(1+normrnd(0,sigma,length(phaseASpotLoadQ),1));
+mphaseBSpotLoadQ=phaseBSpotLoadQ.*(1+normrnd(0,sigma,length(phaseBSpotLoadQ),1));
+mphaseCSpotLoadQ=phaseCSpotLoadQ.*(1+normrnd(0,sigma,length(phaseCSpotLoadQ),1));
+%电压量测量
+mVoltpA=VoltpA.*(1+normrnd(0,sigma,length(VoltpA),1));
+mVoltpB=VoltpB.*(1+normrnd(0,sigma,length(VoltpB),1));
+mVoltpC=VoltpC.*(1+normrnd(0,sigma,length(VoltpC),1));
+%转换得到徐电压的量测量
+fV012=Tp2f*conj([mVoltpA';mVoltpB';mVoltpC']);
+mfV0=conj(fV012(1,:)');
+mfV1=conj(fV012(2,:)');
+mfV2=conj(fV012(3,:)');
+% VoltpA=sparse(ones(busNum,1));
+% VoltpB=sparse(ones(busNum,1)).*exp(1j*-120/180*pi);
+% VoltpC=sparse(ones(busNum,1)).*exp(1j*+120/180*pi);
 %全部转换为负荷电流
-VoltpA=sparse(ones(busNum,1));
-VoltpB=sparse(ones(busNum,1)).*exp(1j*-120/180*pi);
-VoltpC=sparse(ones(busNum,1)).*exp(1j*+120/180*pi);
-CurpA=conj((iterPhaseASpotLoadP+1j*iterPhaseASpotLoadQ)./VoltpA);
-CurpB=conj((iterPhaseBSpotLoadP+1j*iterPhaseBSpotLoadQ)./VoltpB);
-CurpC=conj((iterPhaseCSpotLoadP+1j*iterPhaseCSpotLoadQ)./VoltpC);
+mCurpA=conj((mphaseASpotLoadP+1j*mphaseASpotLoadQ)./mVoltpA);
+mCurpB=conj((mphaseBSpotLoadP+1j*mphaseBSpotLoadQ)./mVoltpB);
+mCurpC=conj((mphaseCSpotLoadP+1j*mphaseCSpotLoadQ)./mVoltpC);
 %转换为序电流
-f012=Tp2f*conj([CurpA';CurpB';CurpC']);
+mf012=Tp2f*conj([mCurpA';mCurpB';mCurpC']);
 %把三序电流分离出来
-If0=conj(f012(1,:)');
-If1=conj(f012(2,:)');%负荷电流
-If2=conj(f012(3,:)');
-%试着算一下正序电流
-%     fsY11*V1;
-%形成负荷序电流的测量值
+mIf0=conj(mf012(1,:)');
+mIf1=conj(mf012(2,:)');%负荷电流
+mIf2=conj(mf012(3,:)');
+mIf0=-mIf0;
+mIf1=-mIf1;%mIf1是注入电流,相当于发电机电流
+mIf2=-mIf2;
 %% 设定量测误差
-sigma=0.03;
-mIf0=-If0;
-mIf1=-If1;%mIf1是注入电流,相当于发电机电流
-mIf2=-If2;
 %% 先算正序的
 %平衡节点电流
 fprintf('正序\n');
@@ -161,19 +174,26 @@ 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);
-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) );
+% mfsY11=fsY11;
+% mfsY11(:,Balance)=0;
+% mfsY11(Balance,:)=0;
+% mfsY11=mfsY11+sparse(Balance,Balance,ones(length(Balance),1),busNum,busNum);
+% rV1=inv(mfsY11)*(mIf1)+1;
+% sigmaV1=normrnd(0,sigma,length(Loadi),1);
+% V1measurement=rV1(Loadi).*(1+sigmaV1);
+V1measurement=mfV1(Loadi);
+% wV1r=abs(real( rV1(Loadi).*sigmaV1 ));
+% wV1i=abs(imag( rV1(Loadi).*sigmaV1 ));
+% sigmaI1=normrnd(0,sigma,length(Loadi),1);
+% I1measurement=mIf1(Loadi).*(1+sigmaI1);%测量值是等效发电机电流
+I1measurement=mIf1(Loadi);
+% wI1r=abs( real(mIf1(Loadi).*sigmaI1) );
+% wI1i=abs( imag(mIf1(Loadi).*sigmaI1) );
+%正序测量误差的第二种形式
+wV1r=abs(real(V1measurement)).*sigma;
+wV1i=abs(imag(V1measurement)).*sigma;
+wI1r=abs(real(I1measurement)).*sigma;
+wI1i=abs(imag(I1measurement)).*sigma;
 % [ V1r,V1i,I1r,I1i ]=IPMLoop(V1measurement,wV1r,wV1i,I1measurement,wI1r,wI1i,BalI1r,BalI1i,busNum,Loadi,fsY11,Balance,1 );
 % f=sum(([real(I1measurement);imag(I1measurement)]-[I1r;I1i]).^2)+sum((real(rV1)-V1r).^2)+sum((imag(rV1)-V1i).^2);
 % fprintf('目标值 %f\n',full(f));
@@ -183,19 +203,26 @@ BalI2r=real(-sum(mIf2));
 BalI2i=imag(-sum(mIf2));
 %电压
 %制作量测量
-mfsY22=fsY22;
-mfsY22(:,Balance)=0;
-mfsY22(Balance,:)=0;
-mfsY22=mfsY22+sparse(Balance,Balance,ones(length(Balance),1),busNum,busNum);
-rV2=inv(mfsY22)*(mIf2);
-sigmaV2=normrnd(0,sigma,length(Loadi),1);
-V2measurement=rV2(Loadi).*(1+sigmaV2);
-wV2r=abs(real( rV2(Loadi).*sigmaV2 ));
-wV2i=abs(imag( rV2(Loadi).*sigmaV2 ));
-sigmaI2=normrnd(0,sigma,length(Loadi),1);
-I2measurement=mIf2(Loadi).*(1+sigmaI2);%测量值是等效发电机电流
-wI2r=abs( real(mIf2(Loadi).*sigmaI2) );
-wI2i=abs( imag(mIf2(Loadi).*sigmaI2) );
+% mfsY22=fsY22;
+% mfsY22(:,Balance)=0;
+% mfsY22(Balance,:)=0;
+% mfsY22=mfsY22+sparse(Balance,Balance,ones(length(Balance),1),busNum,busNum);
+% rV2=inv(mfsY22)*(mIf2);
+% sigmaV2=normrnd(0,sigma,length(Loadi),1);
+% V2measurement=rV2(Loadi).*(1+sigmaV2);
+V2measurement=mfV2(Loadi);
+% wV2r=abs(real( rV2(Loadi).*sigmaV2 ));
+% wV2i=abs(imag( rV2(Loadi).*sigmaV2 ));
+% sigmaI2=normrnd(0,sigma,length(Loadi),1);
+% I2measurement=mIf2(Loadi).*(1+sigmaI2);%测量值是等效发电机电流
+I2measurement=mIf2(Loadi);
+% wI2r=abs( real(mIf2(Loadi).*sigmaI2) );
+% wI2i=abs( imag(mIf2(Loadi).*sigmaI2) );
+%负序误差的第二种形式
+wV2r=abs(real(V2measurement)).*sigma;
+wV2i=abs(imag(V2measurement)).*sigma;
+wI2r=abs(real(I2measurement)).*sigma;
+wI2i=abs(imag(I2measurement)).*sigma;
 % [ V2r,V2i,I2r,I2i ]=IPMLoop(V2measurement,wV2r,wV2i,I2measurement,wI2r,wI2i,BalI2r,BalI2i,busNum,Loadi,fsY22,Balance,0 );
 % f=sum(([real(I2measurement);imag(I2measurement)]-[I2r;I2i]).^2)+sum((real(rV2)-V2r).^2)+sum((imag(rV2)-V2i).^2);
 % fprintf('目标值 %f\n',full(f));
@@ -205,19 +232,26 @@ BalI0r=real(-sum(mIf0));
 BalI0i=imag(-sum(mIf0));
 %电压
 %制作量测量
-mfsY00=fsY00;
-mfsY00(:,Balance)=0;
-mfsY00(Balance,:)=0;
-mfsY00=mfsY00+sparse(Balance,Balance,ones(length(Balance),1),busNum,busNum);
-rV0=inv(mfsY00)*(mIf0);
-sigmaV0=normrnd(0,sigma,length(Loadi),1);
-V0measurement=rV0(Loadi).*(1+sigmaV0);
-wV0r=abs(real( rV0(Loadi).*sigmaV0 ));
-wV0i=abs(imag( rV0(Loadi).*sigmaV0 ));
-sigmaI0=normrnd(0,sigma,length(Loadi),1);
-I0measurement=mIf0(Loadi).*(1+sigmaI0);%测量值是等效发电机电流
-wI0r=abs( real(mIf0(Loadi).*sigmaI0) );
-wI0i=abs( imag(mIf0(Loadi).*sigmaI0) );
+% mfsY00=fsY00;
+% mfsY00(:,Balance)=0;
+% mfsY00(Balance,:)=0;
+% mfsY00=mfsY00+sparse(Balance,Balance,ones(length(Balance),1),busNum,busNum);
+% rV0=inv(mfsY00)*(mIf0);
+% sigmaV0=normrnd(0,sigma,length(Loadi),1);
+% V0measurement=rV0(Loadi).*(1+sigmaV0);
+V0measurement=mfV0(Loadi);
+% wV0r=abs(real( rV0(Loadi).*sigmaV0 ));
+% wV0i=abs(imag( rV0(Loadi).*sigmaV0 ));
+% sigmaI0=normrnd(0,sigma,length(Loadi),1);
+% I0measurement=mIf0(Loadi).*(1+sigmaI0);%测量值是等效发电机电流
+I0measurement=mIf0(Loadi);
+% wI0r=abs( real(mIf0(Loadi).*sigmaI0) );
+% wI0i=abs( imag(mIf0(Loadi).*sigmaI0) );
+%零序误差的第二种形式
+wV0r=abs(real(V0measurement)).*sigma;
+wV0i=abs(imag(V0measurement)).*sigma;
+wI0r=abs(real(I0measurement)).*sigma;
+wI0i=abs(imag(I0measurement)).*sigma;
 % matlabpool local 3
 tic
 for II=1:3