% PPC  Probabalistic production cost simulation for small system.
%    [ETOC,LOLP,EENS] = PPC(MU,SIGMA,HMAX,HAVAIL,WMAX,WAVAIL,GMAX,GAVAIL,GCOST,SEG)
%    simulates a system with a run-of-the-river hydro power plant
%    (installed capacity HMAX, availability HAVAIL), a wind farm 
%    (installed capacity WMAX, available generation proabilities WAVAIL)
%    and) a thermal power plant (installed capacity GMAX, availability 
%    GAVAIL, operation cost GCOST). The load is normally distributed and
%    has the mean MU and the standard deviation SIGMA. The segment length 
%    when approximating the load duration curve is SEG.
%
%    For hydro and thermal power plants the availability is the
%    probability that the power plant can generate its installed
%    capacity, i.e., a number between 0 and 1. For the wind power plant
%    the availability WAVAIL should be a vector. The first element of
%    WAVAIL states the probability that no wind power generation is
%    available at all, the second element states the probability that 50% 
%    of the installed capacity is avaialable and the last element states
%    the probability that the whole installed capacity is available. The
%    sum of the elements of WAVAIL should be equal to 1.

function [etoc,lolp,eens] = ppc(muD,sigmaD,Hmax,Havail,Wmax,Wavail,...
				Gmax,Gavail,Gcost,seg)

currseg = 1; prevval = 1;
durcurveval = 1 - Phi((currseg*seg-muD)/sigmaD); 
while durcurveval > 1e-4
  F0(1,currseg) = (prevval+durcurveval)/2;
  currseg = currseg + 1;
  prevval = durcurveval;
  durcurveval = 1 - Phi((currseg*seg-muD)/sigmaD); 
end
if Havail < 1
  Hqoffset = round(Hmax/seg); 
  F1 = Havail*[F0 zeros(1,Hqoffset)] + (1-Havail)*[ones(1,Hqoffset) F0];
else
  Hqoffset = round(Hmax/seg);
  F1 = F0;
end
if Wavail(1) < 1
  Wqoffset = round(Wmax/seg/2);
  F2 = Wavail(3)*[F1 zeros(1,2*Wqoffset)] + ...
       Wavail(2)*[ones(1,Wqoffset) F1 zeros(1,Wqoffset)] + ...
       Wavail(1)*[ones(1,2*Wqoffset) F1];
  Wqoffset = round(Wmax/seg);
else
  F2 = F1; Wqoffset = 0;
end
Gqoffset = round(Gmax/seg);
F3 = Gavail*[F2 zeros(1,Gqoffset)] + (1-Gavail)*[ones(1,Gqoffset) F2];
lolp = 100*F3(Hqoffset+Wqoffset+Gqoffset+1);
EENS2 = sum(F2((Hqoffset+Wqoffset+1):length(F2)))*seg;
EENS3 = sum(F3((Hqoffset+Wqoffset+Gqoffset+1):length(F3)))*seg;
etoc = Gcost*(EENS2-EENS3);
eens = EENS3;