Simulated Annealing Driven by the Free-Energy: The SAFE Algorithm

1. Introduction

Optimization is one of the fundamental tools in natural sciences (from biological sciences to physical sciences, from economics to engineering sciences) [1,2,3]. It concerns how certain problems can be solved among a set of feasible alternatives based on certain goals or constraints. Mathematically an optimization problem can be formulated as an objective function that could represent a desired outcome to be maximized, such as profit or utility, or minimized, such as cost or error [4,5] and/or finding an optimal solution, that is, a feasible solution for which the optimal value of the objective function is achieved. The nature of the optimization problem mostly depends on the specific characteristics of the problem, such as the nature of the objective function, the constraints, and the required precision. Usually, the optimizer is either local or global, thus there are different kinds of optimization algorithm techniques (e.g., deterministic, stochastic, hybrid, etc.) [6,7,8]. Now, to find global minima researchers use mostly stochastic or genetic algorithms, ranging from Simulated annealing (SA) [9,10,11] to particle swarm optimization [12] and from Differential evolution [13] to stochastic ranking [14]). The question is if the dimension of any objective function becomes too large and symmetric, can we optimize the problem in a faster time? In other words, can we search for a global minimum for this kind of cost function in a faster manner? To solve this we propose a free-energy driven self-adaptive simulated annealing algorithm (SAFE algorithm). This algorithm can help us to minimize the dimension of the problem by choosing the random dimension given a minimum length of the multidimensional objective function and it calculates and compares the free energy cost of the solution and its previous solution. As free energy has two components, on one hand, the energy part minimizes the cost function or reduces the objective function over iterations. On the other hand, the entropy difference minimizes the dimension of the problem. Based on these it later uses the infamous Metropolis-Hastings algorithm [15]. Later, we apply our algorithm in two examples where in the first case, we choose (i) Convex unimodal function and (ii) Rastrigin function (Non-convex Multimodal function). Furthermore, we ask if the optimization speed is better than standard simulated annealing. Interestingly, we investigated these questions and found that free-energy-driven self-adaptive SA is much faster than standard SA. This distinctive approach underpins our algorithm’s enhanced performance, providing a novel perspective on integrating information theory into optimization. At last, we estimate the parameter set of a stochastic exponential equation.

This article is arranged in the following way. In Section 2 we formalize the free-energy-driven self-adaptive SA algorithms. After that, in Section 3 we choose two examples one from the convex function and one from the non-convex function, and search the global minimum values and optimization speed. Furthermore, we showed parameter estimation of a simple exponential decay stochastic equation as an application. At last, we discussed our results in Section 4.

2. Methods

In this section, we shall describe our Free energy-driven Stochastic-Annealing algorithm. So, first, we shall talk about the statistical physics of Free energy.

2.1. The Free Energy in Statistical Physics

In Statistical Physics, free energy is a crucial concept that quantifies the amount of work that a thermodynamic system can perform. It provides insights into the equilibrium properties and spontaneous processes within the system. The free energy combines the system’s internal energy and entropy, reflecting both the energy available and the degree of disorder.

The Helmholtz free energy, F, is particularly important in systems at constant temperature and volume. It is defined as:

$$F=U-TS,$$

(1)

where U is the internal energy, T is the absolute temperature, and S is the entropy of the system. The Helmholtz free energy is minimized at equilibrium, making it a useful function for determining the stability and spontaneity of processes.

Alternatively, the Gibbs free energy, G, is used for systems at constant temperature and pressure. It is given by:

$$G=H-TS,$$

(2)

where $H=U+PV$ is the enthalpy, with P being the pressure and V the volume of the system. The Gibbs free energy is minimized at equilibrium under these conditions, and it is particularly relevant for chemical reactions and phase transitions.

2.2. Simulated Annealing Driven by the Free-Energy: The SAFE Algorithm

Mathematical optimization problems are pretty diverse. It guides us to choose a set of best solutions out of feasible ones based on some inherent goals or constraints, and/or the target can be to make faster optimization time or convergence time. Optimization problems can be discrete or continuous. Depending on the problem we choose optimizers which can be local or global in nature. Simulated Annealing is one of the algorithms or an optimizer that is typically used to search the global minima of a function. It’s based on the probabilistic thermodynamic approach motivated by a metal cooling or annealing process. In this article, we develop a free-energy driven self-adaptive Simulated Annealing (SA) which can search the global minima of a function where the dimension of the function is completely random over time t. To do this, we explore the effectiveness of a modified SA algorithm incorporating entropy as an additional factor. Initially, we need to formulate the problem mathematically. Let us assume a multidimensional function defined as $f\left(\mathbf{x}\right)$, where $\mathbf{x}=(x_1,x_2,..,x_k)\in {\mathrm{R}}^k$ vector of size k. Interestingly, the dimension k can change over time or in each iteration n. The initialization of this problem can be written as

