pyHarmonySearch

LICENSE

This software is licensed under the BSD 3-Clause License. Please refer to the separate LICENSE file for the exact text of the license. You are obligated to give attribution if you use this code.

ABOUT

pyHarmonySearch is a pure Python implementation of the harmony search (HS) global optimization algorithm. HS is a metaheuristic search algorithm that, similar to simulated annealing, tabu, and evolutionary searches, is based on real-world phenomena. Specifically, HS mimics a jazz band improvising together. Courtesy Wikipedia:

In the HS algorithm, each musician (= decision variable) plays (= generates) a note (= a value) for finding a best harmony (= global optimum) all together.

pyHarmonySearch supports both continuous and discrete variables and can take advantage of parallel processing using Python's multiprocessing module.

REQUIREMENTS

This code does not rely on any third-party software. It only requires Python 3.11 or higher.

INSTALL

pyHarmonySearch is available on PyPI at https://pypi.org/project/pyHarmonySearch/.

Install using pip:

pip install pyHarmonySearch

Install from source:

pip install .

USING THIS CODE

You use this code in three steps:

1. Implement your own objective function that inherits from ObjectiveFunctionInterface. It defines the problem: the parameters, their bounds or discrete values, which are fixed, and the fitness function. See objective_function_interface.py for every method and the examples for complete implementations.
2. Tune the input parameters, which your objective function returns. These are problem-specific, and the numbers used in the example implementations might not be appropriate for your problem.
   • hms (harmony memory size): how many harmonies are kept in memory.
   • hmcr (harmony memory considering rate): the probability of taking each value from memory rather than choosing it at random.
   • par (pitch adjusting rate): the probability of adjusting a value taken from memory.
   • mpap (maximum pitch adjustment proportion): the most a continuous value can move, as a proportion of its distance to the bound. Keep it between 0 and 1 so values stay within their bounds.
   • mpai (maximum pitch adjustment index): the most a discrete value can move, in indices. An adjusted discrete value always moves by at least one index, unless mpai is less than 1 or there is only one possible value.
   • max_imp (maximum improvisations): how many new harmonies to generate. This is the stopping criterion. Each search makes hms + max_imp fitness evaluations, counting the initial harmony memory.
   • maximize() chooses maximization or minimization, and an optional random seed makes runs reproducible.
3. Run HS with harmony_search(objective_function, num_processes, num_iterations, initial_harmonies=None). If given, initial_harmonies must contain exactly hms harmonies, each with one value per parameter; otherwise a ValueError is raised.

It will look something like this:

from multiprocessing import cpu_count

from pyharmonysearch import ObjectiveFunctionInterface, harmony_search


class ObjectiveFunction(ObjectiveFunctionInterface):
    ...  # implement every method of ObjectiveFunctionInterface


if __name__ == '__main__':
    obj_fun = ObjectiveFunction()
    num_processes = cpu_count()  # use number of logical CPUs
    num_iterations = num_processes * 5  # each process does 5 iterations
    results = harmony_search(obj_fun, num_processes, num_iterations)
    print(f'Elapsed time: {results.elapsed_time}\nBest harmony: {results.best_harmony}\nBest fitness: {results.best_fitness}')

More documentation is provided in harmony_search.py and objective_function_interface.py and in the examples.

Iterations and multiprocessing

Like many similar optimization algorithms, HS is stochastic, so you will get a slightly different result every time you run it. Because of this, harmony_search() runs num_iterations independent searches on num_processes processes using Python's multiprocessing module, and returns the best solution found across all of them.

Work is sent to a multiprocessing pool, so your objective function must be picklable (define it at the top level of a module), and the call to harmony_search() must be under if __name__ == '__main__':. This is required under the spawn and forkserver start methods: spawn is the default on macOS and Windows, and forkserver is the default on Linux from Python 3.14.

If your objective function already uses multiprocessing itself (e.g., in get_fitness()), use harmony_search_serial(objective_function, num_iterations, initial_harmonies=None) instead, which runs the iterations one after another in the current process. It isn't exported from the package, so import it with from pyharmonysearch.harmony_search import harmony_search_serial.

Results

HarmonySearchResults, a namedtuple, is returned. It has five fields:

• elapsed_time: a datetime.timedelta.
• best_harmony and best_fitness: the best harmony (a list with one value per parameter) found across all iterations, and its fitness.
• harmony_memories: one final harmony memory per iteration. Each is a list of (harmony, fitness) tuples.
• harmony_histories: one history per iteration. Each is a list of {'gen': generation, 'harmonies': harmony_memory} snapshots, taken every hms improvisations.

