Minimizzazione con il metodo del simplesso

Il metodo del simplesso è un metodo che consente di ricercare i minimi di una funzione n-dimensionale, a condizione che la funzione in esame sia sufficientemente regolare.

Una rappresentazione grafica del metodo del simplesso sulla funzione di Himmelblau

Il metodo valuta la funzione in tre punti, e a seconda dei valori applica le seguenti operazioni: riflessione, espansione, contrazione, e contrazione inversa.

Il metodo è iterativo, e raffina a ogni iterazione la bontà della soluzione. Particolare importanza viene data alla scelta delle condizioni di partenza: in alcuni casi, come nella funzione di Himmelblau, a condizioni iniziali differenti corrispondono minimi differenti. Nella nostra implementazione abbiamo voluto utilizzare le funzionalità più moderne di C++11, creando un'interfaccia che accetta funzionali. Tutti i minimi delle funzioni testate (Himmelblau, Goldstein-Price e Rosenbrock) vengono correttamente individuati dalla nostra implementazione.

Le quattro operazioni del simplesso

Risultati[modifica | modifica wikitesto]

I minimi trovati sono qui tabulati.

Funzione Init (x,y) Minimo (x,y) F(x_min)
Himmelblau 2, 1 3, 2 4.01633e-26
Himmelblau 2 4, -1 3.58443, -1.84813 3.87252e-27
Himmelblau 3 -4, -4 -3.77931, -3.28319 1.38828e-25
Rosenbrock 2, 1 1, 1 0
Goldstein-Price 2, 1 2.37168, 5.25664 3.20586e-36

Notiamo come l'algoritmo è fortemente dipendente dalla scelta dei parametri iniziali: per la funzione di Himmelblau a init diversi corrispondono minimi completamente differenti.

Nella nostra ottimizzazione il ciclo di ottimizzazione continua fino a che le varie operazioni applicate sul centroide non permettono uno scostamento molto vicino alla precisione macchina (ovvero 100 * std::numeric_limits<real>::epsilon()).

Errore simplesso ottimizzazione -plotex.png

L'andamento dello scostamento del centroide in funzione del numero di iterazioni mostra come questo algoritmo converga velocemente verso i minimi locali della funzione.

Listing del codice[modifica | modifica wikitesto]

Header[modifica | modifica wikitesto]

 1 #ifndef SIMPLEX_H
 2 #define SIMPLEX_H
 3 
 4 #include <vector>
 5 #include <functional>
 6 
 7 #include "definitions.h"
 8 
 9 
10 class Simplex
11 {
12 public:
13     static std::vector<real> do_simplex(std::vector<real> init, std::function<real (std::vector<real>)> func);
14 };
15 
16 #endif // SIMPLEX_H