$$\begin{array}{cc}\hfill & \mathrm{Initialize}\mathrm{solution}{\mathbf{x}}_0\in {\mathbb{R}}^k\hfill \\ \hfill & \mathrm{Set}\mathrm{initial}\mathrm{temperature}{T}_0>0\hfill \\ \hfill & \mathrm{Set}\mathrm{cooling}\mathrm{rate}\gamma \in (0,1)\hfill \\ \hfill & \mathrm{Set}\mathrm{minimum}\mathrm{length}\mathrm{of}\mathrm{array}l\in (l_{min},l_{rand})\hfill \\ \hfill & \mathrm{Set}\mathrm{maximum}\mathrm{number}\mathrm{of}\mathrm{iterations}{K}_{\max }\hfill \end{array}$$

After initializing the problem, we iterate this multidimensional function using an SA algorithm where first it perturbs the candidate solution i.e., ${\mathbf{x}}^{\prime }=\mathbf{x}+\mathbf{ffi}\mathbf{x}$ using a uniform distribution that follows the range $(-\Delta ,\Delta )$. After that, it evaluates the objective function (or the energy $E\left({\mathbf{x}}^{\prime }\right))$ and calculates the corresponding numerical entropy of the function $S\left({\mathbf{x}}^{\prime }\right)$ followed by the free-energy function $F\left({\mathbf{x}}^{\prime }\right)$. Here $S\left({\mathbf{x}}^{\prime }\right)$ is defined as a joint entropy of multivariate distributions as

$$\begin{array}{c}S\left(x_1^{\prime },x_2^{\prime },\dots ,x_n^{\prime }\right)\hfill \\ \hfill =-\sum _{x_1^{\prime },x_2^{\prime },\dots ,x_n^{\prime }}P(x_1^{\prime },x_2^{\prime },\dots ,x_n^{\prime })logP(x_1^{\prime },x_2^{\prime },\dots ,x_n^{\prime })\end{array}$$

(3)

where $P(x_1^{\prime },x_2^{\prime },\dots ,x_n^{\prime })$ is the joint probability. We use intrinsic temperature $T_{int}$ which is a parameter that tells precisely that it will not change during the comparison stage of choosing the solution. Furthermore, we assess the cost of free energy as the optimization cost as ${(\Delta F)}^2$. Based on this Monte-Carlo criterion we use the free energy difference criterion to accept the new solution or the best solution. At last, we updated the temperature as the standard SA and used it in the Free energy relation as

$$T_{n+1}=\gamma T_n$$

(4)

The process terminates after $K_{max}$ iterations or when the temperature and cost function reach predefined minimum values. The pseudocode of the algorithm is presented below in Algorithm (1).

Algorithm 1: The SAFE algorithm

for $n=1$ to $K_{\max }$ do       Generate a new candidate solution ${\mathbf{x}}^{\prime }={\mathbf{x}}_n+\delta \mathbf{x}$       where $\delta \mathbf{x}\sim \mathcal{U}(-\Delta ,\Delta )$ (perturbation)      Choose a random array size of the function based on a minimum length of $l_{min}$      Evaluate the objective function $E\left({\mathbf{x}}^{\prime }\right)$      Evaluate the entropy $S\left({\mathbf{x}}^{\prime }\right)$      Compute energy difference $\Delta E=E\left({\mathbf{x}}^{\prime }\right)-E\left({\mathbf{x}}_n\right)$      Compute entropy difference $\Delta S=S\left({\mathbf{x}}^{\prime }\right)-S\left({\mathbf{x}}_n\right)$      Compute Free-energy difference: $\Delta F=\Delta E-\Delta S$ ( for $T_{int}$ fixed at 1)      if $\Delta F\le 0$ or $\tilde{u}\in \left[0,1\right]\le exp\left(-{\displaystyle \frac{\Delta F}{T_n}}\right)$ then      if accepted then         ${\mathbf{x}}_{n+1}={\mathbf{x}}^{\prime }$        else         ${\mathbf{x}}_{n+1}={\mathbf{x}}_n$        end if      end if      Reduce the temperature $T_{n+1}=\gamma T_n$   end for

