Skip to content

Commit

Permalink
Add MTES-KG TEVC 2023
Browse files Browse the repository at this point in the history
  • Loading branch information
@intLyc
intLyc committed Nov 2, 2023
1 parent 00b0d02 commit a8f35e2
Show file tree
Hide file tree
Showing 11 changed files with 940 additions and 0 deletions.
202 changes: 202 additions & 0 deletions MTO/Algorithms/Multi-task/Multi-population/MTES-KG/MTES_KG.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,202 @@
classdef MTES_KG < Algorithm
% <Multi-task/Many-task> <Single-objective> <None/Constrained>

%------------------------------- Reference --------------------------------
% @Article{Li2023MTES-KG,
% title = {Multitask Evolution Strategy with Knowledge-Guided External Sampling},
% author = {Li, Yanchi and Gong, Wenyin and Li, Shuijia},
% journal = {IEEE Transactions on Evolutionary Computation},
% year = {2023},
% }
%--------------------------------------------------------------------------

%------------------------------- Copyright --------------------------------
% Copyright (c) 2022 Yanchi Li. You are free to use the MTO-Platform for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "MTO-Platform" and cite
% or footnote "https://github.com/intLyc/MTO-Platform"
%--------------------------------------------------------------------------

properties (SetAccess = private)
tau0 = 2
alpha = 0.5
adjGap = 50
sigma0 = 0.3
end

methods
function Parameter = getParameter(Algo)
Parameter = {'tau0: External Sample Number', num2str(Algo.tau0), ...
'alpha: DoS/SaS Probability', num2str(Algo.alpha), ...
'adjGap: Gap of Adjust tau', num2str(Algo.adjGap), ...
'sigma0', num2str(Algo.sigma0)};
end

function Algo = setParameter(Algo, Parameter)
Algo.tau0 = str2double(Parameter{1});
Algo.alpha = str2double(Parameter{2});
Algo.adjGap = str2double(Parameter{3});
Algo.sigma0 = str2double(Parameter{4});
end

function run(Algo, Prob)
n = max(Prob.D); % dimension
lambda = Prob.N; % sample points number
mu = round(lambda / 2); % effective solutions number
weights = log(mu + 0.5) - log(1:mu);
weights = weights ./ sum(weights); % weights
mueff = 1 / sum(weights.^2); % variance effective selection mass
% expectation of the Euclidean norm of a N(0,I) distributed random vector
chiN = sqrt(n) * (1 - 1 / (4 * n) + 1 / (21 * n^2));
hth = (1.4 + 2 / (n + 1)) * chiN;
for t = 1:Prob.T
% Step size control parameters
cs{t} = (mueff + 2) / (n + mueff + 5);
damps{t} = 1 + cs{t} + 2 * max(sqrt((mueff - 1) / (n + 1)) - 1, 0);
% Covariance update parameters
cc{t} = (4 + mueff / n) / (4 + n + 2 * mueff / n);
c1{t} = 2 / ((n + 1.3)^2 + mueff);
cmu{t} = min(1 - c1{t}, 2 * (mueff - 2 + 1 / mueff) / ((n + 2)^2 + 2 * mueff / 2));
% Initialization
mDec{t} = mean(unifrnd(zeros(lambda, n), ones(lambda, n)));
ps{t} = zeros(n, 1);
pc{t} = zeros(n, 1);
B{t} = eye(n, n);
D{t} = ones(n, 1);
C{t} = B{t} * diag(D{t}.^2) * B{t}';
invsqrtC{t} = B{t} * diag(D{t}.^-1) * B{t}';
sigma{t} = Algo.sigma0;
eigenFE{t} = 0;
for i = 1:lambda
sample{t}(i) = Individual();
end
mStep{t} = 0; % mean sample step
numExS{t} = []; % external sapmle number memory
sucExS{t} = []; % external sapmle success number memory
tau(t) = Algo.tau0; % external sample number
record_tau{t} = tau(t);
end
rank = {};

while Algo.notTerminated(Prob)
%% Sample new solutions
oldsample = sample;
sample{t} = sample{t}(1:lambda);
for t = 1:Prob.T
for i = 1:lambda
sample{t}(i).Dec = mDec{t} + sigma{t} * (B{t} * (D{t} .* randn(n, 1)))';
end
mStep{t} = mean(sqrt(sum((sample{t}.Decs - mDec{t}).^2, 2)));
end

%% Sample external solutions
for t = 1:Prob.T
% Select auxiliary task
idx = 1:Prob.T; idx(t) = [];
k = idx(randi(end));

