function Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,Loadi,notLoadi,PGi,PVi,PD0,QD0)
%% 
%QDcos=textread('D:\Project\最小化潮流\最小潮流算例\仙海919PDQDglys.txt');
%QD(QD~=0)=PD(QD~=0)./tan(QDcos);
%QD(QD_NON_ZERO_IND)=QD_NON_ZERO;
%% 
%PD(Loadi)=QD(Loadi)./tan(acos(0.98));
AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum);
Loadi=find(QD0~=0 | PD0~=0);
PDi=find(PD0~=0);
QDi=find(QD0~=0);
notLoadi=setdiff(1:Busnum,Loadi);
PGOnly=setdiff(PGi,PDi);
QGOnly=setdiff(PVi,QDi);
PDOnly=setdiff(PDi,PGi);
QDOnly=setdiff(QDi,PVi);
PGPD=intersect(PDi,PGi);
QGQD=intersect(QDi,PVi);
noPGPDQGQD=setdiff(1:Busnum,union(union(PDi,PGi),union(QDi,PVi)));
YCos=Y.*cos(AngleIJ);
YSin=Y.*sin(AngleIJ);
dP1=PG(ismember(PGi,PGOnly))-diag(Volt(PGOnly))*YCos(PGOnly,:)*Volt(PGOnly);
dQ1=QG(ismember(PVi,QGOnly))-diag(Volt(QGOnly))*YSin(QGOnly,:)*Volt(QGOnly);
dP2=PD(ismember(PDi,PDOnly))-diag(Volt(PDOnly))*YCos(PDOnly,:)*Volt(PDOnly);
dQ2=QD(ismember(QDi,QDOnly))-diag(Volt(QOnly))*YSin(QDOnly,:)*Volt(QDOnly);
dP3=PG(ismember(PGi,PGOnly))-diag(Volt(PGOnly))*YCos(PGOnly,:)*Volt(PGOnly);
dQ=QG-QD-diag(Volt)*Y.*sin(AngleIJ)*Volt;

Mat_H=[dP;dQ;];

end