Here $\tilde{u}$ is the random number generator taken from uniform distribution lies between $[0,1]$. The performance of this modified algorithm was compared with a standard simulated annealing algorithm.

3. Results

As we know optimization algorithms often help us to figure out the minima (i.e., local or global) of any function based on an optimization time or a convergence time. This optimization/convergence time is crucial because it helps us to find the nature of the full landscape of a cost function. So, in this section, we shall test some multidimensional functions and try to compare them with standard Simulated Annealing (SA) to see which can reach the global minimum in a shorter time.

3.1. Convex Unimodal Function

Convex unimodal functions are the simplest class of multidimensional functions with a single minima in a multidimensional space. These functions can be seen diversely. For example in physics from the prospect of potential energy landscapes of a harmonic trap to any parameter estimation problem of a mathematical model to calculate the cost function, we find these. Here we write the function as

$$f\left(\mathbf{x}\right)=\sum _{i=1}^N{(x_i-A)}^2$$

(5)

where A is a constant or the point where the minima is located. This objective function we use in our algorithm and found that the results are pretty close to A. We initialize our solution from a uniform random distribution of array size 100 with lower and upper bounds following as $[-100,100]$. As the array size changes in each iteration with a minimum length set at 2, we find the optimizer reduces its dimension from 100 to 2 after the end of the simulation. For target A set to 5 and for each number of simulations, we get the following results in Table 1. Here simulate the function 10 times to reach minima. After that, we plot energy, entropy, and cost function change over each iteration as in this Figure 1, we have compared the energy function over different iteration values for standard SA. This tells us that, the optimization time is faster for free-energy-driven self-adaptive SA than standard SA, which eventually takes a longer time to reach the minima. The self-adaptivity can be seen through the evolution of the temperature parameter (T). Interestingly, the convergence time of the standard SA increases with the initial dimension of the problem.

3.2. Rastrigin Function

To test the global minima of a more complicated function, we choose the Rastrigin function [16]. This function was first characterized by Rastrigin in 1974. The nature of the function is non-convex and it contains multiple local minima and the global minima lies at the value 0. The Rastrigin function is written as

$$f\left(\mathbf{x}\right)=BN+\sum _{i=1}^N\left[x_i^2-Bcos\left(2\pi x_i\right)\right]$$

(6)

Here B is a constant. In our simulation, we set the B value to $0.01$. Our results show that the global minima are situated around 0 can be seen in Table 2. Similar to the convex unimodal function, we record for the whole optimization time the energy function, entropy, and cost function over different iterations. Interestingly, we observe similar trends. As we know, this Rastrigin function can take a longer time to search global minima for any optimizing algorithm but our free-energy driven self-adaptive SA can search this global minimum in a shorter time which can be seen in Figure 2.

In Figure 1 and Figure 2 we try to understand how entropy optimization can play a key role in the SAFE algorithm. In the plot, the dimension or the length of the function can be seen in terms of the iterations. Interestingly, we can see that the dimension reduction happens for the SAFE algorithm over the iterations (which can be seen in Figure 1 and Figure 2) while in the case of the standard SA, it is not changing over the iterations. Even if we try with dynamic array size in standard SA over each iteration. So, the dimension of SA is constant.

3.3. Parameter Estimation in a Noisy Process

Parameter estimation [17] is a process where the parameters are estimated based on the observed data. This estimation is important from engineering to systems biology as it helps us to figure out a more realistic state of the system. As an application of our algorithm, we use this to estimate the parameters of an exponential fit as

$$x\left(t\right)=exp(-\alpha t)+\xi \left(t\right)$$

(7)

here noise $\xi$ is Gaussian in nature which has a mean 0 and standard deviation set to $0.1$. This equation includes rigorous applications from Physics (e.g., radioactive decay, the waiting time distribution, etc.) to biology (e.g., protein degradation, inhibition processes, etc.). Here $\alpha$ is a parameter that will be estimated by our free-energy-driven self-adaptive SA algorithm. We choose a 100 dimensional initial value of $\alpha$ and We shall fit the average true value of $\alpha$ and the fitted average value of $\alpha$. so the objective function can be defined as

$$f\left(\mathbf{x}\right)=\sum _{i=1}^N\sum _{j=1}^3{({\alpha }_{ij}-{\alpha }_j^{true})}^2$$

(8)

