% Capacity Model creator
% COPT=cap_model(units, derated_state_units ,load)
%
% units:(number of different units x3 ) (NO DE-RATED STATES SUPPORTED)
% # of units | Gen. capacity(MW) | FOR
%
% load (MW)
% function COPT=cap_model(units_NDS)
units_NDS = [4,20,0.004975124378109;
1,15,0.004975124378109;
7,5, 0.004975124378109];
% 1,60,0.02];
%1,30,0.08];
%1,30,0.08;]
% % status #states states + capacity
% units_DSI = [%1 2 0 (1-0.02) 25 0.02 0 0
% %1 2 0 (1-0.02) 25 0.02 0 0
% 1 3 0 0.86 20 0.06 40 0.08]
nbUnits_NDS=size(units_NDS,1);
number= units_NDS(:,1);
gen= units_NDS(:,2);
pr= units_NDS(:,3);
maxcapacity=sum(number.*gen);
totalpdf=1;
gen=round(gen*1000); %We round to the KWatt. The rest of the calculations use KW!!
interval=gen(1); %Interval in pdf is the maximum common divisor of ge
for i=2:nbUnits_NDS
interval=gcd(interval,gen(i));
end
for i=1:nbUnits_NDS
capOutaged=0:gen(i):number(i)*gen(i);
individualProb=binopdf(0:number(i),number(i),pr(i));
pdf=[];
pdf(capOutaged/interval+1)=individualProb;
totalpdf=conv(totalpdf,pdf);
end
%stem(0:length(totalpdf)-1,totalpdf);grid on;title('Total probability density function');xlabel('Outaged capacity');ylabel('Probability');
[~,outaged,probability]=find(totalpdf);
outaged=(outaged-1)*interval/1000;
nstates=length(probability);
cumulative=zeros(nstates,1);
for i=nstates:-1:1
cumulative(i)=sum(probability(i:nstates));
end
format long;
COPT=[outaged',probability', cumulative];
disp(' Outaged Indv.Prob Cumul.Prob')
disp(COPT);
% nbUnits_DSI=size(units_DSI,1)
% number_DSI= units_DSI(:,1)
% gen_DSI= [units_DSI(:,3) units_DSI(:,5) units_DSI(:,7)]
% pr_DSI= [units_DSI(:,4) units_DSI(:,6) units_DSI(:,8)]
% maxcapacity=maxcapacity + sum(number_DSI.*gen_DSI(:,3))
% Step_DSI=[0 10 20 30 40 50 60 70 80 90 100]
cumulative
% % for j=1:size(Step_DSI,2)
% %
% % Pint = 0;
% %
% % for n=1:3
% %
% % if Step_DSI(j) - units_DSI((2*n)+1) > 40
% % Pint = Pint;
% % elseif Step_DSI(j) - units_DSI((2*n)+1) > 20
% % Pint = Pint + units_DSI(2+2*n)*(cumulative(3));
% % elseif Step_DSI(j) - units_DSI((2*n)+1) > 0
% % Pint = Pint + units_DSI(2+2*n)*(cumulative(2));
% % else
% % Pint = Pint + units_DSI(2+2*n)*(cumulative(1));
% % end
% %
% % end
% %
% % P(j) = Pint;
% % end
%
% for i=0:10
% P(i+1,1) = i*10;
% end
%
% P(1,2) = units_DSI(4)*cumulative(1) + units_DSI(6)*(1) + units_DSI(8)*(1);
% P(2,2) = units_DSI(4)*cumulative(2) + units_DSI(6)*(1) + units_DSI(8)*(1);
% P(3,2) = units_DSI(4)*cumulative(3) + units_DSI(6)*cumulative(1) + units_DSI(8)*cumulative(1);
% P(4,2) = units_DSI(4)*cumulative(4) + units_DSI(6)*cumulative(2) + units_DSI(8)*cumulative(1);
% P(5,2) = units_DSI(4)*cumulative(5) + units_DSI(6)*cumulative(3) + units_DSI(8)*cumulative(1);
% P(6,2) = units_DSI(4)*cumulative(6) + units_DSI(6)*cumulative(4) + units_DSI(8)*cumulative(2);
% P(7,2) = units_DSI(4)*cumulative(7) + units_DSI(6)*cumulative(5) + units_DSI(8)*cumulative(3);
% P(8,2) = units_DSI(4)*(0) + units_DSI(6)*cumulative(6) + units_DSI(8)*cumulative(4);
% P(9,2) = units_DSI(4)*(0) + units_DSI(6)*cumulative(7) + units_DSI(8)*cumulative(5);
% P(10,2) = units_DSI(4)*(0) + units_DSI(6)*(0) + units_DSI(8)*cumulative(6);
% P(11,2) = units_DSI(4)*(0) + units_DSI(6)*(0) + units_DSI(8)*cumulative(7);
%
% P
% end
