diff --git a/FormG.m b/FormG.m
index 7dd6163..045ef32 100644
--- a/FormG.m
+++ b/FormG.m
@@ -1,8 +1,8 @@
-function Mat_G=FormG(Volt,PD,QD,Loadi,bigM,mVolt,rVolt,sigma,Vbi,PDbi,QDbi)
+function Mat_G=FormG(Volt,PD,QD,Loadi,bigM,mVolt,rVolt,sigma,Vbi,PDbi,QDbi,mPD)
 
 Mat_G=[
-    sparse(PD(Loadi));
-    sparse(PD(Loadi))-0.001;
+    sparse(PD(Loadi))-mPD(Loadi)-bigM*PDbi-0.1;
+    sparse(PD(Loadi))-mPD(Loadi)+bigM*PDbi+0.1;
     sparse(QD(Loadi));
     sparse(QD(Loadi))-0.001;
     Volt'-mVolt'-bigM*Vbi-0.03;
diff --git a/FormLw.m b/FormLw.m
index c91cad8..5f36d4f 100644
--- a/FormLw.m
+++ b/FormLw.m
@@ -24,7 +24,7 @@ QDU(indQD(3:12:end))=1.55*realQD(indQD(3:12:end));
 QDU(indQD(9:12:end))=1.05*realQD(indQD(9:12:end));
 % PF=0.85;
 % QDU=1.0*PD(Loadi).*sqrt(1 -PF.^2)./PF;
-t1=([PDU',PDU',QDU',QDU',0*VoltU,100*ones(1,length(VoltU)),1*ones(1,Busnum+length(Loadi)*2)])';
+t1=([0*PDU',100*ones(1,length(PDU)),QDU',QDU',0*VoltU,100*ones(1,length(VoltU)),1*ones(1,Busnum+length(Loadi)*2)])';
 t2=Mat_G+Init_U'-t1;
 Lw=t2;
 
diff --git a/FormLz.m b/FormLz.m
index a87dcde..d664aef 100644
--- a/FormLz.m
+++ b/FormLz.m
@@ -21,7 +21,7 @@ indQD=find(realQD>0);
 QDL(indQD(3:12:end))=0.95*realQD(indQD(3:12:end));
 QDL(indQD(9:12:end))=0.95*realQD(indQD(9:12:end));
 % QDL=0*PD(Loadi).*sqrt((1-PF.^2))./PF;
-t1=([PDL',PDL',QDL',QDL',-100*ones(1,length(VoltL)),0*VoltL,0*ones(1,Busnum+length(Loadi)*2)])';
+t1=([-100*ones(1,length(PDL)),0*PDL',QDL',QDL',-100*ones(1,length(VoltL)),0*VoltL,0*ones(1,Busnum+length(Loadi)*2)])';
 t2=Mat_G-Init_L'-t1;
 Lz=t2;
 
diff --git a/OPF.m b/OPF.m
index 8ee5c33..8d99300 100644
--- a/OPF.m
+++ b/OPF.m
@@ -52,6 +52,8 @@ Precision=Precision/1;
 sigma=0.01;
 PD0(Loadi)=PD0(Loadi).*(1+normrnd(0,sigma,length(Loadi),1));
 QD0(Loadi)=QD0(Loadi).*(1+normrnd(0,sigma,length(Loadi),1));
+mPD=PD0;
+mQD=QD0;
 mVolt=rVolt.*(1+normrnd(0,sigma,length(rVolt),1))';
 bigM=10;
 Vbi=sparse(ones(Busnum,1));
@@ -83,7 +85,7 @@ while(abs(Gap)>Precision)
     %% Ðγɷ½³Ì¾ØÕó
     Luu=Init_U'.*Init_W'+Init_u*ones(RestraintCount,1);
     Lul=Init_L'.*Init_Z'-Init_u*ones(RestraintCount,1);
-    Mat_G=FormG(Volt,PD,QD,Loadi,bigM,mVolt,rVolt,sigma,Vbi,PDbi,QDbi);
+    Mat_G=FormG(Volt,PD,QD,Loadi,bigM,mVolt,rVolt,sigma,Vbi,PDbi,QDbi,mPD);
     Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,Loadi);
     Ly=Mat_H;
     Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD,RealPD,RealQD,Loadi,KK,PF);
diff --git a/func_deltG.m b/func_deltG.m
index fd8dac2..b3e5ac2 100644
--- a/func_deltG.m
+++ b/func_deltG.m
@@ -44,8 +44,8 @@ dg7_dvbi=sparse(eye(Busnum,Busnum));
 dgPDbi_dvbi=sparse(Busnum,sizeLoadi);
 dgQDbi_dvbi=sparse(Busnum,sizeLoadi);
 %%
-dg3_dPDbi=sparse(-zeros(sizeLoadi));
-dg32_dPDbi=sparse(zeros(sizeLoadi));
+dg3_dPDbi=sparse(-eye(sizeLoadi));
+dg32_dPDbi=sparse(eye(sizeLoadi));
 dg4_dPDbi=sparse(length(Loadi),length(Loadi));
 dg42_dPDbi=dg4_dPDbi;
 dg5_dPDbi=sparse(length(Loadi),Busnum);