Constraints

HS has no built-in notion of constraints. Handle them in get_fitness() with a penalty that makes infeasible harmonies score worse than feasible ones, as the pressure vessel and knapsack examples do.

EXAMPLES

Each example is a standalone script that defines an ObjectiveFunction, runs HS, and prints the result. Each one's docstring explains the problem and its known optimum. Most use every CPU and take from seconds to a few minutes; 2-D_continuous_seed.py is the quickest. Run them from the examples/ directory after installing:

> ./2-D_continuous_seed.py
Elapsed time: 0:00:02.525815
Best harmony: [0.027035817109283933, -0.9950109785422445]
Best fitness: 3.9992441742581275

Feature demonstrations

These examples use simple problems to show one feature of the API at a time. The first four are variations on each other that maximize -(x^2 + (y+1)^2) + 4, whose maximum is 4 at (0, -1).

Problems from the literature

These examples solve well-known optimization problems with known optima, so you can judge how close HS gets.

TESTING

Install the development dependencies, then run the linter and the tests (which require 100% branch coverage) from the repository root:

pip install -e ".[dev]"
ruff check
pytest

REFERENCES

Neat Algorithms - Harmony Search (archived) provides a simple introduction on how HS works. Also, see the HS Wikipedia entry.

HS was first introduced by Geem et al. in 2001:

Z. W. Geem, J. H. Kim, and G. V. Loganathan, "A New Heuristic Optimization Algorithm: Harmony Search," Simulation, vol. 76, no. 2, pp. 60–68, Feb. 2001. doi:10.1177/003754970107600201

The problems in the literature examples come from:

• E. Sandgren, "Nonlinear Integer and Discrete Programming in Mechanical Design Optimization," Journal of Mechanical Design, vol. 112, no. 2, pp. 223–229, Jun. 1990. doi:10.1115/1.2912596 (pressure vessel)
• K. S. Lee and Z. W. Geem, "A new meta-heuristic algorithm for continuous engineering optimization: harmony search theory and practice," Computer Methods in Applied Mechanics and Engineering, vol. 194, no. 36–38, pp. 3902–3933, Sep. 2005. doi:10.1016/j.cma.2004.09.007 (pressure vessel solved with HS; the best known cost of about 6059.71 comes from later work)
• H. Mühlenbein, M. Schomisch, and J. Born, "The parallel genetic algorithm as function optimizer," Parallel Computing, vol. 17, no. 6–7, pp. 619–632, Sep. 1991. doi:10.1016/S0167-8191(05)80052-3 (generalized Rastrigin function)
• J. Burkardt, KNAPSACK_01: Data for the 01 Knapsack Problem (instance P07, from Kreher and Stinson)

Some modifications and improvements have been suggested to the original algorithm. None of these are implemented here. This package implements the original algorithm's structure (harmony memory, memory consideration, pitch adjustment, and replacing the worst harmony), with some choices of its own where the 2001 paper is silent or simpler: discrete pitch adjustment can move up to mpai values (the paper moves to a neighboring value) and moves inward from the first or last value; continuous pitch adjustment moves a value by up to mpap of its distance to a bound; and a new harmony identical to one already in memory is not added again.

• Z. W. Geem, "Improved Harmony Search from Ensemble of Music Players", in Knowledge-Based Intelligent Information and Engineering Systems, B. Gabrys, R. Howlett, and L. Jain, Eds. Springer Berlin / Heidelberg, 2006, pp. 86-93.
• M. Mahdavi, M. Fesanghary, and E. Damangir, "An improved harmony search algorithm for solving optimization problems", Applied Mathematics and Computation, vol. 188, no. 2, pp. 1567-1579, May 2007.
• M. G. H. Omran and M. Mahdavi, "Global-best harmony search", Applied Mathematics and Computation, vol. 198, no. 2, pp. 643-656, May 2008.
• Z. W. Geem and K.-B. Sim, "Parameter-setting-free harmony search algorithm," Applied Mathematics and Computation, vol. 217, no. 8, pp. 3881–3889, Dec. 2010.

DIFFERENCES

Another Python optimization package that implements HS:

• pygmo

pyHarmonySearch differs from other optimization packages in that it relies on no third-party software. It is simple and easy-to-use standalone Python code. The primary reason I developed this was because I got tired of having to compile and maintain my own third-party packages on a machine running several-year-old software for which I have only user-level access. SciPy is particularly nasty to build, and many scientific packages rely on it. I skirted this issue by implementing HS myself from scratch.