%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  encost.m - Uses supplied 24-hour load profile to calculate the
%             loss energy cost to operate the load profile at
%             rated output voltage. Also calculates the average
%             or all day efficiency.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
R1=0.20; R2=0.002; X1=0.60; X2=0.006; Xm=130; Rc=4000;
V1=2400; V2=240; KVA=15; n=V1/V2; R2p=n^2*R2; X2p=n^2*X2;
KWHrcost=0.10;                   % $/KWHr cost
LOAD=[55 75 100 90 65 55];       % Percent apparent power load
PF=[0.8 0.85 0.9 0.85 0.77 0.8]; % Power factor, lagging assumed
HRS=[6 6 4.5 3.5 2 2];           % Hours at load point
n=length(HRS); T(1)=HRS(1); for i=2:n; T(i)=T(i-1)+HRS(i); end
if T ~= 24; disp('WARNING - 24 hour period not specified'); end
for i=1:n
   I2=KVA*1000/V2*LOAD(i)/100/PF(i)*exp(-j*acos(PF(i)));
   E1=(R2p+j*X2p)*I2/n+n*V2;
   Im=E1/(j*Xm); Ic=E1/Rc; Io=Ic+Im; I1=I2/n+Io;
   V1=(R1+j*X1)*I1+E1;
   thet1=angle(V1)-angle(I1); PFin=cos(thet1);
   Pout=real(V2*conj(I2)); Pin=real(V1*conj(I1));
   eff(i)=Pout/Pin*100; Ploss(i)=Pin-Pout; Po(i)=Pout/1000;
end
cost=0; adeff=0;
for i=1:n
   cost=cost+Ploss(i)/1000*KWHrcost*HRS(i);
   adeff=adeff+eff(i)*HRS(i)/T(n);
end
subplot(2,1,1); stairs(T,Po); grid
title('Power load profile');
ylabel('Load, kW'); xlabel('Time, hr');
subplot(2,1,2); stairs(T,eff); grid
title('Effieiency profile');
ylabel('Efficiency, %'); xlabel('Time, hr');
disp('  ALL DAY EFFICIENCY STUDY'); disp(' ')
disp([' 24-Hour Energy cost -  $', num2str(cost)]); disp(' ');
disp([' All-Day Efficiency(%) - ', num2str(adeff)]);