Now, using this convex unimodal objective function we have fit the exponential stochastic equation after averaging 10 simulations, where the mean value of ${\alpha }_j$ is $1.5616$, $2.1131$ and $1.0425$ the true value of the parameters are (i.e., ${\alpha }_{true}^j$) are $1.5$, 2 and 1. So, the fit is really close to the original one and we got a decay nature of the entropy, energy, and cost function over time. Also, this method can play a crucial if one has a high dimensional parameter set and initial conditions are completely unknown. So, this algorithm can be used in those scenarios using the difference between the average parameters of the true value and the average parameters which we shall estimate as an objective function or the cost function.

Figure 3. Fit of the exponential solution and corresponding energy, entropy, and cost function over ${10}^4$ iteration after averaging 10 simulations using SAFE optimization.

Figure 3. Fit of the exponential solution and corresponding energy, entropy, and cost function over ${10}^4$ iteration after averaging 10 simulations using SAFE optimization.

4. Discussion

In this article, we developed a free-energy-driven self-adaptive SA algorithm to search for the global minima of a multidimensional function, where the lengths of the variables are random over the iterations. As examples, we used both a convex function and a non-convex function. In both cases, our algorithm not only searched the global minima but also did so in a much shorter time than standard SA. Interestingly, the global minima of the two examples provided above do not depend on the dimension size. This suggests that entropy-driven self-adaptive SA can be employed in scenarios where dimensionality reduction can help to search global minima and can create faster convergence times. Finally, we demonstrated how this algorithm can be used for parameter estimation.

Free energy is a superposition of the system’s dynamics, quantified by the energy term, and information content of the system, quantified by the corresponding entropy. Here, we will use the concept of Helmholtz Free Energy to improve SA performance. Interestingly, one can visualize this aspect from Information thermodynamics [18] as a change in free energy due to the extraction of work after performing optimization. So, if one considers an initial solution of the optimization problem is prepared (and it follows a probability and at the end of the optimization task ( or by a Maxwell’s Demon [19]) it reaches to global minima and has a final probability distribution. In this way, one can write the fluctuations theorem [20] in the solution space and can define heat and work properties which can be translated further to the infamous Jarzynski’s inequality [21] if we consider the initial state prepared in an equilibrium condition and the terminal stage after optimization is in equilibrium as the free energy does not change over time.

There exist different kinds of modified versions of SA, where the temperature schemes are mainly varied [10] such that it can circumvent the local minima and eventually better search the global minima in a faster manner while bypassing the local minima. However, for higher dimensions, it appears to lack the dynamism required for escaping local minima effectively. On the other hand, our SAFE algorithm minimizes the dimension while minimizing the entropy to give liberty to the optimizer such that it can search the global minima more efficiently and conveniently concerning time. To check it, we compared furthermore our model with the modified SA where the temperature update or the scheme follows as: $T_{n+1}={\displaystyle \frac{T_n}{1+\beta T_n}}$ where $\beta >0$ [22]. Interestingly, we did not find a rapid change in the energy landscape for the modified SA which can be comparable to the SAFE algorithm. Instead, we observed that the SAFE algorithm consistently outperformed the modified SA in terms of finding lower energy solutions within a given computational budget. So, this SAFE algorithm can give us efficient robust solutions for a high-dimensional case.

Although SAFE shows promise as an optimization technique, some drawbacks can limit its effectiveness and applicability. First, the parameters depend on the symmetry of the problem. So, the global optimum is problem-specific. Free-energy-driven algorithms extend beyond SA and can be applied to genetic or evolutionary algorithms. The main goal is to find optimal solutions and/or global minima with high convergence rates while searching a broader space. One potential application of free-energy driven algorithms is the food foraging problem [23], where the agent or decision-maker optimizes its time to search for more food while foraging. In decision-making or game theory, when a decision-maker has to choose among multiple options within a constrained time or when decisions depend on the choices of other decision-makers, this algorithm can guide the decision-making process. Another example of a multidimensional objective function is a path optimization problem in network theory [24] . When dealing with a large set of paths in a network and trying to find the least path in a network among any set of points is challenging in nature due to the network topology. It will be interesting to observe how this algorithm can reduce convergence time in such network configurations.

The potential applications of this algorithm are vast, ranging from machine learning and complex system modeling to real-time optimization problems in various scientific and engineering domains. The scalability, robustness, and flexibility of our approach make it suitable for a wide range of challenging optimization problems, promising significant advancements in the field.