diff --git a/MaxDeviation.m b/MaxDeviation.m
index ad6f95a..2304770 100644
--- a/MaxDeviation.m
+++ b/MaxDeviation.m
@@ -1,8 +1,9 @@
-function [output_arg]=MaxDeviation(rVolt,SEVolt,rVAngel,SEVAngel)
-t1=[rVolt;rVAngel];
-t2=[SEVolt;SEVAngel];
+function [output_arg]=MaxDeviation(rVolt,SEVolt,rVAngel,SEVAngel,rPD,rQD,PD,QD)
+t1=[rVolt;rVAngel;rPD;rQD];
+t2=[SEVolt;SEVAngel;PD;QD];
 t3=(t1(t1~=0)-double(t2(t1~=0)))./t1(t1~=0);
 t4=abs(t3);
+
 output_arg=max(t4);
 end
 
diff --git a/Run_YALMIP.m b/Run_YALMIP.m
index 3533c5c..7519747 100644
--- a/Run_YALMIP.m
+++ b/Run_YALMIP.m
@@ -36,6 +36,8 @@ xVolt=Volt;
 xUAngel=UAngel;
 % VMatrix=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum);
 % dP=PG-PD-diag(xVolt)*(Y.*cos(VMatrix))*xVolt';
+rPD=PD;
+rQD=QD;
 [Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wPG,wQG,wPD,wQD,wVolt,PD,PD0,QD,randPDind,Loadi,notLoadi]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0,QD0,QD,PD);
 %% 定义变量
 BalVolt=Volt(Balance);
@@ -91,4 +93,6 @@ rVAngel=xUAngel';
 SEVolt=x(length(PDi)+length(QDi)+1:length(Volt)+length(PDi)+length(QDi));
 SEVAngel=x(length(Volt)+length(PDi)+length(QDi)+1:end);
 fprintf('最大偏差\n')
-MaxDeviation(rVolt,SEVolt,rVAngel,SEVAngel)
+PD=x(1:length(PDi));
+QD=x(length(PDi)+1:length(PDi)+length(QDi));
+MaxDeviation(rVolt,SEVolt,rVAngel,SEVAngel,rPD(PDi),rQD(QDi),PD,QD)
diff --git a/StatDeviation.m b/StatDeviation.m
new file mode 100644
index 0000000..5e5c442
--- /dev/null
+++ b/StatDeviation.m
@@ -0,0 +1,20 @@
+function [ output_args ] = StatDeviation( this,PG0,QG0,PD0,QD0 )%统计误差
+%STATDEVIATION Summary of this function goes here
+%   Detailed explanation goes here
+
+PD0Array=repmat(PD0,1,this.sampleNum);
+QD0Array=repmat(QD0,1,this.sampleNum);
+PDDev=(this.PDArray-PD0Array)/0.05;
+QDDev=(this.QDArray-QD0Array)/0.05;
+% PG0Array=repmat(PG0,this.sampleNum,1);
+% QG0Array=repmat(QG0,this.sampleNum,1);
+% PGDev=(PG0Array-this.PGArray)/0.01;
+% QGDev=(QG0Array-this.QGArray)/0.01;
+wholeMat=[PDDev;QDDev;];
+t1=wholeMat.^2;
+t2=sum(t1,1);
+t3=t2/size(t1,1);
+t4=t3.^.5;
+output_args=sum(t4)/length(t4);
+end
+