for i = lambda + 1:lambda + tau(t)
if Algo.Gen < 2
sample{t}(i).Dec = mDec{t} + sigma{t} * (B{t} * (D{t} .* randn(n, 1)))';
continue;
end
if rand() < Algo.alpha
% Optimal domain knowledge-guided external sampling (DoS)
sample_k = mDec{k} + sigma{k} * (B{k} * (D{k} .* randn(n, 1)))';
vec = (sample_k - mDec{t});
if norm(vec) < mStep{t}
sample{t}(i).Dec = sample_k;
else
uni_vec = vec ./ norm(vec);
sample{t}(i).Dec = mDec{t} + uni_vec * mStep{t};
end
else
% Function shape knowledge-guided external sampling (SaS)
idx = 1:mu; idx(randi(end)) = [];
vec = mean(oldsample{k}(rank{k}(idx)).Decs);
vec = (vec - mDec{k}) ./ sigma{k};
sample{t}(i).Dec = mDec{t} + sigma{t} * (B{t} * (D{t} .* (B{k}' * (D{k}.^-1 .* vec'))))';
end
end
end

%% Update algorithm parameters
for t = 1:Prob.T
rank{t} = Algo.EvaluationAndSort(sample{t}, Prob, t);

% Storage number and success of external samples
numExS{t}(Algo.Gen) = tau(t);
sucExS{t}(Algo.Gen) = length(find(rank{t}(1:mu) > lambda));

% Negative transfer mitigation
if mod(Algo.Gen, Algo.adjGap) == 0
numAll = sum(numExS{t}(Algo.Gen - Algo.adjGap + 1:Algo.Gen));
sucAll = sum(sucExS{t}(Algo.Gen - Algo.adjGap + 1:Algo.Gen));

if (numAll > 0 && sucAll / numAll > 0.5) || (numAll == 0)
tau(t) = min([Algo.tau0, tau(t) + 1]);
else
tau(t) = max(0, tau(t) - 1);
end
end

