%Agnes Optimization 4
%June 19, 2014
%This model is going to account for all the individual chillers.  It is
%also going to look at the actual wet bulb temperatures and loads on the
%evaporator.

%Notes:
%0 = Off, 1 = On, 2 = Check
%Temperatures coming in F

clc
clear all
close all
fileName = 'Y:\Agnes\From Mitch Paulus\rawDataNumberFormatted.csv'
dTfileName = 'Y:\Agnes\From Mitch Paulus\dT model.csv'
RawData = csvread(fileName);

dTdata = csvread(dTfileName);

%DeltaT = -0.304*Tevap + 3.41E-5 * PlantLoad + 0.0608*Twb + 19.5286
DeltaT = -0.304*(44) + 3.41E-5.*dTdata(:,1) + 0.0608.*dTdata(:,3) + 19.5286;

calcTwb = -9.4644E-20.*dTdata(:,1).^4 + 5.5069E-14.*dTdata(:,1).^3 -1.2957E-8.*dTdata(:,1).^2 + 1.5223E-3 .* dTdata(:,1) - 0.8959;

DeltaTwithCalcTwb = -0.304*(44) + 3.41E-5.*dTdata(:,1) + 0.0608.*calcTwb + 19.5286;

%Flow (gpm) = PlantLoad (kBTU/hr) * 1000 / 500 / DeltaT
PlantFlow = dTdata(:,1) .* 1000 ./ 500 ./ DeltaT;



%Loop through to only get times when there are no check values
j =1;
for i = 1:length(RawData)
   if RawData(i,4) ~= 2 && RawData(i,6) ~=2 && RawData(i,8) ~= 2 && RawData(i,10)~=2 && RawData(i,12)~=2
        Tdb(j) = RawData(i,2);
        Twb(j) = RawData(i,3);
        Qevap(j,1) = RawData(i,5);
        Qevap(j,2) = RawData(i,7);
        Qevap(j,3) = RawData(i,9);
        Qevap(j,4) = RawData(i,11);
        Qevap(j,5) = RawData(i,13);
        OnOff(j,1) = RawData(i,4);
        OnOff(j,2) = RawData(i,6);
        OnOff(j,3) = RawData(i,8);
        OnOff(j,4) = RawData(i,10);
        OnOff(j,5) = RawData(i,12);
        ActualPower(j,1) = RawData(i,14);
        ActualPower(j,2) = RawData(i,15);
        ActualPower(j,3) = RawData(i,16);
        ActualPower(j,4) = RawData(i,17);
        ActualPower(j,5) = RawData(i,18);   
        j = j + 1;
   end        
end


ToutEvap = 36:44;              %Degrees F (leaving evaporater)
Vcond = 50000:5000:100000;            %Conderser flow rate (ft^3/hr)
power = zeros(length(Qevap),length(ToutEvap),length(Vcond),5);
TinCond = zeros(length(Qevap),1);

for i = 1:length(Qevap)   
    for j = 1:length(ToutEvap)
       for k = 1:length(Vcond)
           if Twb(i) < 55
               TinCond(i) = 60+460;
           else
               TinCond(i) = Twb(i) + 5 + 460;
           end
           
           for chiller = 1:5
              if OnOff(i,chiller) == 0
                  power(i,j,k,chiller) = 0;
              else
                  power(i,j,k,chiller) = chillerModel1(Qevap(i,chiller),ToutEvap(j)+460,TinCond(i),Vcond(k),chiller);
              end                                                            
           end           
       end                                
    end                                            
end

%find the conditions for the minimum energy value
minPower = zeros(length(Qevap),3);
for i = 1:length(Qevap)   
    minPower(i,1) = 999999;
   for j = 1:length(ToutEvap)
      for k = 1:length(Vcond)                     
        if power(i,j,k,5) < minPower(i,1) && OnOff(i,5) == 1
            minPower(i,1) = power(i,j,k,5);
            minPower(i,2) = ToutEvap(j);
            minPower(i,3) = Vcond(k);                       
        end                                                       
      end
   end                                
end





%Calculate pumping power
pumpPowerPred = calcPumpPower(PlantFlow);
pumpPowerAct = calcPumpPower(dTdata(:,4));

pumpSavings = pumpPowerAct - pumpPowerPred;
pumpSavingsPercent = pumpSavings ./ pumpPowerAct .* 100;

%calculate the total chiller savings
totalChillerSavings = zeros(length(Qevap),1);
totalChillerLoad = zeros(length(Qevap),1);
totalActualChillerPower = zeros(length(Qevap),1);
totalChillerSavingsPercent = zeros(length(Qevap),1);
for i = 1:length(OnOff)
    for chillerNum = 1:5
        if OnOff(i,chillerNum) == 1
            totalChillerSavings(i) = totalChillerSavings(i) + ActualPower(i,chillerNum) - power(i,9,11,chillerNum);
            totalActualChillerPower(i) = totalActualChillerPower(i) + ActualPower(i,chillerNum);
            totalChillerLoad(i) = totalChillerLoad(i) + Qevap(i,chillerNum);
        end                
    end
    totalChillerSavingsPercent(i) = totalChillerSavings(i) / totalActualChillerPower(i) * 100;    
end





save('Y:\Agnes\From Mitch Paulus\PowerData.mat')