Implementazione[modifica | modifica wikitesto]

  1 #include "simplex.h"
  2 #include "definitions.h"
  3 
  4 #include <vector>
  5 #include <iostream>
  6 #include <algorithm>
  7 #include <cmath>
  8 
  9 #include "definitions.h"
 10 
 11 real f1(std::vector<real> in)
 12 {
 13     real x = in[0];
 14     real y = in[1];
 15     return pow(x*x+y-11, 2)+pow(x+y*y-7, 2);
 16 }
 17 real f2(std::vector<real> in)
 18 {
 19     real x = in[0];
 20     real y = in[1];
 21     return exp(0.5*pow(x*x+y*y-25, 2)) * pow(sin(4*x-3*y), 4)+0.5*pow(2*x+y-10, 2);
 22 }
 23 real f3(std::vector<real> in)
 24 {
 25     return 100 * pow(in[1] - pow(in[0], 2), 2) + pow(1 - in[0], 2);
 26 }
 27 
 28 std::vector<real> Simplex::do_simplex(std::vector<real> init, std::function<real (std::vector<real>)> f)
 29 {
 30     real error_min = 1E2 * std::numeric_limits<real>::epsilon(); // criterio di arresto
 31 
 32     //a: riflessione  -> xr
 33     //b: espansione   -> xe
 34     //g: contrazione  -> xc
 35     //h: contrazione verso x1    
 36     const double a = 1.0, b = 1.0, g = 0.5, h = 0.5; //coefficienti moltiplicativi
 37     
 38     std::vector<real> xcentroid_old(2, 0); 
 39     std::vector<real> xcentroid_new(2, 0);
 40     std::vector<real> f_v(3, 0);         //f valutata ai vertici del simplesso
 41     
 42     int x1 = 0, xn = 0, xnp1 = 0;
 43     // x1:   f(x1) = min { f(x1), f(x2), f(x3) }
 44     // xnp1: f(xnp1) = max { f(x1), f(x2), f(x3) }
 45     // xn:   f(xn) < f(xnp1) && f(xn) > tutti gli altri f(x_i)
 46     
 47     int counter = 0; // contatore di iterazioni
 48     int max_iter = 1E5;
 49 
 50 
 51     std::vector<real> delta(init);
 52     delta[0] = delta[0]/20;
 53     delta[1] = delta[1]/20;
 54     
 55     // Costruzione del simplesso
 56     std::vector< std::vector<real> > x =  std::vector<std::vector<real> >();
 57 
 58     for (int i = 0; i < 2; ++i) {
 59         std::vector<real> tmp(init);
 60         tmp[i] +=  delta[i];
 61         x.push_back(tmp);
 62     }
 63     x.push_back(init);
 64     
 65     // Calcolo del centroide
 66     std::transform(init.begin(), init.end(),
 67                    xcentroid_old.begin(), std::bind2nd(std::multiplies<real>(), 3));
 68 
 69     // Ciclo di ottimizzazione
 70     for (counter = 0; counter < max_iter; ++counter) {
 71 
 72         for (int i = 0; i < 3; ++i) {
 73             f_v[i] = f(x[i]);
 74         }
 75 
 76         x1 = 0;
 77         xn = 0;
 78         xnp1 = 0;
 79 
 80         // Classifica x1, xn, xnp1
 81         for (int i = 0; i < 3; ++i) {
 82             if (f_v[i] < f_v[x1]) {
 83                 x1 = i;
 84             }
 85             if (f_v[i] > f_v[xnp1]) {
 86                 xnp1 = i;
 87             }
 88         }
 89         xn = x1;
 90         for (int i = 0; i < 3; ++i) {
 91             if (f_v[i] < f_v[xnp1] && f_v[i] > f_v[xn])
 92                 xn = i;
 93         }
 94 
 95         // calcolo di xg: centroide dei due vertici migliori (x1, xn)
 96         std::vector<real> x_g(2, 0);
 97         for (int i = 0; i < x.size(); ++i) {
 98             if (i != xnp1)
 99                 std::transform(x_g.begin(), x_g.end(), x[i].begin(), x_g.begin(), std::plus<real>());
100         }
101         std::transform(x_g.begin(), x_g.end(),
102                        x[xnp1].begin(), xcentroid_new.begin(), std::plus<real>());
103         x_g[0] = x_g[0]/2;
104         x_g[1] = x_g[1]/2;
105         
106         // verifica la condizione di arresto
107         real delta = 0;
108         for (int i = 0; i < 2; ++i) {
109             delta += fabs(xcentroid_old[i] - xcentroid_new[i]);
110         }
111 
112         if (delta / 2 < error_min) {
113             break;            //se siamo sufficientemente vicini, termina
114         } else {
115             xcentroid_old.swap(xcentroid_new);
116         }
117         
118         // riflessione:
119         std::vector<real> x_r(2, 0);
120         for (int i = 0; i < 2; ++i) {
121             x_r[i] = x_g[i] + a * (x_g[i] - x[xnp1][i]);
122         }
123         
124         real f_xr = f(x_r);
125 
126         if (f_v[x1] <= f_xr && f_xr <= f_v[xn]) {
127             x[xnp1][0] = x_r[0];
128             x[xnp1][1] = x_r[1];
129         }
130         
131         // espansione:
132         else if (f_xr < f_v[x1]) {
133             std::vector<real> x_e(2, 0);
134             for (int i = 0; i < 2; ++i) {
135                 x_e[i] = x_r[i] + b * (x_r[i] - x_g[i]);
136             }
137             if (f(x_e) < f_xr) {
138                 std::copy(x_e.begin(), x_e.end(), x[xnp1].begin());
139             } else {
140                 std::copy(x_r.begin(), x_r.end(), x[xnp1].begin());
141             }
142         }
143 
144         // contrazione:
145         else if (f_xr > f_v[xn]) {
146             std::vector<real> xc(2, 0);
147             for (int i = 0; i < 2; ++i)
148                 xc[i] = x_g[i] + g * (x[xnp1][i] - x_g[i]);
149             if (f(xc) < f_v[xnp1])
150                 std::copy(xc.begin(), xc.end(), x[xnp1].begin());
151             else {
152 
153                 for (int i = 0; i < x.size(); ++i) {
154                     if (i != x1) {
155                         x[i][0] = x[x1][0] + h * (x[i][0] - x[x1][0]);
156                         x[i][1] = x[x1][1] + h * (x[i][1] - x[x1][1]);
157                     }
158                 }
159 
160             }
161         }
162 
163     }
164     
165     if (counter == max_iter) {
166         std::cout << "Massimo numero di iterazioni raggiunta. Il problema non converge." << std::endl;
167     }
168     return x[x1];
169 }
170 
171 
172 #ifdef ESERCIZIO3
173 
174 int main()
175 {
176     std::vector<real> init(2, 1);
177     std::vector<real> result;
178 
179     auto func1 = [&] (std::vector<real> x) { return f1(x); };
180     auto func2 = [&] (std::vector<real> x) { return f2(x); };
181     auto func3 = [&] (std::vector<real> x) { return f3(x); };
182     
183     result = Simplex::do_simplex(init, func1);
184     std::cout  << "### f1(x_min) = " << func1(result) <<std::endl;
185     std::cout << result[0] << ", " <<  result[1] << std::endl;
186 
187     result = Simplex::do_simplex(init, func2);
188     std::cout  << "### f2(x_min) = " << func2(result) <<std::endl;
189     std::cout << result[0] << ", " <<  result[1] << std::endl;
190 
191     result = Simplex::do_simplex(init, func3);
192     std::cout  << "### f3(x_min) = " << func3(result) <<std::endl;
193     std::cout << result[0] << ", " <<  result[1]  << std::endl;
194 
195     std::cout << "Altri minimi di Himmelblau (diverso init):" << std::endl;
196     init[0] = 4;
197     init[1] = -1;
198     result = Simplex::do_simplex(init, func1);
199     std::cout << result[0] << ", " <<  result[1] << " f1(x_min) = " << func1(result) << std::endl;
200     init[0] = -4;
201     init[1] = -4;
202     result = Simplex::do_simplex(init, func1);
203     std::cout << result[0] << ", " <<  result[1] << " f1(x_min) = " << func1(result) << std::endl;
204 
205     return 0;
206 }
207 
208 #endif
 PrecedenteSuccessivo