Hi Matpower community.
please run tlbo.m
this is *meta heuristic algorithm (teacher-learning-based-Optimization)*
*and in matpower 7.1b *
*i study on case37 and DGs on the network .*
*i have question! please!*
*can you help me !why after 3 finally 5 iteration optimization be end ?!!*
*thanks a lot*
clc;
clear;
close all;
%% Problem Definition
% Cost Function
CostFunction = @(x) costfunction(x);
nVar = 5; % Number of Unknown Variables
VarSize = [1 nVar]; % Unknown Variables Matrix Size
VarMin = [0, 0, 0, 0, 0]; % Unknown Variables Lower Bound
VarMax = [15, 1, 1, 1, 1]; % Unknown Variables Upper Bound
%% TLBO Parameters
MaxIt = 10; % Maximum Number of Iterations
nPop = 50; % Population Size
%% Initialization
% Empty Structure for Individuals
empty_individual.Position = [];
empty_individual.Cost = [];
% Initialize Population Array
pop = repmat(empty_individual, nPop, 1);
% Initialize Best Solution
BestSol.Cost = inf;
% Initialize Population Members
for i=1:nPop
pop(i).Position = unifrnd(VarMin, VarMax, VarSize);
pop(i).Cost = CostFunction(pop(i).Position);
if pop(i).Cost < BestSol.Cost
BestSol = pop(i);
end
end
% Initialize Best Cost Record
BestCosts = zeros(MaxIt,1);
%% TLBO Main Loop
for it=1:MaxIt
% Calculate Population Mean
Mean = 0;
for i=1:nPop
Mean = Mean + pop(i).Position;
end
Mean = Mean/nPop;
% Select Teacher
Teacher = pop(1);
for i=2:nPop
if pop(i).Cost < Teacher.Cost
Teacher = pop(i);
end
end
% Teacher Phase
for i=1:nPop
% Create Empty Solution
newsol = empty_individual;
% Teaching Factor
TF = randi([1 2]);
% Teaching (moving towards teacher)
newsol.Position = pop(i).Position ...
+ rand(VarSize).*(Teacher.Position - TF*Mean);
% Clipping
newsol.Position = max(newsol.Position, VarMin);
newsol.Position = min(newsol.Position, VarMax);
% Evaluation
newsol.Cost = CostFunction(newsol.Position);
% Comparision
if newsol.Cost<pop(i).Cost
pop(i) = newsol;
if pop(i).Cost < BestSol.Cost
BestSol = pop(i);
end
end
end
% Learner Phase
for i=1:nPop
A = 1:nPop;
A(i)=[];
j = A(randi(nPop-1));
Step = pop(i).Position - pop(j).Position;
if pop(j).Cost < pop(i).Cost
Step = -Step;
end
% Create Empty Solution
newsol = empty_individual;
% Teaching (moving towards teacher)
newsol.Position = pop(i).Position + rand(VarSize).*Step;
% Clipping
newsol.Position = max(newsol.Position, VarMin);
newsol.Position = min(newsol.Position, VarMax);
% Evaluation
newsol.Cost = CostFunction(newsol.Position);
% Comparision
if newsol.Cost<pop(i).Cost
pop(i) = newsol;
if pop(i).Cost < BestSol.Cost
BestSol = pop(i);
end
end
end
% Store Record for Current Iteration
BestCosts(it) = BestSol.Cost;
% Show Iteration Information
disp(['Iteration ' num2str(it) ': Best Cost = ' num2str(BestCosts(it))]);
end
%% Results
% run withDG
figure;
plot(BestCosts, 'LineWidth', 2);
% semilogy(BestCosts, 'LineWidth', 2);
xlabel('Iteration');
ylabel('Best Cost');
grid on;
hold on;
%% power flow with DG
%power flow case 37 IEEE
clc;
close all;
bus = 37;
Vdeflecty_wo = zeros (bus,1); %Voltage deflected yearly
%constraint
consmaxV = 1.05;
consminV = 0.95;
%constant
Ces = 490000; %cost of sold energy to consumer
490IRR=490000(IRR/MWh)
Ceb = 245000; %cost of buy energy from main grid
245IRR=245000(IRR/MWh)
Ruc = 326000; %cost of energy not supply
326IRR=326000(IRR/MWh)
Cenvi = 245000; % Co2 emission tax rate
Bco2 = 163000; % Co2 emission(Kg/MWh)
Bjmu = 10200000; %5000*0.15= 750*13600=10200000 sercice of
equpment
Cmar = 29400000 ; %490000*60 = 24900000
RGp = 264000; %price purchasing energy from the main grid
Ci = 735000; % investement cost
Ciop = 0.1; %operation and maintenance
sensfactor = zeros(36,6); % sense factor for make desition
mpc = loadcase('w_dg_case37'); %load case37 to mpc
%%load demand
% demloadPD = importdata('demandload.xlsx');%load demandload.xlsx form root
case30-v1(load curve).
% mpc.bus( :,3) = demloadPD.data.LFP( :,4);
% mpc.bus( :,4) = demloadPD.data.LFP( :,5);
%% power flow without DG
mpopt = mpoption('pf.alg','PQSUM');
results = radial_pf(mpc,mpopt);
%% calculat sence factor
sensfactor(:,1:2) = results.branch(:,1:2);% put index brach
sensfactor(:,3) = sum(abs(real(get_losses(results))));%get power loss in the
branch
sensfactor(:,4) = abs(results.branch(:,14));
sensfactor(:,5) = (sensfactor(:,3)./sensfactor(:,4));
sensfactor(:,6) = (real(get_losses(results)));
[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 (1:4,2); % chose 4th of end
%%general
loss = sum(real(get_losses(results)));
Pi_main = results.gen(1,2); %product main generator
%%cost function
f3 = Ceb *loss;
f1 = Pi_main*Ces ;
%% deflected voltage
Lanbada_w =sum((results.bus( :,8)>consmaxV | results.bus( :,8)<consminV));
Vy_w = find (results.bus( :,8)>consmaxV | results.bus( :,8)<consminV);%find
deflect voltage for one year
Vdeflecty_w(Vy_w) = find (results.bus( :,8)>consmaxV | results.bus(
:,8)<consminV);% find deflect voltage
%% loss in branch
Loss_b = (abs(real(get_losses(results))));
%%Energy not supply bus 25 -29 will be fault
penns = sqrt((mpc.bus(25,3)^2+ mpc.bus(25,4)^2))*1+sqrt((mpc.bus(26,3)^2+
mpc.bus(26,4)^2))*1 ...
+sqrt((mpc.bus(27,3)^2+ mpc.bus(27,4)^2))*1;
ENNS_w =Ruc*penns;%Energy not supply without DG* 3dolar
E_w = sqrt((results.gen(1,2)+results.gen(1,3))^2);%product GEN1 (MVA)
Breli_w = ENNS_w;
Bupda_w = sum(1.73*Cmar*( abs( results.branch(:,14))));
Bimpu_w = Lanbada_w*Bjmu;
Benvi_w = E_w*Cenvi*Bco2;
Bploss_w = loss*RGp;
%Bsub = Ceb(E_wo-E_w)-8760(Cdgi*Pi*CFi); Calculat After dession DG
Pi_dg = sum(results.gen(2:5,2));
%% IPS
BANN_w = Pi_main*Ces ; %annual yoeld of IPS in vestment
Cinv = sum(results.gen(2:5,2).*sensfactor(1:4,2))* Pi_main *Ci*Ciop;
Bips_w = BANN_w/Cinv ;
Bsub_w = (Ceb*E_w)-(f3*Pi_dg); %Calculat After dession DG
Bdisco =Bploss_w+Bimpu_w+Breli_w+Bupda_w+Benvi_w+Bsub_w;function mpc = w_dg_case37
%CASE30 Power flow data for 30 bus, 6 generator + 5 DG case.
% MATPOWER
%% MATPOWER Case Format : Version 2
mpc.version = '2';
%%----- Power Flow Data -----%%
%% system MVA base
mpc.baseMVA = 15;
%% bus data
% bus_i type Pd Qd Gs Bs area Vm Va
baseKV zone Vmax Vmin
mpc.bus = [
1 3 0 0 0 0 1 1 0
10 1 1.05 0.95;
2 1 0.630 0.315 0 0 1 1 0 10 1 1.05 0.95;
3 1 0 0 0 0 1 1 0 10 1 1.05 0.95;
4 2 0 0 0 0 1 1 0 10 1 1.05 0.95;
5 1 0.085 0.040 0 0 1 1 0 10 1 1.05 0.95;
6 2 0 0 0 0 1 1 0 10 1 1.05 0.95;
7 2 0 0 0 0 1 1 0 10 1 1.05 0.95;
8 1 0.085 0.040 0 0 1 1 0 10 1 1.05 0.95;
9 2 0.042 0.040 0 0 1 1 0 10 1 1.05 0.95;
10 1 0.140 0.070 0 0 1 1 0 10 1 1.05 0.95;
11 1 0.126 0.061 0 0 1 1 0 10 1 1.05 0.95;
12 2 0 0 0 0 1 1 0 10 1 1.05 0.95;
13 1 0.085 0.040 0 0 1 1 0 10 1 1.05 0.95;
14 2 0 0 0 0 1 1 0 10 1 1.05 0.95;
15 1 0.093 0.044 0 0 1 1 0 10 1 1.05 0.95;
16 1 0.085 0.040 0 0 1 1 0 10 1 1.05 0.95;
17 2 0 0 0 0 1 1 0 10 1 1.05 0.95;
18 2 0 0 0 0 1 1 0 10 1 1.05 0.95;
19 1 0.085 0.040 0 0 1 1 0 10 1 1.05 0.95;
20 2 0.042 0.021 0 0 1 1 0 10 1 1.05 0.95;
21 1 0.085 0.040 0 0 1 1 0 10 1 1.05 0.95;
22 1 0.038 0.018 0 0 1 1 0 10 1 1.05 0.95;
23 1 0.085 0.040 0 0 1 1 0 10 1 1.05 0.95;
24 2 0 0 0 0 1 1 0 10 1 1.05 0.95;
25 1 0.042 0.021 0 0 1 1 0 10 1 1.05 0.95;
26 2 0.166 0.080 0 0 1 1 0 10 1 1.05 0.95;
27 1 0.042 0.021 0 0 1 1 0 10 1 1.05 0.95;
28 1 0.042 0.021 0 0 1 1 0 10 1 1.05 0.95;
29 1 0.042 0.021 0 0 1 1 0 10 1 1.05 0.95;
30 1 0.126 0.063 0 0 1 1 0 10 1 1.05 0.95;
31 2 0 0 0 0 1 1 0 10 1 1.05 0.95;
32 1 0.085 0.040 0 0 1 1 0 10 1 1.05 0.95;
33 1 0.042 0.021 0 0 1 1 0 10 1 1.05 0.95;
34 2 0 0 0 0 1 1 0 10 1 1.05 0.95;
35 1 0.085 0.040 0 0 1 1 0 10 1 1.05 0.95;
36 1 0.042 0.021 0 0 1 1 0 10 1 1.05 0.95;
37 1 0.085 0.040 0 0 1 1 0 10 1 1.05 0.95;
];
%% generator data
% bus Pg Qg Qmax Qmin Vg mBase status Pmax
Pmin Pc1 Pc2 Qc1min Qc1max Qc2min Qc2max ramp_agc ramp_10
ramp_30 ramp_q apf
mpc.gen = [
1 11.52 0 1.50 -0.20 1 15 1 10 0 0 0 0
0 0 0 0 0 0 0 0
20 1 0 1.50 -0.20 1 15 1 10 0 0 0 0 0
0 0 0 0 0 0 0
21 1 0 1.50 -0.20 1 15 1 10 0 0 0 0 0
0 0 0 0 0 0 0
22 1 0 1.50 -0.20 1 15 1 10 0 0 0 0 0
0 0 0 0 0 0 0
37 1 0 1.50 -0.20 1 15 1 10 0 0 0 0 0
0 0 0 0 0 0 0
];
%% branch data
% fbus tbus r x b rateA rateB rateC ratio
angle status angmin angmax
mpc.branch = [
1 2 0.10 0.1 0 6.98 6.98 6.98 0
0 1 -360 360;
2 3 0.52 0.52 0 6.98 6.98 6.98 0 0
1 -360 360;
3 4 0.22 0.22 0 4.83 4.83 4.83 0
0 1 -360 360;
4 5 0.71 0.71 0 2.30 2.30 2.30 0 0
1 -360 360;
5 6 0.13 0.05 0 2.30 2.30 2.30 0 0
1 -360 360;
6 7 0.32 0.05 0 2.30 2.30 2.30 0 0
1 -360 360;
7 8 0.04 0.08 0 2.30 2.30 2.30 0 0
1 -360 360;
8 9 0.43 0.09 0 2.30 2.30 2.30 0 0
1 -360 360;
9 10 0.17 0.06 0 2.30 2.30 2.30 0 0
1 -360 360;
10 11 0.13 0.06 0 2.30 2.30 2.30 0 0
1 -360 360;
11 12 0.15 0.06 0 2.30 2.30 2.30 0 0
1 -360 360;
12 13 0.41 0.06 0 2.30 2.30 2.30 0 0
1 -360 360;
13 14 0.06 0.05 0 2.30 2.30 2.30 0 0
1 -360 360;
14 15 0.17 0.04 0 2.30 2.30 2.30 0 0
1 -360 360;
15 16 0.17 0.05 0 2.30 2.30 2.30 0 0
1 -360 360;
16 17 0.32 0.08 0 2.30 2.30 2.30 0 0
1 -360 360;
17 18 0.17 0.08 0 2.30 2.30 2.30 0 0
1 -360 360;
18 19 0.11 0.12 0 2.30 2.30 2.30 0 0
1 -360 360;
2 35 0.69 0.93 0 1.56 1.56 1.56 0 0 1
-360 360;
35 36 0.22 0.19 0 1.56 1.56 1.56 0 0 1
-360 360;
36 37 0.11 0.03 0 1.56 1.56 1.56 0 0 1
-360 360;
3 31 0.28 0.08 0 4.83 4.83 4.83 0 0
1 -360 360;
31 32 0.28 0.05 0 4.83 4.83 4.83 0 0
1 -360 360;
32 33 0.50 0.08 0 4.83 4.83 4.83 0 0
1 -360 360;
33 34 0.32 0.08 0 4.83 4.83 4.83 0 0
1 -360 360;
6 23 0.11 0.08 0 2.30 2.30 2.30 0 0
1 -360 360;
23 24 0.30 0.14 0 2.30 2.30 2.30 0 0 1
-360 360;
24 25 0.35 0.11 0 2.30 2.30 2.30 0 0 1
-360 360;
25 26 0.28 0.04 0 2.30 2.30 2.30 0 0 1
-360 360;
26 27 0.22 0.04 0 2.30 2.30 2.30 0 0 1
-360 360;
27 28 0.22 0.04 0 2.30 2.30 2.30 0 0 1
-360 360;
28 29 0.11 0.03 0 2.30 2.30 2.30 0 0 1
-360 360;
29 30 0.15 0.03 0 2.30 2.30 2.30 0 0 1
-360 360;
8 22 0.11 0.01 0 1.56 1.56 1.56 0 0 1
-360 360;
9 21 0.11 0.04 0 1.56 1.56 1.56 0 0 1
-360 360;
12 20 0.15 0.09 0 1.56 1.56 1.56 0 0 1
-360 360;
];
%%----- OPF Data -----%%
%% generator cost data
% 1 startup shutdown n x1 y1 ... xn
yn
% 2 startup shutdown n c(n-1) ... c0
% column 1
mpc.gencost = [
2 0 0 3 0.02 2 0;
2 0 0 3 0.02 2 0;
2 0 0 3 0.02 2 0;
2 0 0 3 0.02 2 0;
2 0 0 3 0.02 2 0;
];
%% convert branch impedances from Ohms to p.u.
[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
Vbase = mpc.bus(1, BASE_KV) * 1e3; %% in Volts
Sbase = mpc.baseMVA * 1e6; %% in VA
mpc.branch(:, [BR_R BR_X]) = mpc.branch(:, [BR_R BR_X]) / (Vbase^2 / Sbase);
function f = costfunction(x)
mpc = loadcase('w_dg_case37');
% objective function
mpopt = mpoption('pf.alg','PQSUM');
results = radial_pf(mpc,mpopt);
% assignin('base','results',results);
%constant
Ces = 490000; %cost of sold energy to consumer
490IRR=490000(IRR/MWh)
Ceb = 245000; %cost of buy energy from main grid
245IRR=245000(IRR/MWh)
Ruc = 326000; %cost of energy not supply
326IRR=326000(IRR/MWh)
Cenvi = 245000; % Co2 emission tax rate
Bco2 = 163000; % Co2 emission(Kg/MWh)
Bjmu = 10200000; %5000*0.15= 750*13600=10200000 sercice of
equpment
Cmar = 29400000 ; %490000*60 = 24900000
RGp = 264000; %price purchasing energy from the main grid
Ci = 735000; % investement cost
%Ciop = 0.1; %operation and maintenance
%%
mpc.gen(1,2) = x(1);
mpc.gen(2,2) = x(2);
mpc.gen(3,2) = x(3);
mpc.gen(4,2) = x(4);
mpc.gen(5,2) = x(5);
% % objective function
% assignin('base','results',results);
% % save results;
loss = sum(real(get_losses(results)));
assignin('base','loss',loss);
f3 = Ces *loss*8760; % cost of losses
assignin('base','f3',f3);
g2 =
sum(((results.gen(2,2).^2).*(results.gencost(2,4))+(results.gen(2,2)).*(results.gencost(2,5))+(results.gencost(2,6))));
g3 =
sum(((results.gen(3,2).^2).*(results.gencost(3,4))+(results.gen(3,2)).*(results.gencost(3,5))+(results.gencost(3,6))));
g4 =
sum(((results.gen(4,2).^2).*(results.gencost(4,4))+(results.gen(4,2)).*(results.gencost(4,5))+(results.gencost(4,6))));
g5 =
sum(((results.gen(5,2).^2).*(results.gencost(5,4))+(results.gen(5,2)).*(results.gencost(5,5))+(results.gencost(5,6))));
f2 = (g2+g3+g4+g5)*8760;% cost of operation and maintenance
% f2 = results.f;
assignin('base','f2',f2);
f1 = ((x(2)*Ceb)+(x(3)*Ceb)+(x(4)*Ceb)+(x(5)*Ceb))*8760; % DG owner income
assignin('base','f1',f1);
f = (-f1+f2)+f3;
end
