# Iterative optimization RCO: A "Ruler & Compass" deterministic method

Maurice Clerc

<figure><img src="image_1.png" /></figure>

## To cite this version:

Maurice Clerc. Iterative optimization RCO: A "Ruler & Compass" deterministic method. 2024. hal-04530935

HAL Id: hal-04530935

https://hal.science/hal-04530935

Preprint submitted on 4 Apr 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_1.png" /></figure>

HAL
open science

# Iterative optimization RCO: A "Ruler & Compass" deterministic method

Maurice Clerc

► To cite this version:

Maurice Clerc. Iterative optimization RCO: A "Ruler & Compass" deterministic method. 2024. [hal-04530935v2](https://hal.science/hal-04530935v2)

HAL Id: hal-04530935

https://hal.science/hal-04530935v2

Preprint submitted on 28 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Iterative optimization — RCO: A “Ruler & Compass”
deterministic method

Maurice Clerc*

2024-04-01

## 1 Introduction

Take a piece of paper and a pen, and draw any "reasonable" function (i.e., not entirely random), even one that varies greatly, as shown in the figure 1a. For many algorithms, finding the minimum would be challenging, but the Ruler & Compass Optimizer (RCO) described here requires very few iterations. Furthermore, it is deterministic. How is this possible?

No magic, just elementary geometry.

<footnote>*Maurice.Clerc@WriteMe.com</footnote>
# Iterative optimization — RCO: A “Ruler & Compass”
deterministic method

Maurice Clerc*

2024-09-27 DRAFT

## Abstract

We present here the basic version of an iterative optimization algorithm that is deterministic, needs just one parameter, and often capable of finding a good solution after very few evaluations of the fitness function. We explain its principles using a multimodal one-dimensional problem. For such problems, it could be applied using nothing more than a ruler and a compass. We provide classical examples illustrating the algorithm’s properties, as well examples where it fails. As this version does not address possible stagnation, it is really only useful for low-dimensional problems (typically at most ten), where each evaluation of a position in the search space is very costly in terms of computational resources.

## 1 Introduction

Many optimization problems, especially in engineering, have the following characteristics: they can only be solved by stochastic iterative methods; the number of variables is not very high; each iteration requires very costly computer resources. In addition, the user often has to adjust several parameters and many tests are required. This is why we propose a deterministic algorithm that can provide a good solution in very few iterations even on problems that are usually seen as difficult, and for which just one parameter is needed (and moreover often easy to define).

Several definitions of the theoretical *a priori* difficulty of an optimization problem (minimization here) have been proposed ([3]). For example, one can consider the number of local minima and the relative size of the subspace of acceptable solutions. If this number is high and the size small, then any iterative algorithm functioning as a black box would require numerous samples of positions in the search space and, therefore, many evaluations of positions.

Consider the function with the landscape shown in Figure 1 where the minimum is 0.999507 at position 2.412616. Although it is one-dimensional, its theoretical difficulty is high because it has many local minima.

For the algorithms tested below, we say that a solution is acceptable if its value is less than 0.9997.

For stochastic algorithms, we estimate the probability of success by running them a thousand times.

Table 1 then shows that, for these, many evaluations are needed to have a good chance of finding a satisfactory solution and, therefore, the actual difficulty (large number of evaluations) is in good agreement with the theoretical difficulty.

However, it is possible to define a deterministic (in the non-stochastic sense) algorithm, for which this agreement is challenged. On this problem, the Ruler & Compass Optimizer (RCO) is extremely efficient (see figure 2). More generally, it is also effective on other supposedly difficult problems. And conversely, it can be highly ineffective on problems that are supposed to be easy.

<footnote>*Maurice.Clerc@WriteMe.com</footnote>
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>Figure 1: RCO on a test function.</figcaption></figure>

## 2 Principle

The iterative construction of positions to sample in the search space can be detailed on a very simple example, defined by

$$f(x) = (x - 3)^2 \text{ for } x \in [0, 10]$$

The only parameter is an estimated lower bound. Here -5 to better see what happens, although -0.1 would be sufficient.

Tables 1 to 5 illustrate the process. Of course, in higher dimension D lines become (hyper-)planes but the method is the same. At each time step we consider D planes defined by D + 1 points and their intersection point with the lower bound plane. There are two cases:
<figure><img src="image_3.png" /><figcaption>Figure 1: Test landscape sinCos15. Many local minima, very small set of acceptable solutions.</figcaption></figure>

Thus, the mere existence of such an algorithm largely calls into question the relevance of classical difficulty estimates, which will have to be revised, at least by specifying their field of application in the space of optimizers.

Informally, the No Free Lunch Theorem ([6]) allows us to say that there is no such thing as an efficient optimizer for all problems. If “efficient” is to be interpreted as “finding a solution without too much difficulty”, then what RCO suggests is a dual proposition: there is no a priori measure of difficulty valid for all optimizers.

Analyzing this conjecture is not the purpose of this study. In what follows, we'll look at the detailed operation of RCO and show on other examples how its behavior is atypical. This will allow us to sketch a typology of problems for which it is interesting to use it.

## 2 RCO principle

On our Test landscape RCO finds an acceptable solution for sure after 20 evaluations (see the figure 2).

However to detail the iterative construction it is better to consider a simpler example, defined by

$$f(x) = (x - 3)^2 \quad \text{for } x \in [0, 10] \qquad (1)$$

Its minimum is obviously 0 on $x = 3$.

For this basic version a user-defined parameter is needed: a fixed lower bound. Note that in practice we often know one. Moreover RCO is in fact able to use an adaptive one, automatically estimated during the process (see Appendix 8.5), but the figures are then more difficult to understand.

Here we set this lower bound to -5, to more clearly see what happens, although -0.1 would be sufficient.

Tables 2 to 6 illustrate the process. Of course, in higher dimension $D$ the lines become (hyper-)planes but the method is the same. At each time step we consider $D$ planes defined by $D+1$ points and their intersection point<fnref n="1" /> with the lower bound plane. There are two cases:

<footnote n="1">¹It may happen that there is more than an unique intersection point. If so we consider we are in the case 1.</footnote>
* If its position lies outside the search space (possibly at infinity), then the new position is determined by a weighted combination of the $2^D$ previous positions.

* If its position is within the search space, then it is kept as the new position.

<figure><img src="image_6.png" /><figcaption>Table 2: Principle. Step 2</figcaption></figure>

<figure><img src="image_7.png" /></figure>
<table><caption>Table 1: On the problem in figure 1, metaheuristics, as expected, sometimes struggle to find an acceptable solution with a high probability, estimated over 1000 executions. Since RCO is deterministic, only one run is required.</caption><thead><tr><th>Algorithm</th><th>Stochastic</th><th>Evaluations/run</th><th>Success rate</th></tr></thead><tbody><tr><td>CMA-ES</td><td>*</td><td>10 000</td><td>0.12</td></tr><tr><td>DE</td><td>*</td><td>10 000</td><td>0.73</td></tr><tr><td rowspan="2">ACO</td><td>*</td><td>410</td><td>0.80</td></tr><tr><td></td><td>4010</td><td>0.94</td></tr><tr><td>Jaya</td><td>*</td><td>1000</td><td>0.985</td></tr><tr><td rowspan="2">SPSO</td><td>*</td><td>180</td><td>0</td></tr><tr><td></td><td>190</td><td>1</td></tr><tr><td rowspan="2">GA-MPC</td><td>*</td><td>80</td><td>0</td></tr><tr><td></td><td>90</td><td>1</td></tr><tr><td>RCO</td><td></td><td>20</td><td>1</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>Figure 2: RCO on the test landscape sinCos15. Path followed (the first position is 0). An acceptable solution is found after 20 evaluations.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Table 4: Principle. Step 9, zoom</figcaption></figure>

<figure><img src="image_9.png" /></figure>
1. If its position lies outside the search space (possibly at infinity), then the new position is determined by a barycenter of the $2^D$ previous positions. For the calculation of the barycenter, the weight of each position is inversely proportional to the value of the function at that position. See the detail in Appendix 8.2.

2. If its position is within the search space, then it is kept as the new position.

Then, in the list of positions, the first one is removed and the new one is added.

<table><caption>Table 2: Principle. Step 1</caption><tr><td>An horizontal lower bound is defined. The “corners” are evaluated. The line connecting them intersects the lower bound beyond the search space, this intersection point is rejected and will be replaced with another one (step 2).</td></tr></table>

<figure><img src="image_5.png" /></figure>
<table><caption>Table 5: Principle. Step 10, and zoom</caption>
   <tr>
    <td>
     After Step 9, the line connecting the last two points is nearly horizontal and intersects the lower boundary outside of the search space. Consequently, we generate a new position using the same method as in Step 2. This occurrence becomes increasingly frequent as the positions approach the solution. Hence, the new position often approximates the mean of the last two positions, which is a favorable behavior for better approaching the solution.
    </td>
   </tr>
  </table>

<figure><img src="image_10.png" /></figure>

# 3 Examples

The definitions of the problems are given in the source codes (8.4.1). Some classical problems have been shifted in order to be not too easy.

## 3.1 Six Hump Camel Back

The dimension of the problem is two, so the landscape has four corners. The search space is $[-2, 2] \times [-1, 1]$.
<figure><img src="image_6.png" /><figcaption>Table 4: Principle. Step 3</figcaption></figure>

<figure><img src="image_7.png" /></figure>
The minimum is -1.031628453489877. A sure lower bound is -1.1.

<figure><img src="image_11.png" /><figcaption>Table 6: Six Hump Camel Back</figcaption></figure>
<table><caption>Table 5: Principle. Step 9, zoom</caption><tr><td>...</td><td>and again.</td></tr></table>

<figure><img src="image_8.png" /></figure>
<figure><img src="image_12.png" /><figcaption>Figure 2: Six Hump Camel Back. Lower bound -1.1. After the construction of 4+10 points the result is -0.9575.</figcaption></figure>

## 3.2 Shifted Rastrigin

The minimum is zero so a sure lower bound is -0.1.

<table><caption>Table 7: Shifted Rastrigin</caption><thead><tr><td>Dimension</td><td>Constructed points (evaluations)</td><td>Best fitness</td></tr></thead><tbody><tr><td>1</td><td>52</td><td>0.000557</td></tr><tr><td>2</td><td>1004</td><td>0.001511</td></tr><tr><td>3</td><td>1508</td><td>0.002548</td></tr><tr><td>4</td><td>50016</td><td>0.001735</td></tr><tr><td>5</td><td>250032</td><td>0.004518</td></tr></tbody></table>
<figure><img src="image_9.png" /><figcaption>Table 6: Principle. Step 10, and zoom</figcaption></figure>

# 3 Examples

The definitions of the problems are given in Appendix 8.7. Some of them have been shifted in order to be not too easy.

## 3.1 Six Hump Camel Back

The dimension of the problem is two, so the landscape has four corners. The search space is [-2, 2] x [-1, 1]. The minimum is -1.031628453489877. A sure lower bound is -1.1. Table 7 shows the results for increasing
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>Figure 3: Shifted Rastrigin</figcaption></figure>
number of evaluations and Figure 3 the construction for 14 evaluations.

<table><caption>Table 7: Six Hump Camel Back</caption><thead><tr><td>Constructed points (evaluations)</td><td>Best fitness</td></tr></thead><tbody><tr><td>14</td><td>-0,957541</td></tr><tr><td>24</td><td>-1,030227</td></tr><tr><td>34</td><td>-1,030227</td></tr><tr><td>44</td><td>-1,031227</td></tr><tr><td>54</td><td>-1,031227</td></tr><tr><td>64</td><td>-1,031473</td></tr><tr><td>74</td><td>-1,031473</td></tr><tr><td>84</td><td>-1,031473</td></tr><tr><td>94</td><td>-1,031473</td></tr><tr><td>104</td><td>-1,031473</td></tr></tbody></table>

<figure><img src="image_10.png" /></figure>
## 3.3 Rosenbrock

The minimum is zero so a sure lower bound is -0.1.

<table><caption>Table 8: Rosenbrock</caption><thead><tr><td>Dimension</td><td>Evaluations</td><td>Best fitness</td></tr></thead><tbody><tr><td>2</td><td>252</td><td>0.002788</td></tr><tr><td>3</td><td>1508</td><td>0.001003</td></tr><tr><td>4</td><td>10016</td><td>0.003756</td></tr><tr><td>5</td><td>15032</td><td>0.005719</td></tr></tbody></table>

## 3.4 Pressure Vessel

There are four values to find, the two first ones are discrete. The minimum value is 6059.714335048436, as proved in Yang et al. 2013.

The lower bound used here is 6000. Due to the discrete variables, the algorithm encounters difficulties in converging.

<figure><img src="image_15.png" /><figcaption>Table 9: Pressure Vessel</figcaption></figure>

<table><thead><tr><td>Evaluations</td><td>Best fitness</td></tr></thead><tbody><tr><td>1052</td><td>11056.81</td></tr><tr><td>2049</td><td>6646.71</td></tr><tr><td>3074</td><td>6615.977</td></tr><tr><td>4074</td><td>6615.977</td></tr><tr><td>5074</td><td>6615.977</td></tr><tr><td>6076</td><td>6592.011</td></tr></tbody></table>

## 3.5 Gear Train

There are four values to find, all discrete. The minimum is $2.700857 \times 10^{-12}$.

Here the lower bound is simply set to 0. The algorithm quickly proposes a solution that is far from optimal, and then stagnates without further improvement.
<figure><img src="image_11.png" /><figcaption>Figure 3: Six Hump Camel Back. Lower bound -1.1. After the construction of 4+10 points the result is -0.9575.</figcaption></figure>

## 3.2 Shifted Rastrigin

The minimum is zero, so a sure lower bound is -0.1. Table 8 shows the results for different dimensions and different numbers of evaluations and Figure 4 the constructions for 1D and 2D landscapes after 12 and 14 evaluations respectively.

<table><caption>Table 8: Shifted Rastrigin</caption><thead><tr><td>Dimension</td><td>Constructed points (evaluations)</td><td>Best fitness</td></tr></thead><tbody><tr><td>1</td><td>52</td><td>0.000557</td></tr><tr><td>2</td><td>1004</td><td>0.001511</td></tr><tr><td>3</td><td>1508</td><td>0.002548</td></tr><tr><td>4</td><td>50016</td><td>0.001735</td></tr><tr><td>5</td><td>250032</td><td>0.004518</td></tr></tbody></table>
<figure><img src="image_16.png" /><figcaption>Table 10: Gear Train. Discretization is applied only at the very end to all values, so more search effort does not always imply a better result.</figcaption></figure>

<table><thead><tr><td>Evaluations</td><td>Best fitness</td></tr></thead><tbody><tr><td>116</td><td>3.10E-03</td></tr><tr><td>130</td><td>1.20E-03</td></tr><tr><td>173</td><td>6.65E-09</td></tr><tr><td>223</td><td>6.65E-09</td></tr><tr><td>1753</td><td>6.65E-09</td></tr></tbody></table>

## 4 Comparison

Let's compare RCO with a slightly improved version of Standard PSO (SPSO-Dicho on my technical website (Clerc 2024b)). To slightly favor PSO, I've selected the best run out of five. The results are presented in 11. Also, refer to the figure 4 for Shifted Griewank, where the two methods appear to be equivalent, particularly for dimensions 1 to 5 and the number of evaluations respectively set to 100, 200, 5000, 10000, and 20000.

In lower dimensions (less than 10), RCO performs better and sometimes significantly so. However, for higher dimensions, RCO stagnates while PSO continues to find better solutions, as evidenced by Shifted Rastrigin and Rosenbrock.
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>Figure 4: Shifted Rastrigin</figcaption></figure>
<table><caption>Table 11: RCO vs PSO</caption><thead><tr><th>Problem</th><th>Dimension</th><th>Evaluations</th><th>RCO</th><th>PSO</th></tr></thead><tbody><tr><td rowspan="2">Six Hump Camel Back</td><td rowspan="2">4</td><td>54</td><td><strong>-1.031277</strong></td><td>-1.011030</td></tr><tr><td>504</td><td><strong>-1.031473</strong></td><td>-1.031628</td></tr><tr><td rowspan="4">Shifted Rastrigin</td><td>1</td><td>52</td><td><strong>0.000557</strong></td><td>0.510724</td></tr><tr><td>2</td><td>1004</td><td><strong>0.001511</strong></td><td>0.4397499</td></tr><tr><td>5</td><td>250032</td><td><strong>0.004578</strong></td><td>0.994959</td></tr><tr><td>10</td><td>251024</td><td>35.246</td><td><strong>1.9899</strong></td></tr><tr><td rowspan="3">Rosenbrock</td><td>2</td><td>252</td><td><strong>0.002788</strong></td><td>266.94</td></tr><tr><td>5</td><td>15032</td><td><strong>0.005719</strong></td><td>0.0864</td></tr><tr><td>10</td><td>101024</td><td>8.996</td><td><strong>5.19e-08</strong></td></tr><tr><td rowspan="3">Pressure Vessel</td><td rowspan="3">4</td><td>3074</td><td><strong>6615.977</strong></td><td>6890.347</td></tr><tr><td>5074</td><td><strong>6615.977</strong></td><td>6821.931</td></tr><tr><td>20074</td><td><strong>6583.75</strong></td><td>6820.410</td></tr><tr><td rowspan="3">Gear Train</td><td rowspan="3">4</td><td>173</td><td>6.65e-09</td><td><strong>2.35e-09</strong></td></tr><tr><td>1154</td><td>6.65e-09</td><td><strong>2.35e-09</strong></td></tr><tr><td>199752</td><td>1.18e-09</td><td><strong>1.26e-09</strong></td></tr></tbody></table>

<figure><img src="image_17.png" /><figcaption>Figure 4: Shifted Griewank. RCO and PSO are equivalent.</figcaption></figure>
## 3.3 Rosenbrock

The minimum is zero, so a sure lower bound is -0.1. Table 9 shows the results for dimensions two to five with an increasing number of evaluations.

<table><caption>Table 9: Rosenbrock</caption><thead><tr><td>Dimension</td><td>Evaluations</td><td>Best fitness</td></tr></thead><tbody><tr><td>2</td><td>252</td><td>0.002788</td></tr><tr><td>3</td><td>1508</td><td>0.001003</td></tr><tr><td>4</td><td>10016</td><td>0.003756</td></tr><tr><td>5</td><td>15032</td><td>0.005719</td></tr></tbody></table>

## 3.4 Pressure Vessel

There are four values to find, the two first ones are discrete. The minimum value is 6059.714335048436, as proved in [7].

The lower bound used here is 6000. Table 10 shows the results for an increasing number of evaluations. Due to the discrete variables, the algorithm encounters difficulties in converging.

<figure><img src="image_14.png" /><figcaption>Table 10: Pressure Vessel</figcaption></figure>

<table><thead><tr><td>Evaluations</td><td>Best fitness</td></tr></thead><tbody><tr><td>1052</td><td>11056.81</td></tr><tr><td>2049</td><td>6646.71</td></tr><tr><td>3074</td><td>6615.977</td></tr><tr><td>4074</td><td>6615.977</td></tr><tr><td>5074</td><td>6615.977</td></tr><tr><td>6076</td><td>6592.011</td></tr></tbody></table>

## 3.5 Gear Train

There are four values to find, all discrete. The minimum is $2.700857 \times 10^{-12}$.
# 5 Complexities

To evaluate an algorithm a common practice is to perform a theoretical analysis of its space complexity and of its time complexity. Note that, however, and as explained in Clerc 2024a, this is not always pertinent and a safer approach is to estimate the real memory space that is used and the real computing time.

Here, at each time step, there are two possible cases:

1. solving a linear system of $D$ equations to find a new position with $D$ coordinates. The needed space is $D^2 + D$ and the time complexity $O(D^3)$.

2. computing a linear combination of $2^D$ positions with $D$ coordinates. The needed space is $D \times 2^D$ and the number of multiplications is $2^D$.

To assess an algorithm, it's common practice to conduct a theoretical analysis of both its space complexity and time complexity. However, as explained in Clerc 2024a, this approach isn't always appropriate, and a more reliable method is to estimate the actual memory usage and computational time.

In this context, there are two possible scenarios at each time step:

1. Solving a linear system of $D$ equations to determine a new position with $D$ coordinates. The required space is $D^2 + D$, with a time complexity of $O(D^3)$.

2. Computing a linear combination of $2^D$ positions, each with $D$ coordinates. The required space is $O(D \times 2^D)$, and the number of multiplications is $O(2^D)$.

However such a classical reasoning does not correspond to the behavior of the algorithm in practice. The point is that the second situation usually does not occur very often. Moreover it depends on the landscape of the problem and also on the “budget” (the number of evaluations).

Let's consider, for example, the computing time per evaluation for two problems: Planes and Shifted Rastrigin. We have already seen the Shifted Rastrigin and Planes is defined by

$$f(x) = \sum_{d=1}^{D} (x(d) - 3) \qquad (1)$$

The figure 5 depicts its landscape in two dimensions.
Here the lower bound is simply set to 0. A we can see on Table 11 the algorithm quickly proposes a solution that is far from optimal, and then stagnates without further improvement.

<figure><img src="image_15.png" /><figcaption>Table 11: Gear Train. Discretization is applied only at the very end to all values, so more search effort does not always imply a better result.</figcaption></figure>

<table><thead><tr><td>Evaluations</td><td>Best fitness</td></tr></thead><tbody><tr><td>116</td><td>3.10E-03</td></tr><tr><td>130</td><td>1.20E-03</td></tr><tr><td>173</td><td>6.65E-09</td></tr><tr><td>223</td><td>6.65E-09</td></tr><tr><td>1753</td><td>6.65E-09</td></tr></tbody></table>

# 4 Comparison

Let's compare RCO with a slightly improved version of Standard PSO (SPSO-Dicho on my technical website ([4]). To favor PSO, I've selected the best run out of five. The results are presented in 12. Also, refer to the figure 5 for Shifted Griewank, where the two methods appear to be equivalent, particularly for dimensions 1 to 5 and the number of evaluations respectively set to 100, 200, 5000, 10000, and 20000.

In lower dimensions (less than 10), RCO performs better and sometimes significantly so. However, for higher dimensions, RCO stagnates while PSO continues to find better solutions, as evidenced by Shifted Rastrigin and Rosenbrock.
<figure><img src="image_18.png" /><figcaption>Figure 5: 2D Planes problem landscape.</figcaption></figure>

As we can see on table 12 the real computing time indeed grows exponentially, but more like $1.35^D$ than like $2^D$.
<table><caption>Table 12: RCO vs PSO</caption><thead><tr><th>Problem</th><th>Dimension</th><th>Evaluations</th><th>RCO</th><th>PSO</th></tr></thead><tbody><tr><td>Six Hump Camel Back</td><td>2</td><td>54</td><td><strong>-1.031227</strong></td><td>-1.011030</td></tr><tr><td rowspan="4">Shifted Rastrigin</td><td>1</td><td>52</td><td><strong>0.000557</strong></td><td>0.510724</td></tr><tr><td>2</td><td>1004</td><td><strong>0.001511</strong></td><td>0.4397499</td></tr><tr><td>5</td><td>250032</td><td><strong>0.004578</strong></td><td>0.994959</td></tr><tr><td>10</td><td>251024</td><td>35.246</td><td><strong>1.9899</strong></td></tr><tr><td rowspan="3">Rosenbrock</td><td>2</td><td>252</td><td><strong>0.002788</strong></td><td>266.94</td></tr><tr><td>5</td><td>15032</td><td><strong>0.005719</strong></td><td>0.0864</td></tr><tr><td>10</td><td>101024</td><td>8.996</td><td><strong>5.19e-08</strong></td></tr><tr><td rowspan="3">Pressure Vessel</td><td rowspan="3">4</td><td>3074</td><td><strong>6615.977</strong></td><td>6890.347</td></tr><tr><td>5074</td><td><strong>6615.977</strong></td><td>6821.931</td></tr><tr><td>20074</td><td><strong>6583.75</strong></td><td>6820.410</td></tr><tr><td rowspan="3">Gear Train</td><td rowspan="3">4</td><td>173</td><td>6.65e-09</td><td><strong>2.35e-09</strong></td></tr><tr><td>1154</td><td>6.65e-09</td><td><strong>2.35e-09</strong></td></tr><tr><td>199752</td><td><strong>1.18e-09</strong></td><td>1.26e-09</td></tr></tbody></table>

<figure><img src="image_16.png" /><figcaption>Figure 5: Shifted Griewank. RCO and PSO are equivalent.</figcaption></figure>
<table><caption>Table 12: Computing time per iteration</caption><thead><tr><td>Dimension</td><td>Planes</td><td>Shifted Rastrigin</td></tr></thead><tbody><tr><td>1</td><td>2.40E-05</td><td>1.60E-05</td></tr><tr><td>2</td><td>3.20E-05</td><td>4.80E-05</td></tr><tr><td>3</td><td>4.50E-05</td><td>5.30E-05</td></tr><tr><td>4</td><td>6.09E-05</td><td>6.69E-05</td></tr><tr><td>5</td><td>8.17E-05</td><td>8.47E-05</td></tr><tr><td>6</td><td>1.05E-04</td><td>1.05E-04</td></tr><tr><td>7</td><td>1.37E-04</td><td>1.38E-04</td></tr><tr><td>8</td><td>1.80E-04</td><td>1.78E-04</td></tr><tr><td>9</td><td>2.57E-04</td><td>2.31E-04</td></tr><tr><td>10</td><td>3.91E-04</td><td>3.15E-04</td></tr><tr><td>11</td><td>5.15E-04</td><td>4.67E-04</td></tr><tr><td>12</td><td>1.04E-03</td><td>1.09E-03</td></tr></tbody></table>

<figure><img src="image_19.png" /></figure>

# 6 When it doesn't work

On some problems, even continuous and of low dimension, RCO doesn't work well. Let's consider for example the Frequency-Modulated Sound Waves (FM) from the CEC 2011 competition benchmark (Das and Suganthan 2011). The search space is $[-6.4, 6.35]^6$ and the minimum is zero. Even after 50064 evaluations the best final value is 24.07. This is because many intersection positions (21238) are outside the search space (i.e. unacceptable).

Note that in such case it may be a little bit beneficial to expand the search space, if possible. If it is now $[-500, 500]^6$ there are 17322 unacceptable intersections and the final best value is 17.42. But many classical stochastic algorithms find a far better solution, precisely because the stochasticity can cope with the highly cahotic landscape (see the cross section 6).

In some instances, even when problems are continuous and of low dimension, RCO doesn't perform well. For instance, let's examine the Frequency-Modulated Sound Waves from the CEC 2011 competition benchmark (Das and Suganthan 2011). The search space is $[-6.4, 6.35]^6$ and the
# 5 Complexities

To assess an algorithm, it's common practice to conduct a theoretical analysis of both its space complexity and time complexity. However, as explained in [2], this approach isn't always appropriate, and a more reliable method is to estimate the actual memory usage and computational time.

In this context, there are two possible scenarios at each time step:

1. Solving a linear system of $D$ equations to determine a new position with $D$ coordinates. The required space is $D^2 + D$, with a time complexity of $O(D^3)$.

2. Computing a linear combination of $2^D$ positions, each with $D$ coordinates. The required space is $O(D \times 2^D)$, and the number of multiplications is $O(2^D)$.

However such a classical reasoning does not correspond to the behavior of the algorithm in practice. The point is that the second situation usually does not occur very often except at the end of the process. Moreover it depends on the landscape of the problem and also on the “budget” (the maximum acceptable number of evaluations).

Let's consider, for example, the computing time per evaluation for two problems: Planes and Shifted Rastrigin. We have already seen the Shifted Rastrigin and Planes is defined by

$$f(x) = \sum_{d=1}^{D} (x(d) - 3) \qquad (2)$$

Figure 6 depicts its landscape in two dimensions.

<figure><img src="image_17.png" /><figcaption>Figure 6: 2D Planes problem landscape.</figcaption></figure>

As we can see on table 13 the real computing time indeed grows exponentially, but more like $1.35^D$ than like $2^D$.
minimum value is zero. Despite 50064 evaluations, the best final value achieved is 24.07. This is primarily due to a considerable number of intersection positions (21238) lying outside the search space, rendering them unacceptable.

It's worth noting that in such cases, expanding the search space might provide some benefit if feasible. For example, if the search space is expanded to [-500, 500]⁶, there are 17322 unacceptable intersections, resulting in a final best value of 17.42. However, many classical stochastic algorithms tend to discover significantly better solutions. This is because stochasticity can effectively navigate a highly chaotic landscape, as is the one of FM (see Figure: 6).

<figure><img src="image_20.png" /><figcaption>Figure 6: Frequency-Modulated Sound Waves problem. Cross section on dimensions 2 and 5.</figcaption></figure>

# 7 Conclusion

Even in low dimensions, evaluating the objective function can sometimes be very expensive. In such cases, RCO may be useful as it may propose a sure solution after a few evaluations. However, there are some drawbacks:
<table><caption>Table 13: Computing time per iteration</caption><thead><tr><td>Dimension</td><td>Planes</td><td>Shifted Rastrigin</td></tr></thead><tbody><tr><td>1</td><td>2.40E-05</td><td>1.60E-05</td></tr><tr><td>2</td><td>3.20E-05</td><td>4.80E-05</td></tr><tr><td>3</td><td>4.50E-05</td><td>5.30E-05</td></tr><tr><td>4</td><td>6.09E-05</td><td>6.69E-05</td></tr><tr><td>5</td><td>8.17E-05</td><td>8.47E-05</td></tr><tr><td>6</td><td>1.05E-04</td><td>1.05E-04</td></tr><tr><td>7</td><td>1.37E-04</td><td>1.38E-04</td></tr><tr><td>8</td><td>1.80E-04</td><td>1.78E-04</td></tr><tr><td>9</td><td>2.57E-04</td><td>2.31E-04</td></tr><tr><td>10</td><td>3.91E-04</td><td>3.15E-04</td></tr><tr><td>11</td><td>5.15E-04</td><td>4.67E-04</td></tr><tr><td>12</td><td>1.04E-03</td><td>1.09E-03</td></tr></tbody></table>

<figure><img src="image_18.png" /></figure>

# 6 When it doesn't work

In some instances, even when problems are continuous and of low dimension, RCO doesn't perform well. For instance, let's examine the Frequency Modulated Sound Waves from the CEC 2011 competition benchmark ([5]). The search space is $[-6.4, 6.35]^6$ and the minimum value is zero. Despite 50064 evaluations, the best final value achieved is 24.07. This is primarily due to a considerable number of intersection positions (21238) lying outside the search space, rendering them unacceptable.

It's worth noting that in such cases, expanding the search space might provide some benefit if feasible. For example, if the search space is expanded to $[-500, 500]^6$, there are 17322 unacceptable intersections, resulting in a final best value of 17.42. However, many classical stochastic algorithms tend to discover significantly better solutions. This is because stochasticity can effectively navigate a highly chaotic landscape, as is the one of FMSW (see Figure 7). For example the slightly improved PSO already mentioned needs just two runs to find $2.218602 \times 10^{-27}$.
* It requires a reasonably good lower bound. On the other hand, compared to most other methods, this is the only user-defined parameter. In practice, for real problems, such a bound is often known.

* It does not work well on some problems, even of low dimension.

* It does not perform well on discrete problems, particularly when all variables are discrete. However, it still seems to be usable when only some variables are discrete.

* Its computing time per iteration grows exponentially with the dimension of the search space. It does not grow as quickly as theoretically expected, but nonetheless, in practice, it can be difficult to use on a laptop for high-dimensional problems.

# 8 Appendix

## 8.1 A bit of geometry

For the "Ruler & Compass" construction in dimension 1 and as seen on the figure 2 we have to be able to define two segments **a** and **b** so that $\frac{\mathbf{a}}{\mathbf{b}} = \mathbf{A}/\mathbf{B}$. The method is more than 2600 years old (Thales of Miletus) so you may have forgotten it ...

The figure 7 shows how to do.

To perform the "Ruler & Compass" construction in one dimension, as depicted in Figure 2, it is necessary to define two segments, denoted as **a** and **b**, such that $\frac{\mathbf{a}}{\mathbf{b}} = \mathbf{A}/\mathbf{B}$. This method, dating back over 2600 years to Thales of Miletus, may have slipped from your memory over time...

Figure 7 illustrates the procedure.

<figure><img src="image_21.png" /><figcaption>Figure 7: Ruler and Compass method to define proportional segments.</figcaption></figure>

## 8.2 Sensitivity

There is only one user-defined parameter: the lower bound. Therefore, it is crucial to understand the algorithm's sensitivity to this parameter. As depicted in the figures 8, we observe that efficiency
<figure><img src="image_19.png" /><figcaption>Figure 7: Frequency Modulated Sound Waves problem. Cross section on dimensions 5 and 6.</figcaption></figure>

# 7 Conclusions and future works

Even in low dimensions, evaluating the objective function can sometimes be very costly. In such cases, RCO may be beneficial as it may propose a good solution after only a few evaluations. However, there are some disadvantages:

* It requires a reasonably good lower bound. That said, compared to most other methods, this is the only user-defined parameter. In practice, for real-world problems, such a bound is often known. Moreover RCO is able to automatically define an adaptive lower bound (see Appendix 8.5).

* It performs poorly on some problems, even in low dimensions, when the landscape is highly chaotic.

* It does not perform well on discrete problems, particularly when all variables are discrete. However, it still appears to be usable when only some of the variables are discrete.
is indeed sensitive to changes in the lower bound. Though the sensitivity is not extreme, it could still be beneficial to experiment with various values. It's worth noting that setting the lower bound to the exact minimum value, if known, is not the optimal choice.

<figure><img src="image_22.png" /><figcaption>Figure 8: Efficiency is somewhat sensitive to the user-defined lower bound, but not excessively so.</figcaption></figure>
* Its computation time per iteration increases exponentially with the dimension of the search space. Although it does not grow as quickly as theoretically predicted, in practice, it can still be challenging to use on a laptop for high-dimensional problems.

* As with many initial presentations of new optimization algorithms, such as Ant Colony Optimization (1991), Particle Swarm Optimization (1995), and Differential Evolution (1995), a formal convergence analysis has not yet been provided. Specifically, while it is clear that there is no "explosion" effect, as seen in the original PSO version ([1]), it would be beneficial to more precisely identify in which cases stagnation occurs and how it can be avoided. There is, in fact, an experimental RCO version that attempts to address the issue of stagnation, though it is not particularly convincing. For instance, in the case of the 10D Rosenbrock function (see Table 12), it finds a value of 2.3 instead of 8.996, which is still significantly worse than PSO ($5.19 \times 10^{-8}$).

# 8 Appendix

## 8.1 A bit of geometry

To perform the “Ruler & Compass” construction in one dimension, as depicted in Figure 3, it is necessary to define two segments, denoted as **a** and **b**, such that $\mathbf{a}/\mathbf{b} = \mathbf{A}/\mathbf{B}$. This method, dating back over 2600 years to Thales of Miletus, may have slipped from your memory over time...

Figure 8 illustrates the procedure.

<figure><img src="image_20.png" /><figcaption>Figure 8: Ruler & Compass method to define proportional segments.</figcaption></figure>

## 8.2 Barycenter

More generally we define the barycenter $xNew$ of the $2^D$ current positions $x$ thanks to the values of the function $f$ on these positions. Here is the Octave/Matlab© code of this calculation:
A clear example of sensitivity can be observed when attempting to solve the optimal control problem of the bifunctional catalyst blend in the CEC 2011 competition (Das and Suganthan 2011).

As depicted in Figure 9, even a slight modification to the lower bound leads to significantly different behavior due to the very small slope. Conversely, due to the same reason, a position far from the optimal one already exhibits a value close to the minimum.

<figure><img src="image_23.png" /><figcaption>Figure 9: CEC 2011 bifunctional catalyst blend optimal control problem. 12 evaluations.</figcaption></figure>

If you are unsure about the lower bound, you can utilize an adaptive one (refer to the parameter *adaptOption* in the source code). In fact, the results can occasionally be quite satisfactory, as evidenced by Table 13. Thus, it may be worthwhile to try different methods: user-defined and the adaptive ones. In the source code, there are three potential formulas for adaptation, but they are rather rudimentary. Certainly, there is room for improvement.

<table><caption>Table 13: Adaptive lower bound may improve the efficiency ... or the contrary!</caption><thead><tr><th></th><th>Dimension</th><th>Evaluations</th><th>Lower bound</th><th>Non adaptive</th><th>Adaptive 3 coeff. 0.5</th><th>Adaptive 3 coeff. 0.1</th></tr></thead><tbody><tr><td>Six Hump Camel Back</td><td>2</td><td>104</td><td>-1.1</td><td>-1.031473</td><td><strong>-1.03157</strong></td><td>-0.947</td></tr><tr><td rowspan="2">Shifted Rastrigin</td><td>3</td><td>1508</td><td>-0.1</td><td>0.002548</td><td>0.06729</td><td><strong>3.27E-05</strong></td></tr><tr><td>4</td><td>50016</td><td>-0.1</td><td>0.00173</td><td><strong>1.78E-15</strong></td><td><strong>0</strong></td></tr><tr><td rowspan="2">Rosenbrock</td><td>3</td><td>1508</td><td>-0.1</td><td>0.001</td><td>0.1866</td><td>0.238</td></tr><tr><td>4</td><td>10016</td><td>-0.1</td><td>0.003756</td><td><strong>7.52E-22</strong></td><td><strong>8.94E-24</strong></td></tr><tr><td>Pressure Vessel</td><td>4</td><td>6076</td><td>6000</td><td>6592</td><td>6955</td><td>1.29E+05</td></tr><tr><td>Gear Train</td><td>4</td><td>1753</td><td>0</td><td>6.65E-09</td><td>1.38E-06</td><td><strong>6.60E-10</strong></td></tr></tbody></table>
```R
% Barycenter
minf=min(f);
if minf>=0 w=sum(f)-f; else w=sum(f-minf)-(f-minf); end
if w>0 xNew=sum(w'.*x)/sum(w);
else xNew=mean(x); end % Particular case, w may be null
```

Note: Of course, other formulas are possible—particularly those that use the differences between the $f$-values and the lower bound—but this one seems to produce good results.

## 8.3 Number of minima

In black box optimization of a non-discrete function it is theoretically impossible to know the number of (local) minima, because the set of positions contains an infinity of points. On a digital computer it is theoretically possible by exhaustive search for, because the $\varepsilon$-machine, the number of positions is finite. However it is huge (typically $2^{53D}$ for a definition space of dimension $D$).

So, we have to make a compromise: an exhaustive search on a grid whose spacing is not too large (it may generate many fake minima) and not too small (too much computing time, numerical instability).

Then the following algorithm for a function $f$ may be useful:

Define a grid with $N$ positions along each dimension
For each position $x$ evaluate $f(x)$, and also the values $f_i = f(x_i)$ of the neighbors $x_i$ (usually 2D except on the frontier).
If $f(x) < f_i$ for all $i$, then assume that $x$ is a (local) minimum.

This algorithm requires about $2^{N+1}$ evaluations but does not need a large memory size. We could evaluate just once each position and save the result, but then the memory size would be big. And also the computing time anyway: for each position we would have to find if it has been already evaluated.

When two neighbors on the grid have the same value the algorithm may assume they are both “minima”, as we can see on the figure 9.

By applying this algorithm we find that the landscape of sinCos15 has 28 minima in dimension one (and 818 in dimension two) including just one global minimum (see Figure 1). Six Hump Camel Back has six minima including two global ones. And of course Parabola is unimodal. For the Frequency-modulated Sound Waves, we find that there are at least 99,325 minima. Further refinement would be too time consuming.

## 8.4 Difficulty measures

Detailed descriptions, examples and references can be found in [3]. As explained in it the $\delta_{\neg NisB}$ measure is interesting, for it implicitly takes into account not only the number of local minima but also the sizes of their attraction basins. Using a tolerance threshold is usually less discriminant and the normalized roughness is too rudimentary.

As we can see, when considering the sizes of the attraction basins, the difficulty of FMSW is estimated to be less than that of Six Hump Camel Back. This is because, although there are many local minima, each basin is very small. Therefore, many stochastic algorithms can easily escape these minima.

## 8.5 About the lower bound

As mentioned the basic RCO needs a user-defined lower bound. Therefore, it is interesting to understand the algorithm's sensitivity to this parameter. As depicted in the figures 10, we observe that efficiency is indeed
8.3 More examples

For amusement, here are a few 1D and 2D figures illustrating how rapidly the algorithm approaches
the solution, thanks to steps defined purely by geometry. Initially, they are quite large and gradually
diminish in size as they approach the solution's position.

<figure><img src="image_24.png" /><figcaption>Figure 10: Multiparaboloid. 14 evaluations, final value 0.38 (real minimum 0).</figcaption></figure>
<figure><img src="image_21.png" /></figure>

Figure 9: Looking for minima. The algorithm may find fake minima when two neighbors (here 2 and 4) on the grid have the same value..

<table><caption>Table 14: Difficulty measures in [0,1]. The ones of Parabola Six Hump Camel Back are given for comparison. The tolerance thresholds are relative to the global minimum value. For sinCos15 it corresponds to a value at most equal to 0.9997, as used in the Introduction.</caption><thead><tr><th></th><th>sinCos15</th><th>Parabola</th><th>Six Hump Camel Back</th><th>Frequency Modulated Sound Waves</th></tr></thead><tbody><tr><td>Normalized roughness</td><td>0.9643</td><td>0</td><td>0.3333</td><td>&gt; 0.9999</td></tr><tr><td>Tolerance threshold 1.93 × 10<sup>-4</sup></td><td>0.9996</td><td>0.9972</td><td>0.99997</td><td>1.0</td></tr><tr><td>δ<sub>NisB</sub></td><td>0.4706</td><td>0.0711</td><td>0.6258</td><td>0.5263</td></tr></tbody></table>
<figure><img src="image_25.png" /><figcaption>(a) 12 evaluations, final value -1.899584 (real minimum -1.8996).</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>(b) 16 evaluations, final value -1.732.</figcaption></figure>
sensitive to changes in the lower bound. Though the sensitivity is not extreme, it could still be beneficial to experiment with various values. It's worth noting that setting the lower bound to the exact minimum value, if known, may be not the optimal choice.

<figure><img src="image_22.png" /><figcaption>Figure 10: Efficiency is somewhat sensitive to the user-defined lower bound, but not excessively so.</figcaption></figure>

A clear example of sensitivity can be observed when attempting to solve the optimal control problem of the Bifunctional Catalyst Blend in the CEC 2011 competition ([5]).

As depicted in Figure 11, even a slight modification to the lower bound leads to significantly different behavior due to the very small slope. Conversely, due to the same reason, a position far from the optimal one already exhibits a value close to the minimum.

<figure><img src="image_23.png" /><figcaption>Figure 11: CEC 2011 Bifunctional Catalyst Blend optimal control problem. 12 evaluations.</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>Figure 12: Branin. 44 evaluations, final value 0.3988 (real minimum 0.397887).</figcaption></figure>
If you are unsure about the lower bound, you can utilize an adaptive one. It is evaluated as follows after initialization and at each iteration.

```python
if min_f>0 lower_bound=min_f*coeff; else lower_bound=min_f*(2-coeff);
```

Where `min_f` is the minimum function value found so far and `coeff` an user-defined parameter in `[0, 1[].`

It appears that the results can occasionally be quite satisfactory, as evidenced by Table 15. Thus, it may be worthwhile to try different methods: user-defined and adaptive ones.

<table><caption>Table 15: Adaptive lower bound may improve the efficiency ... or the contrary!</caption><thead><tr><th></th><th>Dimension</th><th>Evaluations</th><th>Lower bound</th><th>Non adaptive</th><th>Adaptive coeff. 0.5</th><th>Adaptive coeff. 0.1</th></tr></thead><tbody><tr><td>Six Hump Camel Back</td><td>2</td><td>104</td><td>-1.1</td><td>-1.031473</td><td><strong>-1.03157</strong></td><td>-0.947</td></tr><tr><td rowspan="2">Shifted Rastrigin</td><td>3</td><td>1508</td><td>-0.1</td><td>0.002548</td><td>0.06729</td><td><strong>3.27E-05</strong></td></tr><tr><td>4</td><td>50016</td><td>-0.1</td><td>0.00173</td><td><strong>1.78E-15</strong></td><td><strong>0</strong></td></tr><tr><td rowspan="2">Rosenbrock</td><td>3</td><td>1508</td><td>-0.1</td><td>0.001</td><td>0.1866</td><td>0.238</td></tr><tr><td>4</td><td>10016</td><td>-0.1</td><td>0.003756</td><td><strong>7.52E-22</strong></td><td><strong>8.94E-24</strong></td></tr><tr><td>Pressure Vessel</td><td>4</td><td>6076</td><td>6000</td><td>6592</td><td>6955</td><td>1.29E+05</td></tr><tr><td>Gear Train</td><td>4</td><td>1753</td><td>0</td><td>6.65E-09</td><td>1.38E-06</td><td><strong>6.60E-10</strong></td></tr></tbody></table>

## 8.6 More examples

For amusement, here are a few 1D and 2D figures illustrating how rapidly the algorithm approaches the solution, thanks to steps defined purely by geometry. Initially, they are quite large and gradually diminish in size as they approach the solution's position.
<figure><img src="image_28.png" /><figcaption>Figure 13: Shifted Griewank. 14 evaluations, final value 0.34 (real minimum 0).</figcaption></figure>

## 8.4 Source codes

Written for Octave/Matlab®. Depending on the versions a few instructions (especially those for plotting) may not work with Octave.

Here is a typical run for Pressure Vessel:

```matlab
>> RCO(11,4,3000)
====================================================
Best after init 4.636519e+07
function 11, dimension 4, budget 3000
lowerBound 6000.000000
Evaluations 3074 (including 58 for discretisation)
noWayCount 12
outsideCount 341
fBest 6.615977e+03
on
  1.0000    0.5000    50.2674    97.8799
>>
```
<figure><img src="image_24.png" /><figcaption>Figure 12: Multiparaboloid. 14 evaluations, final value 0.38 (real minimum 0).</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Figure 13: Branin. 44 evaluations, final value 0.3988 (real minimum 0.397887).</figcaption></figure>
### 8.4.1 Problems

See also in the attached cec.zip file:

c_bifunc_data.mat

data6Bus.m

diffsolv.m

cec2011_TNEP.m

that are needed for some problems of the CEC 2011 competition.

```matlab
function [xMin, xMax, lowerBound, quantis]=get_func(func,D)
One=ones(1,D);
quantis=[];
switch func
    case -3 % abs(sin(x) + ax)
        xMin=zeros(1,D);
        xMax=10*One;
        lowerBound=-0.1;
    case -2 % Test
        xMin=zeros(1,D);
        xMax=10*One;
        lowerBound=-0.1;
    case {-1,0,3,4,5,8}
        xMin= zeros(1,D);
        xMax=10*One;
        lowerBound=-0.1;
    case 6
        xMax=10*One;
        xMin=-xMax;
        lowerBound=-0.1;
    case 7
        xMin=-One;
        xMax=One;
        lowerBound=-0.1;
    case 9 % Rosenbrock
        xMax=100*One;
        xMin=-xMax;
        lowerBound=-0.1;
        %lowerBound=0;
    case 10 % Six Hump Camel Back
        xMin=[-2,-1];
        xMax=-xMin;
```
<figure><img src="image_26.png" /><figcaption>(a) 12 evaluations, final value -1.899584 (real minimum -1.8996).</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>(b) 16 evaluations, final value -1.732.</figcaption></figure>

Figure 14: Combination of sinus.
lowerBound=-1.1;
%lowerBound=-1.031628453489877;

case 11 % Pressure Vessel (D=4)
    q=0.0625;
    quantis=[q,q,0,0];
    xMin=[q,q,10,10];
    xMax=[99,99,200,200];
    lowerBound=6000;
    %lowerBound=6059.714335048436;

case 12 % Gear Train (D=4)
    D=4;
    ones4=ones(1,4);
    xMin=12*ones4;
    xMax=60*ones4;
    quantis=[1,1,1,1];
    lowerBound=0;
    %lowerBound=2.700857e-12;

case 13 % Griewank
    xMax=100*One;
    xMin=-xMax;
    lowerBound=0;

case 14 % Sin Sin
    xMin=-2.7*One;
    xMax=7.5*One;
    lowerBound=-2*D;

case 15 % CEC 2011. Frequency-Modulated (FM) Sound Waves
    % COUNTER EXAMPLE
    D=6;
    xMin=-6.4*One;
    xMax=6.35*One;
    lowerBound=-0.0001;

case 16 % CEC 2011. The bifunctional catalyst blend optimal control problem
    D=1;
    xMin=0.6;
    xMax=0.9;
    lowerBound=-0.00001;

case 17 % CEC 2011. Transmission Network Expansion Planning
    % Difficult. Finds only 231 after 50000 evaluations
    % as most heuristics find the minimum 220.
<figure><img src="image_28.png" /><figcaption>Figure 15: Shifted Griewank. 14 evaluations, final value 0.34 (real minimum 0).</figcaption></figure>

## 8.7 Problem definitions

<table><thead><tr><th>Name</th><th>Definition</th><th>Search space</th><th>Minimum</th></tr></thead><tbody><tr><td>sinCos15</td><td>f = 2 + &prod;<sub>d=1</sub><sup>D</sup> (cos (x<sub>d</sub>) + |sin (x<sub>d</sub>) cos (x<sub>d</sub>)| sin (15x<sub>d</sub>) + x<sub>d</sub>/10)</td><td>[0, 10]<sup>D</sup></td><td>0.999507 for D = 1</td></tr><tr><td>Six Hump Camel Back</td><td>(4 - 2.1x<sub>1</sub><sup>2</sup> + x<sub>1</sub><sup>4</sup>/3) x<sub>1</sub><sup>2</sup> + x<sub>1</sub>x<sub>2</sub> + 4 (x<sub>2</sub><sup>2</sup> - 1) x<sub>2</sub><sup>2</sup></td><td>[-2, 1]<sup>2</sup></td><td>-1.031628453489877</td></tr><tr><td>Shifted Rastrigin</td><td><span style="font-size: 2em; vertical-align: middle;">{</span> u<sub>d</sub> = x<sub>d</sub> - 10d/(D + 1) f = 10D + &sum;<sub>d=1</sub><sup>D</sup> u<sub>d</sub><sup>2</sup> - 10 cos (2&pi;u<sub>d</sub>)<span style="font-size: 2em; vertical-align: middle;">}</span></td><td>[-10, 10]<sup>D</sup></td><td>0</td></tr><tr><td>Rosenbrock</td><td>f = &sum;<sub>d=1</sub><sup>D-1</sup> (100 (x<sub>d+1</sub> - x<sub>d</sub><sup>2</sup>)<sup>2</sup> + (1 - x<sub>d+1</sub><sup>2</sup>))</td><td>[-100, 100]<sup>D</sup></td><td>0</td></tr><tr><td>Pressure Vessel</td><td><span style="font-size: 2em; vertical-align: middle;">{</span> 0.0193x<sub>3</sub> - x<sub>1</sub> &le; 0 0.00954x<sub>3</sub> - x<sub>2</sub> &le; 0 1296000 - &pi;x<sub>3</sub><sup>2</sup> (x<sub>4</sub> + 4x<sub>3</sub>/3) &le; 0 f = 0.6224x<sub>1</sub>x<sub>3</sub>x<sub>4</sub> + 1.7781x<sub>2</sub>x<sub>3</sub><sup>2</sup> + x<sub>1</sub><sup>2</sup> (3.1661x<sub>4</sub> + 19.84x<sub>3</sub>)<span style="font-size: 2em; vertical-align: middle;">}</span></td><td><span style="font-size: 2em; vertical-align: middle;">{</span> q = 0.0625 {q, 2q, ..., 99}<sup>2</sup> &times; [10, 200]<sup>2</sup><span style="font-size: 2em; vertical-align: middle;">}</span></td><td>6059.714335048436</td></tr><tr><td>Gear Train</td><td>(f = 1.0/6.931 - x<sub>1</sub>x<sub>2</sub>/(x<sub>3</sub>x<sub>4</sub>))<sup>2</sup></td><td>{12, 60}<sup>4</sup></td><td>2.700857 &times; 10<sup>-12</sup></td></tr><tr><td>Frequency Modulated Sound Waves</td><td><span style="font-size: 2em; vertical-align: middle;">{</span> y<sub>t</sub> = &sum;<sub>t=1</sub><sup>100</sup> x<sub>1</sub> sin (x<sub>2</sub>t&theta; + x<sub>3</sub> sin (x<sub>4</sub>t&theta; + x<sub>5</sub> sin (x<sub>6</sub>t&theta;))) y<sub>0</sub> = &sum;<sub>t=1</sub><sup>100</sup> sin (5t&theta; - 1.5 sin (4.8t&theta; + 2 sin (4.9t&theta;))) f = (y<sub>t</sub> - y<sub>0</sub>)<sup>2</sup><span style="font-size: 2em; vertical-align: middle;">}</span></td><td>[-6.4, 6.35]<sup>6</sup></td><td>0</td></tr></tbody></table>
References

[1] M. Clerc and J. Kennedy. "The particle swarm - explosion, stability, and convergence in a multidimensional complex space." In: *IEEE Transactions on Evolutionary Computation* 6.1 (Feb. 2002), pp. 58-73.

[2] Maurice Clerc. "Iterative optimization -Complexity and Efficiency are not antinomic." Mar. 2024. DOI: 10.13140/RG.2.2.26318.43847. URL: https://hal.science/hal-04487869 (visited on 03/11/2024).

[3] Maurice Clerc. *Iterative Optimizers - Difficulty Measures and Benchmarks*. Wiley, 2019. URL: http://www.iste.co.uk/book.php?id=1477 (visited on 04/29/2019).

[4] Maurice Clerc. *PSO technical site*. URL: http://clerc.maurice.free.fr/pso/ (visited on 03/27/2024).

[5] Swagatam Das and P. N. Suganthan. *Problem Definitions and Evaluation Criteria for CEC 2011 Competition on Testing Evolutionary Algorithms on Real World Optimization Problems*. Tech. rep. 2011.

[6] David H. Wolpert and William G. Macready. "No Free Lunch Theorems for Optimization." In: *IEEE Transactions on Evolutionary Computation* 1.1 (1997), pp. 67-82.

[7] Xin-She Yang et al. "True global optimality of the pressure vessel design problem: a benchmark for bio-inspired optimisation algorithms." In: *International Journal of Bio-Inspired Computation* 5.6 (Jan. 2013), pp. 329-335. ISSN: 1758-0366. DOI: 10.1504/IJBIC.2013.058910. URL: http://www.inderscienceonline.com/doi/abs/10.1504/IJBIC.2013.058910 (visited on 11/05/2017).
plotPoints=false;
plotLine=true;
plotDash=false; % Only for D=1
plotRun=(plotPoints || plotLine || plotDash) && D<=2;
fprintf('\n====================================================================')
% Define the problem
[xMin, xMax, lowerBound, quantis]=get_func(func,D);
tStart=cputime;
% ----------------------------- Initialisation
% Define all "corners"
x=corners(xMin,xMax);
% Evaluate them
[nPoints, ~]=size(x);
for n=1:nPoints
    f(n)=fit(x(n,:),func);
end
% Save the best
[fBest, Ind]=min(f);
fprintf('\n Best after init %e',fBest)
xBest=x(Ind,:);
fBestQ=fBest;
if plotRun
    line=zeros(tMax+nPoints,D+1); % Just to speed up
    for n=1:nPoints
        line(n,1:D)=x(n,:); line(n,D+1)=f(n);
    end
end
% -----------------------------
tp=nPoints;
for t=1:tMax
    if adaptLBcoeff ~=0 % Adaptive lower bound
        if ~isempty(quantis) && qAlways
            fB=fBestQ;
        else
            fB=fBest;
        end
    switch adaptOption
        case 1 % Method 1
            if fB>=0
                LB=adaptLBcoeff*fB;
            else
                LB=(1+adaptLBcoeff)*fB;
            end
        case 2 % Method 2
    end
end
```pascal
% LB=(lowerBound+fB)*adaptLBcoeff;
LB=adaptLBcoeff*lowerBound +(1-adaptLBcoeff)*fB;

case 3 % Method 3
  if fB>0
    LB=fB*adaptLBcoeff;
  else
    % LB=(lowerBound+fB)*adaptLBcoeff;
    LB=adaptLBcoeff*lowerBound +(1-adaptLBcoeff)*fB;
end

end
else % User defined lower bound
  LB=lowerBound;
end
[xn,fn]=new(x,f,xMin,xMax,func,LB,t);

if plotRun % Keep positions for future plots
  tp=tp+1;
  line.tp,1:D)=xn; line.tp,D+1)=fn;
end

if ~isempty(quantis) && qAlways
  if fn<fBestQ % Keep the best solution after discretisation
    xQ=quantify(xn,quantis);
    fQ=fit(xQ,func);
    tQuantis=tQuantis+1;
  if fQ<fBestQ
    fBestQ=fQ;
    xBest=xQ;
  end
end
else
  if fn<fBest
    fBest=fn;
    xBest=xn;
  end
end

% Define another line (D=1) or plane (D=2) or hyperplane (D>2)
for p=1:nPoints-1
  x(p, :)=x(p+1, :);
  f(p)=f(p+1);
end
x(nPoints, :)=xn;
```
```matlab
fPlot(func,D,xMin,xMax,lowerBound);
hold on
end
if plotLine
    switch D
        case 1
            plot(line(:,1), line(:,2),'r-', 'LineWidth',2);
            %axis([xMin,xMax,lowerBound, max(line(:,2))]);
            hold on
        case 2
            plot3(line(:,1), line(:,2),line(:,3),'r-', 'LineWidth',2);
            hold on
        end
    end
if plotPoints
    switch D
        case 1
            plot(line(:,1), line(:,2),'r.', 'MarkerSize',12);
            %axis([xMin,xMax,lowerBound, max(line(:,2))]);
            hold on
        case 2
            plot3(line(:,1), line(:,2),line(:,3),'r.', 'MarkerSize',12);
            hold on
        end
    end
if plotDash && D==1
    [tMax,~]=size(line);
    for t=1:tMax-2
        X=[line(t,1),line(t+2,1)];
        Y=[line(t,2),lowerBound];
        plot(X,Y,'r--','LineWidth',1);
        hold on
    end
    for t=1:tMax % Vertical lines
        X=[line(t,1),line(t,1)];
        Y=[lowerBound,line(t,2)];
        plot(X,Y,'r--','LineWidth',1);
        hold on
    end
    %axis([xMin,xMax,lowerBound, max(line(:,2))]);
end
end
```
### 8.4.2 Algorithm

```matlab
function RCO(func,D,tMax) % Ruler and Compass Optimisation
    rng(123456789) % Only for test with function -2
    %{
        Maurice.Clerc@WriteMe.com
        2024-04-03
    }
```

A simple deterministic optimisation method.
On dimension 1 you just need a ruler and a compass to apply it.
And also on higher dimension D ... assuming you have
a D-ruler and a D-compass!
```matlab
    %}
    warning('off') % Don't display warnings (in particular for singular matrix)
    global noWayCount outsideCount
    noWayCount=0;
    outsideCount=0;
    qAlways=true; % For (partly) discrete problem
    % false (experimental) => discretisation only at the very end, on xBest
    % true => at each time step
    tQuantis=0; % For (partly) discrete problems. Number of improvements
    adaptLBcoeff=0.5; %0.1; % ~=0 => adaptive lower bound
    % 0 => only user defined lower bound (see get_func.m)
    adaptOption=2; % 1, 2 or 3 Adaptation method.It may be useful to
                    % try the three.
    % Possible plots for D=1 or 2
    plotLandscape=true;
    plotPoints=false;
    plotLine=true;
    plotDash=false; % Only for D=1
    plotRun=(plotPoints || plotLine || plotDash) && D<=2;
    fprintf('\n=============================================================')
    % Define the problem
    [xMin, xMax, lowerBound, quantis]=get_func(func,D);
    tStart=cputime;
    % ------------------------------- Initialisation
    % Define all "corners"
    x=corners(xMin,xMax);
    % Evaluate them
    [nPoints, ~]=size(x);
    for n=1:nPoints
        f(n)=fit(x(n,:),func);
    end
    % Save the best
    [fBest, Ind]=min(f);
    fprintf('\n Best after init %e',fBest)
    xBest=x(Ind,:);
```
fBestQ=fBest;
if plotRun
    line=zeros(tMax+nPoints,D+1); % Just to speed up
    for n=1:nPoints
        line(n,1:D)=x(n,:); line(n,D+1)=f(n);
    end
end
% ------------------------------------
tp=nPoints;
for t=1:tMax
    if adaptLBcoeff ~=0 % Adaptive lower bound
        if ~isempty(quantis) && qAlways
            fB=fBestQ;
        else
            fB=fBest;
    end

    switch adaptOption
    case 1 % Method 1
        if fB>=0
            LB=adaptLBcoeff*fB;
        else
            LB=(1+adaptLBcoeff)*fB;
        end
    end

    case 2 % Method 2
        % LB=(lowerBound+fB)*adaptLBcoeff;
        LB=adaptLBcoeff*lowerBound +(1-adaptLBcoeff)*fB;
    end

    case 3 % Method 3
        if fB>0
            LB=fB*adaptLBcoeff;
        else
            % LB=(lowerBound+fB)*adaptLBcoeff;
            LB=adaptLBcoeff*lowerBound +(1-adaptLBcoeff)*fB;
        end
    end

    else % User defined lower bound
        LB=lowerBound;
    end

    [xn,fn]=new(x,f,xMin,xMax,func,LB,t);

    if plotRun % Keep positions for future plots
        tp=tp+1;
        line.tp,1:D)=xn; line.tp,D+1)=fn;
    end
