% Delivered to the U.S. Government in Accordance with Contract FA8650-10-C-7080.
% 
% (c) 2011 United Technologies Corp.
% Developed with the sponsorship of the Defense Advanced Research Projects Agency (DARPA).
% Permission is hereby granted, free of charge, to any person obtaining a copy of this data,
% including any software or models in source or binary form, as well as any drawings, 
% specifications, and documentation (collectively "the Data"), to deal in the Data without 
% restriction, including without limitation the rights to use, copy, modify, merge, publish,
% distribute, sublicense, and/or sell copies of the Data, and to permit persons to whom the Data
% is furnished to do so, subject to the following conditions:
%  
% The above copyright notice and this permission notice shall be included in all copies or
% substantial portions of the Data.
%  
% THE DATA IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING
% BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
% NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS, SPONSORS, DEVELOPERS, CONTRIBUTORS, OR
% COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
% OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE DATA OR THE
% USE OR OTHER DEALINGS IN THE DATA.

% Authors:  Fernando Rodriguez                      John Finn
%           Sr Engineer/Scientist                   Contractor (Summer Intern 2011)
%           Power Electronics Systems               Embedded Systems & Networks
%           United Technologies Research Center     United Technologies Research Center
%           411 Silver Lane                         411 Silver Lane
%           East Hartford,CT 06108                  East Hartford,CT 06108
%           rodrigf@utrc.utc.com                    ccfinnje@utrc.utc.com
%           tel: (860) 610-7102                     tel: +1 860-610-7696
%
%% EPS_Load_Flow.m
% This file executes load flow analysis of an EPS given a Netlist.
clear all; close all; clc
% Set Newton Raphson solver parameters
accuracy = 0.0001; maxiter = 10;

%% ----------------------------Load NETLIST--------------------------------
[Source_data, TL_data, SW_data, Load_data, Conv_data, EPS_Code] = topologyRev03_36L;
% find total number of buses "max_bus"
max_bus1 = max([Source_data(:,1); Load_data(:,1)]);
max_bus2 = max([TL_data(:,3:4); SW_data(:,3:4); Conv_data(:,2:3)]);
max_bus = max([max_bus1; max_bus2(1); max_bus2(2)]);
clear max_bus1 max_bus2
%
%% -------Graph Reduction: Eliminate components not power by source--------
[LNK_data, Source_data, Load_data, Conv_data, EPS_Code]= Graph_reduction...
    (Source_data, TL_data, SW_data, Load_data, Conv_data, EPS_Code);
%
%% -------Bus re-index: renumbers "n" powered buses from 1 thru n----------
[Source_data, TL_data, SW_data, Load_data, Conv_data, EPS_Code, Mapping]...
    = reorder(LNK_data, Source_data, Load_data, Conv_data, EPS_Code);
%
%% ---------EPSdata_process: Create EPS subsytems for analysis-------------
[ALL_linedata, ALL_busdata, n_min, n_max, a, subs, EPS_Code, Conv_data, n_cd]...
    = EPSdata_process(Source_data, TL_data, SW_data, Load_data, Conv_data, EPS_Code);
%
%% Initialize Output variables to zero
nbus = length(ALL_busdata(:,1));
Vm_=zeros(1,nbus); deltad_=Vm_; Pd_=Vm_; Qd_=Vm_; Pg_=Vm_; Qg_=Vm_;
Pdt_=0; Qdt_=0; Pgt_=0; Qgt_=0;
%
%% Evaluate EPS subsystems
i=subs; ii=1; iii=1; iv=1;
while i ~=0
    r=sum(a(1:i)-1);
    q=r-a(i)+2;
    linedata=ALL_linedata(q:r,2:4); % Defines line data for subsystem i
    if i < subs;
        for ii=1:length(Conv_data(:,1))
            if Conv_data(ii,5) == i+1
                iii=Conv_data(ii,2);
                iv=Conv_data(ii,3);
                eff = Conv_data(ii,6);
                ALL_busdata(iii,5:6) = [Pg_(iv)/eff 0];
            end
        end
    end
    busdata = ALL_busdata(n_min(i):n_max(i),:); % Defines code i busdata
    for r=1:length(busdata(:,1))   % Re-index busdata and set Vm
        busdata(r,1)=busdata(r,1)-n_min(i)+1;
        busdata(r,3) = EPS_Code(2,i);
    end
    vh = [i ii iii iv];                     % variable hold
    [Ybus, nbr, nl, nr] = lfybus(linedata); % Create Ybus
    %lfnewton(Ybus, nbr, nl, nr, busdata, linedata)   % Run load flow solver
    lfnewton
    busoutprep                              % Save results
    i=vh(1);ii=vh(2);iii=vh(3); iv=vh(4);   % restore variable values
    i=i-1;
end

% Subtract Converter Power and Load used for intermidiate solution
bus_in = Conv_data(:,2);
bus_out = Conv_data(:,3);
if n_cd ~= 0
    for i = 1:length(Conv_data(:,1))
        Pg_(bus_out(i)) = 0;
        Qg_(bus_out(i)) = 0;
        Pd_(bus_in(i)) = 0;
        Qd_(bus_in(i)) = 0;
    end
end
Pdt_=sum(Pd_);
Qdt_=sum(Qd_);
Pgt_=sum(Pg_);
Qgt_=sum(Qg_);
busout           % Display calculated bus values
I_ = (conj(S_)./conj(V_))';

%% ----Output structure for Power Flow reults in busout.m-----
%         Bus Bus  Vltg  Angle  -----Load-----   -Generator- 
%         No  code Mag   Degree   watt    kvar    watt     var
% busdata=[1   1    Vm   0.0      0.00    0.00    0.00    0.00
%
