# Harmony Search in MATLAB

URL: https://yarpiz.com/harmony-search-ypea117
Author: Mostapha Kalami Heris
Published: 2015-09-01
Updated: 2026-10-30
Topics: Metaheuristics

> Implement Harmony Search in MATLAB with memory sampling, pitch adjustment, a shrinking bandwidth, and a bounded Sphere test problem.

Harmony Search (HS) is a global optimization algorithm inspired by the way musicians improvise harmonies. It was proposed by Zong Woo Geem, Joong Hoon Kim and G. V. Loganathan in 2001. Each solution is called a Harmony, and the algorithm keeps an archive of promising solutions called the Harmony Memory (HM). For more information, see the [Harmony Search article on Wikipedia](https://en.wikipedia.org/wiki/Harmony_search).

The analogy is simple. A musician who improvises can play a note from memory, change a remembered note a little, or play a new random note. HS builds each new solution in the same way, one decision variable at a time.

## How it works

At every iteration, new harmonies are generated. Each component $x_j$ of a new harmony is set by two rules.

**Memory consideration.** With probability HMCR (harmony memory consideration rate), the component is taken from a harmony chosen at random from the HM. With probability $1 - \mathrm{HMCR}$, it is a random value within the bounds of the variable:

$$
x_j =
\begin{cases}
x_j^{(k)},\ k \text{ chosen at random from the HM} & \text{with probability } \mathrm{HMCR} \\
\text{a uniform random value in } [x_j^{\min}, x_j^{\max}] & \text{with probability } 1 - \mathrm{HMCR}
\end{cases}
$$

**Pitch adjustment.** Then, with probability PAR (pitch adjustment rate), the component is shifted by a small random amount. The size of the shift is set by the Fret Width (FW), also called the bandwidth, like the small changes of pitch that guitar and violin players make:

$$
x_j \leftarrow x_j + \mathrm{FW}\cdot \delta
$$

Here $\delta$ is a random number. In this code, $\delta$ is drawn from the standard normal distribution.

Pitch adjustment does not change every component. Each component is adjusted only with probability PAR, which is 0.1 in the code. The new harmonies are then evaluated. If they are better than the worst harmonies in the HM, they take their place.

## About this implementation

This post shares an open-source implementation of Harmony Search in MATLAB. The algorithm is implemented in a structured way. If you know MATLAB, you should find the code easy to use in your research projects.

### What is in the download

The `Harmony Search` folder holds three files:

- `hs.m`: the HS script, with the problem definition, the parameters, the main loop and the results.
- `Sphere.m`: the sample cost function, $f(x) = \sum_{i=1}^{n} x_i^2$.
- `main.m`: a one line entry point that runs `hs.m`.

To run the code, open the `Harmony Search` folder in MATLAB and run `main.m`. The program prints the best cost of every iteration and then plots the best cost against the iteration number on a logarithmic scale.

The script creates random values with `unifrnd`, a function of the Statistics and Machine Learning Toolbox, so it probably needs that toolbox.

### Parameters in the code

| Parameter | Meaning | Value |
| --- | --- | ---: |
| `nVar` | Number of decision variables | 5 |
| `VarMin`, `VarMax` | Lower and upper bounds of the variables | $-10$, $10$ |
| `MaxIt` | Maximum number of iterations | 5000 |
| `HMS` | Harmony memory size | 25 |
| `nNew` | Number of new harmonies in each iteration | 20 |
| `HMCR` | Harmony memory consideration rate | 0.9 |
| `PAR` | Pitch adjustment rate | 0.1 |
| `FW` | Fret width (bandwidth) | $0.02\,(\text{VarMax} - \text{VarMin})$ |
| `FW_damp` | Damping ratio of the fret width | 0.995 |

### Details of this implementation

- The classic HS creates one new harmony per iteration. This code creates `nNew` harmonies, merges them with the HM, sorts them by cost and keeps the best `HMS`.
- The code applies pitch adjustment to every component with probability PAR, whether the component came from the HM or was chosen at random. In the classic HS, only components taken from the HM are pitch adjusted.
- The fret width is multiplied by `FW_damp` after every iteration, so the adjustments get smaller over time. A uniform shift, $\delta$ in $[-1, 1]$, is also in the code as a commented line.
- Positions outside the bounds are set back to the nearest bound.

To solve your own problem, write a cost function that takes a row vector and returns a number, and set `CostFunction`, `nVar`, `VarMin` and `VarMax` in the problem definition section of `hs.m`.

## Check the result

Start with the Sphere objective and verify that the best coordinates approach the origin while respecting the bounds. Keep the variable ranges in mind when setting fret width, since it controls the absolute size of pitch adjustments. This implementation creates multiple new harmonies per iteration, so compare evaluation budgets rather than iteration counts with a one-harmony version. Repeat runs before judging a parameter change.

## Related posts

- [YPEA](https://yarpiz.com/yarpiz-evolutionary-algorithms-ypea) is a MATLAB toolbox that includes HS and 14 other algorithms that share one problem definition.
- [Binary and Real-Coded Genetic Algorithms in MATLAB](https://yarpiz.com/genetic-algorithms-ypea101) is another population based method with the same code structure.

## How to cite

Mostapha Kalami Heris, Harmony Search in MATLAB (YPEA117), Yarpiz, 2015. https://yarpiz.com/harmony-search-ypea117