Free cookie consent management tool by TermsFeed Policy Generator

source: misc/publications/2012/journals/IEEE_TII/EvalSolutionTC4.m @ 13390

Last change on this file since 13390 was 9157, checked in by shuttere, 12 years ago
File size: 7.4 KB
Line 
1errorBranch=zeros(108,24);
2numberSamples=Samples;
3evaluatedSamples=Samples;
4penaltyFactors=[5,5,1,0.1];
5costPenaltyGrid=zeros(numberSamples,24);
6costPenaltyTotal=zeros(numberSamples,1);
7cost=zeros(numberSamples,24);
8costTemp=zeros(3,1);
9resultingCost=zeros(numberSamples,1);
10resultingCostLog=zeros(numberSamples,1);
11
12%%sample ziehen und bewerten
13for s=1:numberSamples   
14   
15    %Schedules mit Zeitfenster übelagern, für jeden Zeitschritt PF lösen
16    n=0; %spaltenindex in der Kostenmatrix, Zeitstempel im Schedule   
17    for j=1:24
18        casexx=eval('ieeecase118');   
19 
20        n=n+1;
21         
22            r = 1 + 0.04.*randn(1,1);
23            %Lasten ins Stromnetz einfügen
24            %Lastprofile setzen und randomisieren
25            %Busse 1-30 häuslich
26            casexx.bus(1:30,3)=casexx.bus(1:30,3)*scheduleBaseLoad1(j)*r;
27            casexx.bus(1:30,4)=casexx.bus(1:30,4)*scheduleBaseLoad1(j)*r;
28            %Busse 31-60 häuslich/gewerblich
29            casexx.bus(31:60,3)=casexx.bus(31:60,3)*(scheduleBaseLoad2(j)+scheduleBaseLoad4(j))/2*r;
30            casexx.bus(31:60,4)=casexx.bus(31:60,4)*(scheduleBaseLoad2(j)+scheduleBaseLoad4(j))/2*r;
31            %Busse 61-90 häuslich/industriell
32            casexx.bus(61:90,3)=casexx.bus(61:90,3)*(scheduleBaseLoad3(j)+scheduleBaseLoad7(j))/2*r;
33            casexx.bus(61:90,4)=casexx.bus(61:90,4)*(scheduleBaseLoad3(j)+scheduleBaseLoad7(j))/2*r;
34            %Busse 91-100 industriell/gewerblich
35            casexx.bus(91:100,3)=casexx.bus(91:100,3)*(scheduleBaseLoad8(j)+scheduleBaseLoad5(j))/2*r;
36            casexx.bus(91:100,4)=casexx.bus(91:100,4)*(scheduleBaseLoad8(j)+scheduleBaseLoad5(j))/2*r;
37            %Busse 101-110 industriell
38            casexx.bus(101:110,3)=casexx.bus(101:110,3)*scheduleBaseLoad9(j)*r;
39            casexx.bus(101:110,4)=casexx.bus(101:110,4)*scheduleBaseLoad9(j)*r;
40            %Busse 111-118 gewerblich           
41            casexx.bus(111:118,3)=casexx.bus(111:118,3)*scheduleBaseLoad6(j)*r;
42            casexx.bus(111:118,4)=casexx.bus(111:118,4)*scheduleBaseLoad6(j)*r;
43
44            %Erzeugung aus Windkraft setzen
45            %zufällige Windkraft ziehen
46            %%Erzeuger randomisieren
47            %Weibull- Verteilung wie in "probabilistic constrained load flow
48            %considering integration of wind power generation and electric vehicles"
49            %jedoch direkt über wbl- funktion und mit geänderten parametern
50            v=wblrnd(1.11,3)*14*scheduleWindSpeed(j); % gibt im mittel 1; 14m/s=windspeed rated
51            %shape=1.11, scale=3
52            %PW1 aus Leistungskurve berechnen
53            PW=polyval(p1,v);         
54           
55            if(PW<0)
56                PW1=0;
57            end
58            if(PW>5)
59                PW=5;
60            end         
61           
62     
63      %PV u Wind setzen
64            indexBus=1;
65            for i=1:108
66                if(mod(i,2)==0)
67                    casexx.gen(indexBus,2)=casexx.gen(indexBus,2)+(scheduleIrradiance(j)*Solution(i)*0.5);
68                    indexBus=indexBus+1;
69                end
70                if(mod(i,2)==1)
71                    casexx.gen(indexBus,2)=casexx.gen(indexBus,2)+(Solution(i)*PW);                 
72                end               
73            end
74            %PV u. Wind mit halber Leistung
75            indexBus=1;
76            for i=109:216
77                if(mod(i,2)==0)
78                    casexx.gen(indexBus,2)=casexx.gen(indexBus,2)+(scheduleIrradiance(j)*Solution(i)*0.25);
79                    indexBus=indexBus+1;
80                end
81                if(mod(i,2)==1)
82                    casexx.gen(indexBus,2)=casexx.gen(indexBus,2)+(Solution(i)*PW*0.5);                 
83                end               
84            end       
85     
86            %LF berechnen mit Fast Decoupled XB (PF_ALG=2)
87            opt = mpoption('PF_ALG', 2,'VERBOSE',0,'OUT_ALL',0);
88            results=runpf(casexx,opt);
89           
90            %Restriktionen überprüfen
91                % max Leitungsscheinleistung
92                results.branchMVA=sqrt(results.branch(:,14).^2+results.branch(:,15).^2);
93               
94                %anpassen, falls notwendig
95                 k=1;
96                 deltaMVA=0;
97                 while k<=length(MaxBranchMVA)
98                    if results.branchMVA(k,1)>MaxBranchMVA(k,1)
99                       deltaMVA=deltaMVA+(results.branchMVA(k,1)-MaxBranchMVA(k,1))^2;
100                       errorBranch(k,j)=1;
101                       %dlmwrite('errorBranch.txt',errorBranch,'-append');
102                    end
103                    k=k+1;   
104                 end
105
106                % max Generatorblindleistung               
107                 deltaQ=0;
108%                 k=1;
109%                 while k<=size(results.gen,1)
110%                     if results.gen(k,3)>results.gen(k,4)
111%                         deltaQ=deltaQ+(results.gen(k,3)-results.gen(k,4))^2;
112%                     end
113%                     if results.gen(k,3)<results.gen(k,5)
114%                         deltaQ=deltaQ+(results.gen(k,3)-results.gen(k,5))^2;
115%                     end
116%                     k=k+1;
117%                 end
118               
119                % max Spannung
120                k=1;
121                deltaV=0;
122                while k<=size(results.bus,1)
123                    if results.bus(k,8)>results.bus(k,12)
124                        deltaV=deltaV+(results.bus(k,8)-results.bus(k,12))^2;
125                    end
126                    if results.bus(k,8)<results.bus(k,13)
127                        deltaV=deltaV+(results.bus(k,8)-results.bus(k,13))^2;
128                    end
129                    k=k+1;
130                end
131               
132            %Kosten ausdrücken
133                       
134            %Kostenfunktion: min. Verlustleistung
135            sumgenerated=sum(results.gen(:,2));
136            sumload=sum(results.bus(:,3));
137            losses=sumgenerated-sumload;
138           
139            costPenaltyGrid(s,n)=deltaV*penaltyFactors(3)+deltaMVA*penaltyFactors(4);
140            cost(s,n)=losses;     
141    end
142            %Anzahl d. verbauten Kraftwerke soll minimiert werden
143            maxPlants=0;
144            lengthSolution=(sum(Solution)-18);       
145            if(lengthSolution<0)
146                lengthSolution=0;
147            end
148            if(sum(Solution)>27)
149                maxPlants=sum(Solution)-27;
150            end
151           
152            %an jedem Bus entweder Wind oder PV
153            config=0;
154            for i=1:2:108
155                if((Solution(i)+Solution(i+1))>1)
156                config = config +1;
157                end
158            end
159           
160   
161    costPenaltyTotal(s)=penaltyFactors(1)*lengthSolution+penaltyFactors(2)*maxPlants+sum(costPenaltyGrid(s,:));   
162   
163    if(s==4 && Generations>1)
164       
165        for i=1:3   
166        costTemp(i)=sum(cost(i,:))+costPenaltyTotal(i)/Generations;
167        end
168       
169        fitnessTemp=sum(costTemp(:))/size(costTemp,1);
170       
171        if(fitnessTemp>threshold)
172            evaluatedSamples=3;
173            break;
174        end
175    end
176end
177
178%% Kosten ermitteln
179for i=1:evaluatedSamples
180     
181    if(Generations < 1)
182        Generations=1;
183    end
184   
185    resultingCost(i)=sum(cost(i,:))+costPenaltyTotal(i);
186end
187
188%dlmwrite('logSamples.txt',evaluatedSamples,'-append','newline', 'pc');
189
190Fitness=sum(resultingCost(1:evaluatedSamples))/evaluatedSamples;
191Fitness
Note: See TracBrowser for help on using the repository browser.