end

if ~isempty(quantis) && qAlways
    if fn<fBestQ % Keep the best solution after discretisation
        xQ=quantify(xn,quantis);
        fQ=fit(xQ,func);
        tQuantis=tQuantis+1;
        if fQ<fBestQ
            fBestQ=fQ;
            xBest=xQ;
        end
    end
end

else
    if fn<fBest
        fBest=fn;
        xBest=xn;
    end
end

% Define another line (D=1) or plane (D=2) or hyperplane (D>2)
for p=1:nPoints-1
    x(p,:)=x(p+1,:);
    f(p)=f(p+1);
end
x(nPoints,:)=xn;
f(nPoints)=fn;

end % end of the run

% ---------------------------- Results

if ~isempty(quantis) % Some problems may be (partly) discrete
    if ~qAlways
        xBest=quantify(xBest,quantis);
        fBest=fit(xBest,func);
        tQuantis=tQuantis+1;
    else
        fBest=fBestQ;
    end
end

tEnd=cputime;
fprintf(\n function %i, dimension %i, budget %i',func,D,tMax);
if adaptLBcoeff~=0

switch adaptOption
    case 1

fprintf(\n Adaptation 1')

fprintf(\n Totally adaptive lower bound, coefficient %f',adaptLBcoeff);
case 2
    fprintf('\n Adaptation 2');
    fprintf('\n Mix user defined lower bound %e + adaptive, coefficient %f',... lowerBound,adaptLBcoeff);
    fprintf('\n (no matter the sign of fBest)');

case 3
    fprintf('\n Adaptation 3');
    fprintf('\n Mix user defined lower bound %e + adaptive, coefficient %f',... lowerBound,adaptLBcoeff);
    fprintf('\n (depending on the sign of fBest)');

end
else
    fprintf('\n lowerBound %f',lowerBound);
end
fprintf('\n Evaluations %i', tMax+nPoints+tQuantis);
if tQuantis>0
    fprintf(' (including %i for discretisation)',tQuantis)
end
fprintf('\n noWayCount %i',noWayCount);
fprintf('\n outsideCount %i',outsideCount);
fprintf('\n fBest %e',fBest)
fprintf('\n on \n'); disp(xBest)
fprintf('\n CPU time %e \n',tEnd-tStart)
%========================================== Plots
if plotRun figure; end
if plotLandscape && D<=2
    fPlot(func,D,xMin,xMax,lowerBound);
    hold on
end
if plotLine
    switch D
        case 1
            plot(line(:,1), line(:,2),'r-', 'answidh',2);
            %axis([xMin,xMax,lowerBound, max(line(:,2))]);
            hold on
        case 2
            plot3(line(:,1), line(:,2),line(:,3),'r-', 'answidh',2);
            hold on
    end
end
if plotPoints
    switch D
        case 1
            plot(line(:,1), line(:,2),'r.', 'MarkerSize',12);
            %axis([xMin,xMax,lowerBound, max(line(:,2))]);
            hold on
    end
```matlab
case 2
    plot3(line(:,1), line(:,2),line(:,3),'r.', 'MarkerSize',12)
    hold on
end
end
if plotDash && D==1
    [tMax,~]=size(line);
    for t=1:tMax-2
        X=[line(t,1),line(t+2,1)];
        Y=[line(t,2),lowerBound];
        plot(X,Y,'r--','LineWidth',1);
        hold on
    end
    for t=1:tMax % Vertical lines
        X=[line(t,1),line(t,1)];
        Y=[lowerBound,line(t,2)];
        plot(X,Y,'r--','LineWidth',1);
        hold on
    end
    %axis([xMin,xMax,lowerBound, max(line(:,2))]);
end
end
```

```matlab
function C=corners(xMin,xMax)
D=numel(xMin);
C=combinRepet([0,1],D);
[nPoints,~]=size(C);
for n=1:nPoints
    for d=1:D
        if C(n,d)==0
            C(n,d)=xMin(d);
        else
            C(n,d)=xMax(d);
        end
    end
end
end
```

```matlab
function [xn,fn]=new(x,f,xMin,xMax,func,lowerBound,t)
% Construct the new point and evaluate it
global noWayCount outsideCount
```
% At the very beginning the intersection would necessarily be
% either "noWay" or "outside".
% We don't count it and don't check it.
if t==1
    xn out2in(f,x,lowerBound);
    fn=fit(xn,func);
return
end
% ----------------------
[~,D]=size(x);
% Define (hyper)planes (or lines if D=1)
for n=1:D
    P=[];
    F=[];
    for d=1:D+1
        rank=n+d-1;
        P=[P;x(rank,:)];
        F=[F;f(rank)];
    end
end
Plane(n,1:D+1)=plane(P,F);
end
P=[];
for n=1:D
    P=[P;Plane(n,:)];
end
[xn,noWay]=intersect(P,lowerBound);
if noWay % Weighted combination
    noWayCount=noWayCount+1; % Just for information
    xn=out2in(f,x,lowerBound);
else % Intersection but maybe outside the search space
    outside=false;
    for d=1:D
        if xn(d)<xMin(d) || xn(d)>xMax(d)
            outside=true;
            outsideCount=outsideCount+1;
            break;
        end
    end
    % Note: For simplicity we apply the same modification for
    % "noWay" and "outside". You may try something else.
    if outside
        xn=out2in(f,x,lowerBound);
    end
end
fn=fit(xn,func);
end

function xq=quantify(x,quantis)
D=numel(x);
for d=1:D
    q=quantis(d);
    if q>0
        xq(d)=q*nearest(x(d)/q);
    else
        xq(d)=x(d);
    end
end
end

function [xn,noWay]=intersect(P,lowerBound)
[nbPlans,n]=size(P);
D=n-1;
P3=[];
for np=1:nbPlans
    P3=[P3;lowerBound-P(np,n)];
end
M=[];
for np=1:nbPlans
    M=[M;P(np,1:D)];
end
%fprintf('\n M, P3 \n'); disp([M,P3])
xn=linsolve(M,P3);
% No solution
noWay=false;
check find(xn==Inf);
%fprintf('\n check1 '); disp(check)
if ~isempty(check)
    noWay=true;
else
    check find(xn==Inf);
    if ~isempty(check)
        noWay=true;
    else
        for d=1:D
            if isnan(xn(d))
                noWay=true; break;
            end
        end
    end
end
end
end

function C = combinRepet(list, K)
% Combinations with repetition
% We assume K >0
nL = length(list);
if K == 1 % Just the column vector
    C = list(:);
else
    [L{K:-1:1}] = ndgrid(1:nL);
    rL = reshape(cat(K+1, L{:}), [], K);
    C = list(rL);
end
end

function xn out2in(f,x,lowerBound)
[~,D] = size(x);
if min(f) <0
deltaf = fLowerBound;
else
deltaf = f;
end
w = sum deltaf - deltaf;
sw = sum(w);
for d = 1:D
xn(d) = sum(w' * x(:,d)) / sw;
end
end

function A=plane(P,F)
[n,~]=size(P);
M=[P,ones(n,1)]; % Add a column of 1
A=linsolve(M,F);'
end

function fPlot(func,D,xMin,xMax,lowerBound)
a=0.01; % For highly variable functions (like -2) use 0.001
step=a*max(xMax-xMin);
```matlab
switch D
    case 1 % D=1
        X=xMin:step:xMax;
        lX=lengh(X);
        for i=1:lX
            Y(i)=fit(X(i),func);
        end
        plot(X,Y)
        axis([xMin,xMax,lowerBound, max(Y)]);
        %axis([xMin,xMax]);
    case 2 % D=2
        X=xMin(1):step:xMax(1);
        Y=xMin(2):step:xMax(2);
        lX=lengh(X);
        lY=lengh(Y);
        [Xx,Yy]=meshgrid(X,Y);
        for lx=1:lX
            x=X(lx);
            for ly=1:ly
                y=Y(ly);
                Z(lx,ly)=fit([x,y],func);
            end
        end
        %surf(Xx,Yy,Z');
        mesh(Xx,Yy,Z');
       xlabel("x");
        ylabel("y");
        alpha 0.01 % For transparency
    end
end
```

## References

Clerc, Maurice (Mar. 2024a). "Iterative optimization -Complexity and Efficiency are not antinomic." DOI: 10.13140/RG.2.2.26318.43847. URL: https://hal.science/hal-04487869 (visited on 03/11/2024).

— (2024b). *PSO programs*. URL: http://clerc.maurice.free.fr/pso/ (visited on 03/27/2024).

Das, Swagatam and P. N. Suganthan (2011). *Problem Definitions and Evaluation Criteria for CEC 2011 Competition on Testing Evolutionary Algorithms on Real World Optimization Problems*. Tech. rep.
Yang, Xin-She et al. (Jan. 2013). "True global optimality of the pressure vessel design problem: a benchmark for bio-inspired optimisation algorithms." In: *International Journal of Bio-Inspired Computation* 5.6, pp. 329-335. ISSN: 1758-0366. DOI: 10.1504/IJBIC.2013.058910. URL: http://www.inderscienceonline.com/doi/abs/10.1504/IJBIC.2013.058910 (visited on 11/05/2017).