0 Daumen
951 Aufrufe

Hallo Community,

ich hätte gern ein Rechenbeispiel zum Thema Donwhill-Simplex Verfahren oder aber auch Nelder-Mead Simplex Verfahren.

Laut Aufgabenstellung ist das Minimum für folgende Funktion zu finden:

x^2 - 4x + y^2 - y - xy

Wäre über einen Ansatz sehr dankbar. Eine Literatur empfehlung nehme ich auch dankend entgegen!

Avatar von

2 Antworten

0 Daumen

Unbenannt.PNG

Unbenannt1.PNG

Text erkannt:

3D plot:

Avatar von 36 k
0 Daumen
Eine Literatur empfehlung nehme ich auch dankend entgegen!

Inzwischen gibt es einen Wikipedia-Artikel zum Thema, der auch IMHO so gut ist, dass man den Algorithmus leicht nachvollziehen kann (2013 war das noch nicht so!). An Hand der in der Frage gegebenen Funktion kann man auch gut zeigen, wie er abläuft. Ich habe dafür ein C++-Programm geschrieben (s.u.) und habe versucht den Ablauf mit Geoknecht3D zu visualisieren:

blob.png

Die Funktion \(f(x,y) = x^{2} - 4x + y^{2} - y - xy\) hat bekanntermaßen (s. Antwort von Moliets) ein Minimum bei \((3|\,2)\), was ich oben als Punkt markiert habe. Jeder Iterationsschritt besteht aus drei Punkten, die ich als Dreieck dargestellt habe, die sich im Verlauf immer weiter um dieses Optimum zuziehen. Im Laufe der Iterationen färben sich die Dreiecke von hellgrün bis dunkelrot.

Hier noch ein Beispiel ...

blob.png

... bei dem ich das 'Startdreieck' abseits des Optimums gesetzt habe. Man kann sehr schön sehen, wie sich der Algorithmus Richtung Optimum bewegt und es schließlich auch 'einkreist'.

Das Programm:

// --   Beispiel für den Downhill Simplex Algorithmus in C++
#include <algorithm>
#include <numeric>
#include <iostream>
#include <vector>

using xi_type = CVector;
// -- liefert diese Funktion 'true', so bricht der Algorithmus ab
bool genau_genug( const std::vector< xi_type >& xi )
{
  using std::abs;
  auto u = xi[1] - xi[0];
  auto v = xi[2] - xi[0];
  auto w = xi[1] - xi[2];
  return abs(u)*abs(v)*abs(w)/(2*abs(u % v)) <= 0.1; // -> Umkreisradius <= 1/10
}

int main()
{
  using namespace std;
  auto const alpha = 1.;
  auto const beta = 0.5;
  auto const gamma = 2.;
  auto const sigma = 0.5;

  // -- die zu optimierende Funktion aus https://www.mathelounge.de/54579
  auto const func = []( const xi_type& p ) {
      double x = X(p); double y = Y(p);
      return x*x -x*y + y*y - 4*x - y;
  };

  // == das Startdreieck
  vector< xi_type > xi = {xi_type(0,0), xi_type(5,1), xi_type(1,5) };
//    vector< xi_type > xi = {xi_type(-4,-4), xi_type(-4,-4.5), xi_type(-4.5,-4.2) };
  auto const N = size_t(xi.size()-1);

  cout.precision(4);
  auto farbe = CVector(0xb,0xf,0xb); // ... alles für den Farbverlauf
  auto const t = 0.2;
  auto Z_factor = 0.2;

  auto const ist_besser = [func]( const xi_type& a, const xi_type& b ) {
          return func(a) < func(b);  // Vergleich: a ist besser als b, wenn func(a) < func(b)
      };
  int k=0;
  for( ; !genau_genug(xi); ++k, farbe = t*CVector(0xf,0x4,0x4) + (1-t)*farbe)
  {
      cout << "dreieck("; // Ausgabe im Format für Geoknecht3D
      for(auto x: xi)
          cout << X(x) << "|" << Y(x) << "|" << func(x)*Z_factor << " ";
      cout << ")" << gz::farbe(farbe) << endl;

// sortiere die Punkte: x[0] ist der beste, x[N] der schlechteste
      sort( begin(xi), end(xi), ist_besser );
// m: Mittelpunkt der 'besseren' Punkte
      auto m = accumulate( begin(xi), end(xi)-1, xi_type() ) / N; 
// reflektiere den schlechtesten Punkt xi[N] am Mittelpunkt m, -> r
      auto r = (1+alpha) * m - alpha * xi[N];
      if( ist_besser(r, xi[0]) )
      {  // wenn r beser als x[0], bestimme den expandierten Punkt e und ersetze x[N] durch
// den besseren von r und e
          auto e = (1+gamma)*m - gamma*xi[N];
          xi[N] = min( r, e, ist_besser );
          continue;
      }
      if( ist_besser(r, xi[N-1]) )
      {  // wenn r besser als der zweitschlechteste Punkt x[N-1] ist, ersetze x[N] durch r
          xi[N] = r;
          continue;
      }
      // sei h der bessere der beiden Punkte x[N] und r. Bestimme den
// kontrahierten Punkt c = (1+ß)m - ßh
      auto h = min( r, xi[N], ist_besser );
      auto c = (1+beta) * m - beta * h;
      if( ist_besser(c, xi[N]) )
      {  // wenn c besser als x[N] ersetze x[N] durch c
          xi[N] = c;
          continue;
      }
      for( auto& x: xi ) // komprimiere den Simplex in Richtung des besten Punktes x[0]
          x = sigma * xi[0] + (1. - sigma) * x;
  }
  cout << "punkt(3|2|" << func(xi_type(3,2))*Z_factor << ")" << endl;  // für Geoknecht3D
  auto opt = *min_element( begin(xi), end(xi), ist_besser );
  cout.precision(6);
  cout << "Optimum[" << k << "]: x=" << X(opt) << " y=" << Y(opt)
<< " f(x,y)=" << func(opt) << endl;
  return 0;
}
Avatar von 48 k

Ein anderes Problem?

Stell deine Frage

Willkommen bei der Mathelounge! Stell deine Frage einfach und kostenlos

x
Made by a lovely community