function [jaco]=Jacobi(Balance,busNum,QGi,Volt,YdotSinVolt,YdotCosVolt,diag_Volt_YdotSin,diag_Volt_YdotCos)
diag_YdotSinVolt_=diag(YdotSinVolt);
diag_YdotCosVolt_=diag(YdotCosVolt);
dPdTyta=diag_Volt_YdotSin*diag(Volt)-diag_YdotSinVolt_*diag(Volt); % 简化第三次
dQdTyta=-diag_Volt_YdotCos*diag(Volt)+diag_YdotCosVolt_*diag(Volt);%dQ/dThyta
dPdV=diag_YdotCosVolt_+diag_Volt_YdotCos;%dP/dV
dQdV=diag_YdotSinVolt_+diag_Volt_YdotSin;%dQ/dV
%% 平衡节点相角不变
dPdTyta(Balance,:)=0;
dPdTyta(:,Balance)=0;
dPdTyta=dPdTyta+sparse(Balance,Balance,ones(length(Balance),1),busNum,busNum);
dQdTyta(:,Balance)=0;
dPdV(Balance,:)=0;
%% PV节点电压不变
dQdV(QGi,:)=0;
dQdV(:,QGi)=0;
dQdV=dQdV+sparse(QGi,QGi,ones(length(QGi),1),busNum,busNum);
dQdTyta(QGi,:)=0;
dPdV(:,QGi)=0;
% jaco=[    
%     dPdV,dPdTyta;
%     dQdV,dQdTyta;
% ];
jaco=cat(1,cat(2,dPdV,dPdTyta),cat(2,dQdV,dQdTyta));
end