errorBranch=zeros(108,24);
numberSamples=Samples;
evaluatedSamples=Samples;
penaltyFactors=[5,5,1,0.1];
costPenaltyGrid=zeros(numberSamples,24);
costPenaltyTotal=zeros(numberSamples,1);
cost=zeros(numberSamples,24);
costTemp=zeros(3,1);
resultingCost=zeros(numberSamples,1);
resultingCostLog=zeros(numberSamples,1);

%%sample ziehen und bewerten
for s=1:numberSamples    
    
    %Schedules mit Zeitfenster übelagern, für jeden Zeitschritt PF lösen
    n=0; %spaltenindex in der Kostenmatrix, Zeitstempel im Schedule    
    for j=1:24
        casexx=eval('ieeecase118');    
  
        n=n+1;
         
            r = 1 + 0.04.*randn(1,1);
            %Lasten ins Stromnetz einfügen
            %Lastprofile setzen und randomisieren
            %Busse 1-30 häuslich
            casexx.bus(1:30,3)=casexx.bus(1:30,3)*scheduleBaseLoad1(j)*r;
            casexx.bus(1:30,4)=casexx.bus(1:30,4)*scheduleBaseLoad1(j)*r;
            %Busse 31-60 häuslich/gewerblich
            casexx.bus(31:60,3)=casexx.bus(31:60,3)*(scheduleBaseLoad2(j)+scheduleBaseLoad4(j))/2*r;
            casexx.bus(31:60,4)=casexx.bus(31:60,4)*(scheduleBaseLoad2(j)+scheduleBaseLoad4(j))/2*r;
            %Busse 61-90 häuslich/industriell
            casexx.bus(61:90,3)=casexx.bus(61:90,3)*(scheduleBaseLoad3(j)+scheduleBaseLoad7(j))/2*r;
            casexx.bus(61:90,4)=casexx.bus(61:90,4)*(scheduleBaseLoad3(j)+scheduleBaseLoad7(j))/2*r;
            %Busse 91-100 industriell/gewerblich
            casexx.bus(91:100,3)=casexx.bus(91:100,3)*(scheduleBaseLoad8(j)+scheduleBaseLoad5(j))/2*r;
            casexx.bus(91:100,4)=casexx.bus(91:100,4)*(scheduleBaseLoad8(j)+scheduleBaseLoad5(j))/2*r;
            %Busse 101-110 industriell
            casexx.bus(101:110,3)=casexx.bus(101:110,3)*scheduleBaseLoad9(j)*r;
            casexx.bus(101:110,4)=casexx.bus(101:110,4)*scheduleBaseLoad9(j)*r;
            %Busse 111-118 gewerblich           
            casexx.bus(111:118,3)=casexx.bus(111:118,3)*scheduleBaseLoad6(j)*r;
            casexx.bus(111:118,4)=casexx.bus(111:118,4)*scheduleBaseLoad6(j)*r;

            %Erzeugung aus Windkraft setzen
            %zufällige Windkraft ziehen
            %%Erzeuger randomisieren
            %Weibull- Verteilung wie in "probabilistic constrained load flow
            %considering integration of wind power generation and electric vehicles"
            %jedoch direkt über wbl- funktion und mit geänderten parametern
            v=wblrnd(1.11,3)*14*scheduleWindSpeed(j); % gibt im mittel 1; 14m/s=windspeed rated
            %shape=1.11, scale=3
            %PW1 aus Leistungskurve berechnen
            PW=polyval(p1,v);          
            
            if(PW<0)
                PW1=0;
            end
            if(PW>5)
                PW=5;
            end          
            
     
			%PV u Wind setzen
            indexBus=1;
            for i=1:108
                if(mod(i,2)==0)
                    casexx.gen(indexBus,2)=casexx.gen(indexBus,2)+(scheduleIrradiance(j)*Solution(i)*0.5);
                    indexBus=indexBus+1;
                end
                if(mod(i,2)==1)
                    casexx.gen(indexBus,2)=casexx.gen(indexBus,2)+(Solution(i)*PW);                  
                end                
            end
            %PV u. Wind mit halber Leistung
            indexBus=1;
            for i=109:216
                if(mod(i,2)==0)
                    casexx.gen(indexBus,2)=casexx.gen(indexBus,2)+(scheduleIrradiance(j)*Solution(i)*0.25);
                    indexBus=indexBus+1;
                end
                if(mod(i,2)==1)
                    casexx.gen(indexBus,2)=casexx.gen(indexBus,2)+(Solution(i)*PW*0.5);                  
                end                
            end 			
     
            %LF berechnen mit Fast Decoupled XB (PF_ALG=2)
            opt = mpoption('PF_ALG', 2,'VERBOSE',0,'OUT_ALL',0);
            results=runpf(casexx,opt);
            
            %Restriktionen überprüfen
                % max Leitungsscheinleistung
                results.branchMVA=sqrt(results.branch(:,14).^2+results.branch(:,15).^2);
                
                %anpassen, falls notwendig
                 k=1;
                 deltaMVA=0;
                 while k<=length(MaxBranchMVA)
                    if results.branchMVA(k,1)>MaxBranchMVA(k,1)
                       deltaMVA=deltaMVA+(results.branchMVA(k,1)-MaxBranchMVA(k,1))^2;
                       errorBranch(k,j)=1;
                       %dlmwrite('errorBranch.txt',errorBranch,'-append');
                    end 
                    k=k+1;   
                 end

                % max Generatorblindleistung                
                 deltaQ=0;
