What's wrong with my decoupled power flow program

I have write a decoupled power flow program. The code has no problem. However, my results about the reactive power generation are negative. This should be positive. It seems that unreasonable. Can anyone who are experts in power engineering help me,thank you very much?

• m = find(Bus_type ~= 1);
• n = find(Bus_type == 3);
• p = find(Bus_type == 1);
• %Build the B' and B'' matrix
• B1 = inv(imag(Y(m,m)));
• B2 = inv(imag(Y(n,n)));
• %Find the real and reactive power mismatch
• theta_est = zeros(length(ID_bus),1); %Initial Estimate of Angle
• iteration = 200; %Number of iteration
• iter = 1;
• %dP = zeros(length(ID_bus),1);
• %dQ = zeros(length(ID_bus),1);
• %Start iteration
• for iter = 1:iteration
• P = zeros(length(ID_bus),1);
• Q = zeros(length(ID_bus),1);
• for i=1:length(ID_bus)
• for j=1:length(ID_bus)
• P(i) = P(i) + Vhold(i)*Vhold(j)*abs(Y(i,j))*cos(angle(Y(i,j))-theta_est(i)+theta_est(j));
• Q(i) = Q(i) - Vhold(i)*Vhold(j)*abs(Y(i,j))*sin(angle(Y(i,j))-theta_est(i)+theta_est(j));
• end
• end
• dP = P_sch(m) - P(m);
• dQ = Q_sch(n) - Q(n);
• %Find the mismatch Angle
• dtheta = -B1*dP./Vhold(m);
• %Update the angle
• theta_est(m) = theta_est(m) + dtheta;
• %Find the mismatch Voltage
• dV = -B2*dQ./Vhold(n);
• %Update the Voltage
• Vhold(n) = Vhold(n) + dV;
• %Convergence tolerance
• Tolerance = max(max(abs(dP)),max(abs(dQ)));
• end