function [line_flow, bus]=loadlflow(line, bus, tolerance, max_iter, min_condition)
% [line_flow, bus_out]=loadlflow(line, bus, tolerance, max_iter, min_conditiion)
% Solve a load flow for a given electrical system (buses and lines)
%
% Usage:
%
%  [line_flow, bus_out]=loadflow(line, bus)
%  [line_flow, bus_out]=loadflow(line, bus, tolerance)
%  [line_flow, bus_out]=loadflow(line, bus, tolerance, max_iter)
%  [line_flow, bus_out]=loadflow(line, bus, tolerance, max_iter, min_condition)
%
% If un-specified defaults are:
%    tolerance = 1e-5
%    max_iter = 5
%    min_condition = 1e-5
%
% Inputs:
%   line / Line Matrix:
%      Column  Name          Meaning
%      1       from          ID of bus connected to one end
%      2       to            ID of bus on the other end
%      3       resistance    Resistance of the line in per unit (R)
%      4       reactance     reactance of the line in per unit (X)
%      5       susceptance   Susceptance of the line in per unit (line charging) 
%                            Models the overall ground capacitative effects of
%                            a medium length line (50-150 miles). (B)
%                            Modeled as a pi-connected line See [1]
%      6       tap           Transformer tap ratio (Step up and down ratio)
%      7       tapAngle      Transformer phase shift in degrees
%      8       Pmax          Max Real Power (Magnitude)
%      9       Smax          Max Apparent Power (Magnitude)
%      10      factsSkip     Skip line in FACTS Device placement 
%      11      outageSkip    Skip line outages (never outage)
%  
%   bus / Bus Matrix:
%      Column  Name          Meaning
%      1       number        ID of the bus (Sequential, with no Gaps 1:1:N)
%      2       voltage       Voltage on the bus (initial)
%      3       theta         Voltage Angle in degrees (initial)
%      4       Pgen          P generated on the bus
%      5       Qgen          Q generated on the bus
%      6       Pload         P load on the bus
%      7       Qload         Q load on the bus
%      8       G             Shunt Conductance ("resistance" to ground)
%      9       B             Shunt Susceptance ("reactance" to ground)
%      10      type          Type of bus Matlab: 1=Swing, 2=Generator, 3=Normal
%                            Swing: theta and V are known (The reference bus)
%                                   (P & Q are unknown)
%                            Generator: V and P are known
%                                       (Q & theta are unknown)
%                            Normal/Load: P & Q are known
%                                         (V & theta are unknown)
%
%   tolerance:     Convergence tolerance for power mis-matches
%
%   max_iter:      Maximum number of iterations before failure
%
%   min_condition: Minimum acceptable reciprical-condition number estimate on Jacobian
%
% Return Values:
%   line_flow / Line Flows:
%      Column  Name            Meaning
%      1       from            ID of bus connected to one end
%      2       to              ID of bus on the other end
%      3       P_from2to       Real Power going from "from bus" to "to bus"
%      4       P_to2from       Real Power going from "to bus" to "from bus"
%      5       Q_from2to       Reactive Power going from "from bus" to "to bus"
%      6       Q_to2from       Reactive Power going from "to bus" to "from bus"
%      7       S_mag           Maximum apparent power magnitued
%                              (max(|S(from->to)|, |S(to->From)|))
%  Error Conditions:
%      Message ID                Error Message            
%      'loadflow:ill_condition'  'Ill conditioned Jacobian'
%      'loadflow:max_iterations' 'Max Iterations Exceeded'
%
%    Errors can be caught with "try/catch" blocks
%    Ex: 
%     try
%        [line_flow, bus_out]=loadflow(line,bus);
%        ... % Analyze / record results
%     catch
%       err=lasterror;
%       if err.identifier == 'loadflow:ill_condition'
%         disp('Ill Conditioned Jacobian. Disregarding results.');
%       end
%     end
%
%  Bus Model:
%        
%           -------------------------   <- The Bus (Voltage, Theta, P, Q)
%            |   |   ^     |
%          B - G \   |     |
%            -   /   |     |
%            |   \   |     |
%            |   /   |     |
%            |   |   |     V
%            Gnd Gnd Pgen  Pload
%                    Qgen  Qload
%    Each bus has a Voltage, Phase Angle (theta), Real Power, and Reactive power
%    as well as a bus type. There are 3 types of buses:
%      Type 1 - The "swing bus"
%         The swing bus is assumed to be the "reference" for the phase angle.
%         (I.e. it's theta is 0 by definition). This is always a form of 
%         Generator bus, and usually chosen as the generator with the largest
%         output in the system. Since it's a generator, it's Voltage is also 
%         assumed to be known. Available real and reactive power are unknown.
%         Type of bus Matlab: 1=Swing, 2=Generator, 3=Normal
%         (theta and V are known / fixed, P & Q are unknown / solved for.
%          Initial P and Q are "guesses" of correct values)
%
%      Type 2 - The "PV buses" / Generator buses
%         The PV buses represent all the remaining generators in the system
%         Since they are generating and contributing power, the voltage and 
%         available power are assumed to be known and the phase angle and available
%         reactive power are unknown (V and P are known / fixed, Q & theta are
%         unknown / solved for. Initial P and Q are "guesses" of correct values)
%
%      Type 3 - The "PQ buses" / Load buses
%         The load buses represent the normal consumers of power (the loads).
%         At the loads the real and reactive power needs are known, but the 
%         bus voltage and phase angle are unknown (P & Q are known / fixed, 
%         V & theta are unknown / solved for. Initial V and theta are "guesses"
%         of correct values). At these buses, Pload is greater than Pgen
%
%      B - The "susceptance" to ground. (Reciprical of reactance)
%          This represents the losses due to the interaction of the ground
%          and the AC signal on the bus.
%
%      G - The "conductance" to ground. (Reciprical of resistance)
%          This represents the losses due to AC signal to ground.
%
%      Pgen / Qgen - These represent the real and reactive power being 
%          injected into the bus.
%
%      Pload / Qload - These represent the real and reactive power being
%          consumed at the bus.
%
%  Line Model:
%              
%   tap @ tapAngle  
%       |                                                   |
%       --->   3---------/\/\/\/\-^^^^^^^-------------------|
%   |          3    |        R       X                   |  |
%   |----------3    -                                    -  |
%   |          3    - B/2                            B/2 -  |
%   |          3    |                                    |  To Bus
%   |          |    |                                    |
%   From Bus   Gnd  Gnd                                  Gnd
%
%      R - Resistance of the line in per unit
%
%      X - Reactance of the line in per unit
%
%      B - The susceptance (1/capacitance) of the line to ground (shunt)
%          This is a pi-modeled line where the total capacitive effect 
%          to ground is modeled by a capacitance of half the total 
%          effective capacitance is placed on each side of the line.
%          This is an accurate for a medium length lines (50-150 miles).
%
%      tap - Multiplier to step up (or down) voltage for a transformer
%
%      tapAngle - The phase shift imposed by a phase shifting transformer.
%
% References: 
%   1. Power System Analysis: Short-Circuit Load Flow and Harmonics
%      J.C. Das. Dekker. 2002.
%   2. Computational Methods for Electric Power Systems
%      Mariesa Crow. CRC Press. 2003
%   3. Computer Analysis of Power System
%      J. Arrillaga and C.P. Arnold. Wiley. 1994
% Copyright 2005 William M. Siever All Rights Reserved

  % I. Initialize constants: 
    % A. Minimum allowed Jacobian reciporical condition number estimate
    if nargin<5
      min_condition=1e-5; 
    end
    % B. Maximum allowed iterations
    if nargin<4
      max_iter=5;
    end
    % C. Tolerance of P and Q mismatch
    if nargin<3
      tolerance=1e-5;
    end

   % II. Do some simple data conversions
     % A. Convert the bus angle to radians
     bus(:,3)=bus(:,3)*pi/180;
     % B. Convert the line angle to radians
     line(:,7)=line(:,7)*pi/180;
     % C. If transformer tap is 0, make it a 1
     line(:,6)=line(:,6)+(line(:,6)==0);
	      	   	  
  % III. Setup for Newton-Raphson Iteration:
    % A. Get the admittance Matrix
    [y_mag, y_ang]=ymat(line,bus);
    % B. Find the buses that will have a deriviative with respect to 
    %    phase angle (Everything but the reference buses)
    angleMap=find(bus(:,10)~=1);
	      	   	  
    % C. Find the buses that will have a deriviative with respect to 
    %    voltage (the load buses)    
    vMap=find(bus(:,10)==3);

    % D. Find the buses with generators
    genMap=find(bus(:,10)~=3);
    
    % E. Find the swing buses
    swingMap=find(bus(:,10)==1);
	      	   	  
  % IV. Iteration
    for i=1:max_iter
      % A. Compute the Jacobian:
      [jac, deltaP, deltaQ]=jmat(bus, y_mag, y_ang);
      % B. Find complete P mismatch:
      %    Pmismatch = P(generation) - P(load) - deltaP
      Pmismatch=bus(:,4)-bus(:,6)-deltaP;
      % C. Find complete Q mismatch:
      %    Qmismatch = Q(generation) - Q(load) - deltaQ
      Qmismatch=bus(:,5)-bus(:,7)-deltaQ;

      % Update the Qgen on all generator buses
      bus(genMap,5)=bus(genMap,5) + -Qmismatch(genMap);
      % Update the Pgen on the swing bus
      bus(swingMap,4)=bus(swingMap,4) + -Pmismatch(swingMap);

      % D. Determine if current tolerance is acceptable
      if all(abs(Pmismatch) < tolerance) & ...
         all(abs(Qmismatch) < tolerance)
	      	   	  
        % A. Get From Voltages and Angles
        Vf = bus(line(:,1),2);
        Vf_ang = bus(line(:,1),3);
        % B. Get To Voltages and Angles
        Vt = bus(line(:,2),2);
        Vt_ang = bus(line(:,2),3);
        % C. Compute the From Transformer phase shift voltage multiplier
        Vf_tps = Vf./line(:,6);
         
        % D. Get each line's Admittance Magnitude and Angle
        y_line_mag=(1./sqrt(line(:,3).^2+line(:,4).^2));
        y_line_ang=-atan2(line(:,4),line(:,3));
        % E. Compute the Sin and Cos terms of the Admittance Angles
        cos_y_line_ang = cos(-y_line_ang);
        sin_y_line_ang = sin(-y_line_ang);
         
        % F. Compute the from and to angles
        Vf_delta_ang = Vf_ang - line(:,7) - Vt_ang - y_line_ang;
        Vt_delta_ang = Vt_ang + line(:,7) - Vf_ang - y_line_ang;
	      	   	  
        % G. Compute the real power flow
        Pflow_f2t=Vf_tps.*y_line_mag.* ...
                  (Vf_tps.*cos_y_line_ang - Vt.*cos(Vf_delta_ang));
        Pflow_t2f=Vt.*y_line_mag.* ...
                  (Vt.*cos_y_line_ang - Vf_tps.*cos(Vt_delta_ang));
	      	   	  
        % H. Compute the reactive power flow
        Qflow_f2t=Vf_tps.*(y_line_mag.* ...
                  (Vf_tps.*sin_y_line_ang - Vt.*sin(Vf_delta_ang)) - ...
                          Vf_tps.*line(:,5)./2);
        Qflow_t2f=Vt.*(y_line_mag.* ...
                  (Vt.*sin_y_line_ang - Vf_tps.*sin(Vt_delta_ang)) - ...
                       Vt.*line(:,5)./2);
         
        % I. Consolidate the Real and Reactive power flows
        Pflow=[Pflow_f2t Pflow_t2f];
        Qflow=[Qflow_f2t Qflow_t2f];
	      	   	  
        % J. Compute the active power flow
        Sflow=max(sqrt(Pflow.^2+Qflow.^2)')';
	      	   	  
        % K. Prepare the matrix of all lines and flows
        line_flow=[line(:,1:2) Pflow Qflow Sflow];

        % M. Convert the bus angle to degrees
        bus(:,3)=bus(:,3)*180/pi;
        % N. Convert the line angle to degrees
        line(:,7)=line(:,7)*180/pi;
	      	   	  
        % Loadflow complete - return
        return; % Only valid exit point
      end
	      	   	  
      % E. Perform LU Decomposition & verify Jacobian condition
      %    (This uses an estimate of the reciprocal condition number)
      [l,u,p,q] = lu(jac);
      % [l,u,p] = lu(jac); %% BSIEVER: For Matlab Versions prior to 7
      condition=full(min(abs(diag(u)))/max(abs(diag(u))));
      if condition<min_condition
        error('loadflow:ill_condition',sprintf('Ill conditioned Jacobian (Condition Number: %f)',condition));
      end
	      	   	  
      % F. Solve for the deltaAngle and deltaV
      deltaAngdeltaV=q*(u \ (l \ (p*[Pmismatch(angleMap); Qmismatch(vMap)])));
      % deltaAngdeltaV=(u \ (l \ (p*[Pmismatch(angleMap); Qmismatch(vMap)]))); %% BSIEVER: For Matlab Versions prior to 7
      % G. Update the angle and V estimates
      bus(angleMap,3)=bus(angleMap,3) + ...
                      deltaAngdeltaV(1:length(angleMap));
      % Note that voltages are updates using the (dP/dV)V and (dQ/dQ)*V
      % to unify the components of the Jacobian (and reduce computation some)
      % The calculated value is delta V / V, 
      % hence the multiplication by V to find the delta See [3]
      bus(vMap,2)=bus(vMap,2) + ...
                  bus(vMap,2).*deltaAngdeltaV(length(angleMap)+1:end);

    end
	      	   	  
    % V. Error - Maximum number of iterations reached
    error('loadflow:max_iterations','Max Iterations Exceeded')
    return;
	      	   	  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [j, deltaP, deltaQ]=jmat(bus, y_mag, y_ang)
% Form the jacobian matrix from the  bus, and admittance matrices
% and the power mismatch partial equations (the last component of mismatch)
% deltaP is the adjustment to real power (excluding load and generation)
% deltaQ is the adjustment to reactive power (excluding load and generation)
	      	   	  
  % I. Create the delta matrix ( v_ang(i)-v_vang(j)-y_theta(i,j) )
    % A. For each non-zero i,j, get the corresponding delta(i)-delta(j)
    [m, n] = find(y_mag~=0|y_ang~=0);
    delta = bus(m,3)-bus(n,3);
    % B. Sum them: delta(i)-delta(j)-theta(i,j)
    ang = sparse(m,n,delta)-y_ang;
    
    % C. Compute the cos of the angle
    cos_ang = spfun(@cos,ang);
    sin_ang = spfun(@sin,ang);
	      	   	  
  % II. Compute Common Terms:
    % A. Compute |y|*cos(delta(i)-delta(j)-theta(i,j))
    y_cos = y_mag.*cos_ang;
    % B. Compute |y|*sin(delta(i)-delta(j)-theta(i,j))
    y_sin = y_mag.*sin_ang;
	      	   	  
    % C. Get the diagonalized voltages
    busSz = size(bus,1);
    diag_v = spdiags(bus(:,2),0,busSz,busSz);
	      	   	  
    % D. Compute |v|*|y|*cos(delta(i)-delta(j)-theta(i,j))
    v_y_cos_v = diag_v*y_cos*diag_v;
    % E. Compute |v|*|y|*sin(delta(i)-delta(j)-theta(i,j))
    v_y_sin_v = diag_v*y_sin*diag_v;
	      	   	  
    % F. Compute |v|*(|y|*cos(delta(i)-delta(j)-theta(i,j)))*|v|) 
    %    This is the P mismatch delta
    deltaP = bus(:,2).*(y_cos*bus(:,2));
    % G. Compute |v|*(|y|*sin(delta(i)-delta(j)-theta(i,j)))*|v|)
    %    This is the Q mismatch delta
    deltaQ = bus(:,2).*(y_sin*bus(:,2));
	      	   	  
    % H. Generate the Diagonal sparse versions:
    deltaP_diag = spdiags(deltaP,0,busSz,busSz);
    deltaQ_diag = spdiags(deltaQ,0,busSz,busSz);
	      	   	  
  % III. Compute the components of the Jacobian [j11 j12; j21 j22]   
    % dP/ddelta = Non-diagonal terms and diagonal terms
    j11 =  v_y_sin_v-deltaQ_diag; 
    % dQ/ddelta = Non-diagonal terms and diagonal terms
    j21 = -v_y_cos_v+deltaP_diag; 
    % dP/dV*V = Non-diagonal terms and diagonal terms
    % (based on simplification in [3] - calculate delta V / V
    %  instead to unify the elements of the Jacobian)
    j12 =  v_y_cos_v+deltaP_diag;
    % dQ/dV*V = Non-diagonal terms and diagonal terms
    % (based on simplification in [3] - calculate delta V / V
    %  instead to unify the elements of the Jacobian)
    j22 =  v_y_sin_v+deltaQ_diag;
	      	   	  
  % VI. Remove the swing bus from all sub-matrices and
  %     the PV buses from corresponding matrices
    % A. Find the swing (reference) buses
    swings=find(bus(:,10)==1);
    % B. Find everything but the load buses (The PVs and the swings)
    non_load=find(bus(:,10)~=3);  
	      	   	  
    % C. Remove swing buses from dP/ddelta
    j11(swings,:)=[];
    j11(:,swings)=[];
    % D. Remove the Generator buses and swing from dQ/ddelta:
    j21(non_load,:)=[];
    j21(:,swings)=[];
    % E. Remove the Generator buses and swing from dP/dV:
    j12(:,non_load)=[];
    j12(swings,:)=[];
    % F. Remove the Generator buses and swing from dQ/dV:
    j22(:,non_load)=[];
    j22(non_load,:)=[];
	      	   	  
  % V. Return the jacobian:
    j=[j11 j12; j21 j22];  
    return;