% Update CMA parameters
% Update mean decision variables
oldDec = mDec{t};
mDec{t} = weights * sample{t}(rank{t}(1:mu)).Decs;
% Update evolution paths
ps{t} = (1 - cs{t}) * ps{t} + sqrt(cs{t} * (2 - cs{t}) * mueff) * invsqrtC{t} * (mDec{t} - oldDec)' / sigma{t};
hsig = norm(ps{t}) / sqrt(1 - (1 - cs{t})^(2 * (ceil((Algo.FE - lambda * (t - 1)) / (lambda * Prob.T)) + 1))) < hth;
pc{t} = (1 - cc{t}) * pc{t} + hsig * sqrt(cc{t} * (2 - cc{t}) * mueff) * (mDec{t} - oldDec)' / sigma{t};
% Update covariance matrix
artmp = (sample{t}(rank{t}(1:mu)).Decs - repmat(oldDec, mu, 1))' / sigma{t};
delta = (1 - hsig) * cc{t} * (2 - cc{t});
C{t} = (1 - c1{t} - cmu{t}) * C{t} + c1{t} * (pc{t} * pc{t}' + delta * C{t}) + cmu{t} * artmp * diag(weights) * artmp';
% Update step size
sigma{t} = sigma{t} * exp(cs{t} / damps{t} * (norm(ps{t}) / chiN - 1));
% Check distribution correctness
if (Algo.FE - lambda * (t - 1)) - eigenFE{t} > (lambda * Prob.T) / (c1{t} + cmu{t}) / n / 10 % to achieve O(N^2)
eigenFE{t} = Algo.FE;
restart = false;
if ~(all(~isnan(C{t}), 'all') && all(~isinf(C{t}), 'all'))
restart = true;
else
C{t} = triu(C{t}) + triu(C{t}, 1)'; % enforce symmetry
[B{t}, D{t}] = eig(C{t}); % eigen decomposition, B==normalized eigenvectors
if min(diag(D{t})) < 0
restart = true;
else
D{t} = sqrt(diag(D{t})); % D contains standard deviations now
end
end
if restart
ps{t} = zeros(n, 1);
pc{t} = zeros(n, 1);
B{t} = eye(n, n);
D{t} = ones(n, 1);
C{t} = B{t} * diag(D{t}.^2) * B{t}';
sigma{t} = min(max(2 * sigma{t}, 0.01), 0.3);
end
invsqrtC{t} = B{t} * diag(D{t}.^-1) * B{t}';
end
record_tau{t} = [record_tau{t}; tau(t)];
end
end
end

function rank = EvaluationAndSort(Algo, sample, Prob, t)
%% Boundary constraint handling
boundCVs = zeros(length(sample), 1);
for i = 1:length(sample)
% Boundary constraint violation
tempDec = sample(i).Dec;
tempDec(tempDec < -0.05) = -0.05;
tempDec(tempDec > 1.05) = 1.05;
boundCVs(i) = sum((sample(i).Dec - tempDec).^2);
end
sample = Algo.Evaluation(sample, Prob, t);
boundCVs(boundCVs > 0) = boundCVs(boundCVs > 0) + max(sample.CVs);
[~, rank] = sortrows([sample.CVs + boundCVs, sample.Objs], [1, 2]);
end
end
end
130 changes: 130 additions & 0 deletions MTO/Problems/Real-world Applications/Optimal Power Flow/Case_IEEE_30.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,130 @@
function mpc = Case_IEEE_30
%CASE30 Power flow data for 30 bus, 6 generator case.
% Please see CASEFORMAT for details on the case file format.
%
% Based on data from ...
% Alsac, O. & Stott, B., "Optimal Load Flow with Steady State Security",
% IEEE Transactions on Power Apparatus and Systems, Vol. PAS 93, No. 3,
% 1974, pp. 745-751.
% ... with branch parameters rounded to nearest 0.01, shunt values divided
% by 100 and shunt on bus 10 moved to bus 5, load at bus 5 zeroed out.
% Generator locations, costs and limits and bus areas were taken from ...
% Ferrero, R.W., Shahidehpour, S.M., Ramesh, V.C., "Transaction analysis
% in deregulated power systems using game theory", IEEE Transactions on
% Power Systems, Vol. 12, No. 3, Aug 1997, pp. 1340-1347.
% Generator Q limits were derived from Alsac & Stott, using their Pmax
% capacities. V limits and line |S| limits taken from Alsac & Stott.

% MATPOWER

%% MATPOWER Case Format : Version 2
mpc.version = '2';

%%----- Power Flow Data -----%%
%% system MVA base
mpc.baseMVA = 100;

%% 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 135 1 1.1 0.95;
2 2 21.7 12.7 0 0 1 1 0 135 1 1.1 0.95;
3 1 2.4 1.2 0 0 1 1 0 135 1 1.05 0.95;
4 1 7.6 1.6 0 0 1 1 0 135 1 1.05 0.95;
5 2 94.2 19 0 0 1 1 0 135 1 1.1 0.95;
6 1 0 0 0 0 1 1 0 135 1 1.05 0.95;
7 1 22.8 10.9 0 0 1 1 0 135 1 1.05 0.95;
8 2 30 30 0 0 1 1 0 135 1 1.10 0.95;
9 1 0 0 0 0 1 1 0 135 1 1.05 0.95;
10 1 5.8 2 0 0 3 1 0 135 1 1.05 0.95;
11 2 0 0 0 0 1 1 0 135 1 1.10 0.95;
12 1 11.2 7.5 0 0 2 1 0 135 1 1.05 0.95;
13 2 0 0 0 0 2 1 0 135 1 1.10 0.95;
14 1 6.2 1.6 0 0 2 1 0 135 1 1.05 0.95;
15 1 8.2 2.5 0 0 2 1 0 135 1 1.05 0.95;
16 1 3.5 1.8 0 0 2 1 0 135 1 1.05 0.95;
17 1 9 5.8 0 0 2 1 0 135 1 1.05 0.95;
18 1 3.2 0.9 0 0 2 1 0 135 1 1.05 0.95;
19 1 9.5 3.4 0 0 2 1 0 135 1 1.05 0.95;
20 1 2.2 0.7 0 0 2 1 0 135 1 1.05 0.95;
21 1 17.5 11.2 0 0 3 1 0 135 1 1.05 0.95;
22 1 0 0 0 0 3 1 0 135 1 1.05 0.95;
23 1 3.2 1.6 0 0 2 1 0 135 1 1.05 0.95;
24 1 8.7 6.7 0 0 3 1 0 135 1 1.05 0.95;
25 1 0 0 0 0 3 1 0 135 1 1.05 0.95;
26 1 3.5 2.3 0 0 3 1 0 135 1 1.05 0.95;
27 1 0 0 0 0 3 1 0 135 1 1.05 0.95;
28 1 0 0 0 0 1 1 0 135 1 1.05 0.95;
29 1 2.4 0.9 0 0 3 1 0 135 1 1.05 0.95;
30 1 10.6 1.9 0 0 3 1 0 135 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 99.211 -3.99 150.0 -20 1.0 100 1 200 50 0 0 0 0 0 0 0 0 0 0;
2 80.00 50.0 60.0 -20 1.0 100 1 80 20 0 0 0 0 0 0 0 0 0 0;
5 50.00 37.0 62.5 -15 1.0 100 1 50 15 0 0 0 0 0 0 0 0 0 0;
8 20.00 37.3 48.7 -15 1.0 100 1 35 10 0 0 0 0 0 0 0 0 0 0;
11 20.00 16.2 40.0 -10 1.0 100 1 30 10 0 0 0 0 0 0 0 0 0 0;
13 20.00 10.6 44.7 -15 1.0 100 1 40 12 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.0192 0.0575 0.0528 130 130 130 0 0 1 -360 360;
1 3 0.0452 0.1652 0.0408 130 130 130 0 0 1 -360 360;
2 4 0.057 0.1737 0.0368 65 65 65 0 0 1 -360 360;
3 4 0.0132 0.0379 0.0084 130 130 130 0 0 1 -360 360;
2 5 0.0472 0.1983 0.0418 130 130 130 0 0 1 -360 360;
2 6 0.0581 0.1763 0.0374 65 65 65 0 0 1 -360 360;
4 6 0.0119 0.0414 0.009 90 90 90 0 0 1 -360 360;
5 7 0.046 0.116 0.0204 70 70 70 0 0 1 -360 360;
6 7 0.0267 0.082 0.017 130 130 130 0 0 1 -360 360;
6 8 0.012 0.042 0.009 32 32 32 0 0 1 -360 360;
6 9 0.00 0.208 0.00 65 65 65 0 0 1 -360 360;
6 10 0.00 0.556 0.00 32 32 32 0 0 1 -360 360;
9 11 0.00 0.208 0.00 65 65 65 0 0 1 -360 360;
9 10 0.00 0.11 0.00 65 65 65 0 0 1 -360 360;
4 12 0.00 0.256 0.00 65 65 65 0 0 1 -360 360;
12 13 0.00 0.14 0.00 65 65 65 0 0 1 -360 360;
12 14 0.1231 0.2559 0.00 32 32 32 0 0 1 -360 360;
12 15 0.0662 0.1304 0.00 32 32 32 0 0 1 -360 360;
12 16 0.0945 0.1987 0.00 32 32 32 0 0 1 -360 360;
14 15 0.221 0.1997 0.00 16 16 16 0 0 1 -360 360;
16 17 0.0524 0.1923 0.00 16 16 16 0 0 1 -360 360;
15 18 0.1073 0.2185 0.00 16 16 16 0 0 1 -360 360;
18 19 0.0639 0.1292 0.00 16 16 16 0 0 1 -360 360;
19 20 0.034 0.068 0.00 32 32 32 0 0 1 -360 360;
10 20 0.0936 0.209 0.00 32 32 32 0 0 1 -360 360;
10 17 0.0324 0.0845 0.00 32 32 32 0 0 1 -360 360;
10 21 0.0348 0.0749 0.00 32 32 32 0 0 1 -360 360;
10 22 0.0727 0.1499 0.00 32 32 32 0 0 1 -360 360;
21 22 0.0116 0.0236 0.00 32 32 32 0 0 1 -360 360;
15 23 0.10 0.202 0.00 16 16 16 0 0 1 -360 360;
22 24 0.115 0.179 0.00 16 16 16 0 0 1 -360 360;
23 24 0.132 0.27 0.00 16 16 16 0 0 1 -360 360;
24 25 0.1885 0.3292 0.00 16 16 16 0 0 1 -360 360;
25 26 0.2544 0.38 0.00 16 16 16 0 0 1 -360 360;
25 27 0.1093 0.2087 0.00 16 16 16 0 0 1 -360 360;
28 27 0.00 0.396 0.00 65 65 65 0 0 1 -360 360;
27 29 0.2198 0.4153 0.00 16 16 16 0 0 1 -360 360;
27 30 0.3202 0.6027 0.00 16 16 16 0 0 1 -360 360;
29 30 0.2399 0.4533 0.00 16 16 16 0 0 1 -360 360;
8 28 0.0636 0.20 0.0428 32 32 32 0 0 1 -360 360;
6 28 0.0169 0.0599 0.013 32 32 32 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
mpc.gencost = [
2 0 0 3 0.0 2.00 0.00375;
2 0 0 3 0.0 1.75 0.0175;
2 0 0 3 0.0 1.00 0.0625;
2 0 0 3 0.0 3.25 0.00834;
2 0 0 3 0.0 3.00 0.025;
2 0 0 3 0.0 3.00 0.025;
];
Loading

0 comments on commit a8f35e2

Please sign in to comment.