%                 k=1;
%                 while k<=size(results.gen,1)
%                     if results.gen(k,3)>results.gen(k,4)
%                         deltaQ=deltaQ+(results.gen(k,3)-results.gen(k,4))^2;
%                     end
%                     if results.gen(k,3)<results.gen(k,5)
%                         deltaQ=deltaQ+(results.gen(k,3)-results.gen(k,5))^2;
%                     end
%                     k=k+1;
%                 end
                
                % max Spannung
                k=1;
                deltaV=0;
                while k<=size(results.bus,1)
                    if results.bus(k,8)>results.bus(k,12) 
                        deltaV=deltaV+(results.bus(k,8)-results.bus(k,12))^2;
                    end
                    if results.bus(k,8)<results.bus(k,13)
                        deltaV=deltaV+(results.bus(k,8)-results.bus(k,13))^2;
                    end
                    k=k+1;
                end
                
            %Kosten ausdrücken
                        
            %Kostenfunktion: min. Verlustleistung
            sumgenerated=sum(results.gen(:,2));
            sumload=sum(results.bus(:,3));
            losses=sumgenerated-sumload;
            
            costPenaltyGrid(s,n)=deltaV*penaltyFactors(3)+deltaMVA*penaltyFactors(4);
            cost(s,n)=losses;      
    end
            %Anzahl d. verbauten Kraftwerke soll minimiert werden
            maxPlants=0;
            lengthSolution=(sum(Solution)-18);       
            if(lengthSolution<0)
                lengthSolution=0;
            end
            if(sum(Solution)>27)
                maxPlants=sum(Solution)-27;
            end
            
            %an jedem Bus entweder Wind oder PV
            config=0;
            for i=1:2:108
                if((Solution(i)+Solution(i+1))>1)
                config = config +1;
                end
            end
            
    
    costPenaltyTotal(s)=penaltyFactors(1)*lengthSolution+penaltyFactors(2)*maxPlants+sum(costPenaltyGrid(s,:));   
    
    if(s==4 && Generations>1)
        
        for i=1:3    
        costTemp(i)=sum(cost(i,:))+costPenaltyTotal(i)/Generations;
        end
        
        fitnessTemp=sum(costTemp(:))/size(costTemp,1);
        
        if(fitnessTemp>threshold)
            evaluatedSamples=3;
            break;
        end
    end
end

%% Kosten ermitteln
for i=1:evaluatedSamples
     
    if(Generations < 1)
        Generations=1;
    end
    
    resultingCost(i)=sum(cost(i,:))+costPenaltyTotal(i);
end

%dlmwrite('logSamples.txt',evaluatedSamples,'-append','newline', 'pc');

Fitness=sum(resultingCost(1:evaluatedSamples))/evaluatedSamples;
Fitness