Harmony Search Algorithm : A music-inspired optimizer

When a musician composes or improvises the music he might choose from 3 strategies:

1. Reuse the line of notes from his memory (Selection/Exploitation).

2. Reuse the notes and tweak some of notes in the line (Exploration/Mutation).

3. Create the new line of notes (Exploration phase).

The function for each choice of operator can be defined in the parentheses.

Harmony Search (HS) was proposed by Geem (2001 ) and has been used for various optimization problems in engineering designs and in broder fields. HS is simple to code and suitable for both continuous and discrete problems. It should be noted that although HS is a population-based, the searching process is performed by a single agent but with long memory list i.e. memory size. It's opposed to other population based methods like GA, PSO, ACO, FA, etc.

In HS, we need few control parameters like

HMS = The maximum number of memory size (i.e. 20-30)

HMCR = Harmony memory considering rate (i.e. 0.9)

PAR = Pitch adjustment rate (i.e. 0.3)

bw = Bandwidth (i.e. 0.1(ub-lb))

bw_damp (to be changed for each iteration) = 0.99bw

The step to perform HS is explained by the pseudo code:

Given: f(x), lb[1..n], ub[1..n], HMS, HMCR, PAR, bw, maxIter

# --- Initialize Harmony Memory (HM) ---

for i = 1..HMS

X[i,1..n] = lb + (ub - lb) .* rand(1,n)

F[i] = f(X[i])

end

Sort HM so F ascending (best at top)

for iter = 1..maxIter

# --- Improvise one new harmony ---

for j = 1..n

if rand < HMCR

r = randi(HMS)

X_new[j] = X[r, j] # memory consideration

if rand < PAR

X_new[j] = X_new[j] + (2*rand-1)*bw[j] # pitch adjust

end

else

X_new[j] = lb[j] + rand*(ub[j]-lb[j]) # randomization

end

# optional: enforce integer vars here by rounding X_new[j]

# always clip bounds:

X_new[j] = min(max(X_new[j], lb[j]), ub[j])

end

F_new = f(X_new)

# --- Greedy replacement of worst ---

find iworst with F[iworst] = max(F)

if F_new < F[iworst]

X[iworst, 1..n] = X_new

F[iworst] = F_new

# optional: keep HM sorted (resort F and X together)

end

# optional: bw = bw * bwDamp

end

# best solution is at index of min(F)

---

Numerical Example

minimize x^2+y^2

s.t. x+y > 1,

-2 < x,y < 2

prob = struct();
prob.f  = @(x) x(1)^2 + x(2)^2;
prob.lb = [-2 -2];
prob.ub = [ 2  2];

% inequality g(x) <= 0; want x+y >= 1  -->  1 - (x+y) <= 0
prob.nonlcon = @(x) deal( 1 - (x(1)+x(2)), [] );
prob.isInt = [false,false];   % try [true,true] to make it integer

opts = struct('HMS',40,'HMCR',0.9,'PAR',0.3,'bw',0.5,...
              'bwDamp',0.997,'maxIter',3000,'penalty',1e6,'seed',2);

out = hs_opt(prob, opts);
disp(out.best.x), disp(out.best.f)

Answer:

Notes:

1. The main function setup the problem statement and constraints.

2. The HM function hs_opt(prob,opts) is called from the main program.

3. The program gives the optimal solution (x,y) and f_min.

function out = hs_opt(prob, opts)

% HARMONY SEARCH (continuous + integer, optional constraints)

% prob.f : @(x) -> scalar objective (lower is better)

% prob.lb, ub : [1 x n] vectors of bounds

% prob.isInt : logical [1 x n] (true for integer vars), optional

% prob.nonlcon : @(x) -> [c, ceq] with c<=0, ceq=0 (optional)

%

% opts.HMS : harmony memory size (default 30)

% opts.HMCR : memory considering rate (default 0.9)

% opts.PAR : pitch adjust rate (default 0.3)

% opts.bw : initial bandwidth (scalar or [1 x n]) (default 0.1*(ub-lb))

% opts.maxIter : iterations (default 2000)

% opts.bwDamp : bandwidth damping factor per iter (default 0.995)

% opts.penalty : base penalty coefficient (default 1e6)

% opts.seed : rng seed (optional)

.

.

.

Penalty function 1e6*max(0,c) is applied to each constraint.

% ---------- helpers ----------
clip = @(x) min(max(x, prob.lb), prob.ub);

    function fp = penalized_f(x)
        fx = prob.f(x);
        if isfield(prob,'nonlcon') && ~isempty(prob.nonlcon)
            [c,ceq] = prob.nonlcon(x);
            vio = 0;
            if ~isempty(c),   vio = vio + sum(max(0,c).^2); end
            if ~isempty(ceq), vio = vio + sum(ceq.^2);      end
            fp = fx + opts.penalty*vio;
        else
            fp = fx;
        end
    end

This hs_opt code can be applied to the real design problem which usually contains several complicated constraints.

---

Reference: Xin-She Yang (Ed.), Nature-Inspired Computation in Engineering (Studies in Computational Intelligence, 637), 2016

AE Kayabekir et al., Harmony Search Algorithm for Structural Engineering Problems, Springer, 2021