%% power flow without DG
%power flow case 30 IEEE
clc;
clear;
close all;
bus = 41;
Vdeflecty_wo = zeros (bus,1);  %Voltage deflected yearly
consmaxV = 1.05;
consminV = 0.95;
Ces = 3492;                %cost of sold energy to consumer 0.57*6126=3491.82 (IRR/kWh) 
Ceb = 2141;                %cost of buy energy from main grid 0.35*6126=2144.1 (IRR/kWh)
Ruc = 13.3; %cost of energy not supply
Cenvi = 0.96;% Co2 emission tax rate
B = 0.3;% Co2 emission(Kg/kWh)
Bjmu = 0.7;% sercice of equpment 
Cmar = 13.3 ;%0.52*1000: MWh *6938 = 13.3 dalar
sensfactor = zeros(41,6);           % sense factor for make desition
mpc = loadcase('case30');%load case30 to mpc 
%%load demand
demloadPD = importdata('demandload.xlsx');%load demandload.xlsx form root case30-v1(load curve).
mpc.bus( :,3) = demloadPD.LFP( :,1);
%% power flow without DG
results = runopf(mpc);
% results1 = runopf(mpc);
%% calculat sence factor
sensfactor(:,1:2) = results.branch(:,1:2);% put index brach
sensfactor(:,3) = (abs(real(get_losses(results))));%get power loss in the branch
sensfactor(:,4) = abs(results.branch(:,14));
sensfactor(:,5) = (sensfactor(:,3)./sensfactor(:,4));
[row1] = find(isnan (sensfactor(:,5)));
[row2] = find(isinf (sensfactor(:,5)));
sensfactor(row1,:) = [];% make empty row
sensfactor(row2,:) = [];% make empty row
sensfactor = sortrows(sensfactor,5);        % sort sense factor
locatdg = sensfactor (37:40,2);             % chose 4th of end 
%%cost function
% f1=results1.f;                  % cost of operation and maintenance 
% loss = sum(real(get_losses(results)));
% f3=7*loss;          % cost of losses
%% deflected voltage
Lanbada_wo =sum((results.bus( :,8)>consmaxV | results.bus( :,8)<consminV));
Vy_wo = find (results.bus( :,8)>consmaxV | results.bus( :,8)<consminV);%find deflect voltage for one year
Vdeflecty_wo(Vy_wo) = find (results.bus( :,8)>consmaxV | results.bus( :,8)<consminV);% find deflect voltage 
ENNS_wo =Ruc*((results.gen(4,2)*1)*(mpc.gencost(4,5))...
    +(results.gen(5,2)*2)*(mpc.gencost(5,5))+(results.gen(6,2)*2)*(mpc.gencost(6,5)));%Energy not supply without DG* 3dolar
E_wo = sum(results.gen(:,2));%product GEN1-6
Bupda_wo = sum(1.73*Cmar*( abs( results.branch(:,14))));
Bimpu_wo = Lanbada_wo*Bjmu;
Bevin = E_wo*Cenvi*B;
% Bsub = Ceb(E_wo-E_w)-8760(Cdgi*Pi*CFi); Calculat After dession DG 