Minimierung einer Abweichung: Mein Solver vs. Excel Solver

KeineAhnung78

Hallo,

ich habe einen Solver mit Gradientenmethode in C++ geschrieben, der leider deutlich schlechter ist als der Solver in Excel. Ich frage mich allerdings, woher die hohe Abweichung kommt und nehme zunächst mal an, dass ich noch einen elementaren Fehler in meinem Solver habe.

Um den Rahmen hier nicht zu sprengen, nur kurz zur eigentlichen Problemstellung:
Ein Rechteckstab mit Stromverdrängung (Skin-Effekt) kann analytisch berechnet werden, es ergeben sich hyperbolische Funktionen für Widerstand und Induktivität. Das Ergebnis kann in einem Kettenbruch entwickelt werden, der dann auch Grundlage für ein elektrisches Ersatzschaltbild mit z.B. 4 Widerständen und 4 Induktivitäten ist. Die Größen ergeben sich durch die Faktoren im Kettenbruch, der allerdings unendlich viele Glieder hat. Daher handelt es sich um eine Näherung. An einigen diskreten Frequenzpunkten soll nun das Ersatznetzwerk möglichst gut die analytischen Ergebnisse approximieren.

Also:

Anzahl Eingangsparameter = 8 (4 Widerstände, 4 Induktivitäten)
Ziel = Minimierung der Fehlerabweichung (Abweichung wird über quadratischen Fehler berechnet)
Die Funktion, die aus den Eingangsparametern die Fehlerabweichung berechnet kann nicht analytisch beschrieben werden!

Nun ja, Excel löst das Problem recht schnell und zuverlässig...mein Solver dagegen nicht .

Ich möchte kurz mein Vorgehen schildern und hoffe, dass dabei schon ein Vorgehensfehler erkennbar wird.
Als Startwerte der Optimierung verwende ich die ersten 8 Glieder des Kettenbruchs.

Gradienten berechnen:

const CVektor<double> CErsatznetzwerk::berechne_gradienten(const CVektor<double> & apunkt, const double & delta_x) 
{
	CVektor<double> grad;

	//aktuelle Werte merken
	vector<double> akt_R_zu_R0 = fR_zu_R0;
	vector<double> akt_L_zu_L0 = fL_zu_L0;

	//Punkt bertrachten
	//Fehler fuer delta_y/2 vor und nach dem punkt berechnen		
	for(unsigned i=0; i<apunkt.size(); i++)
	{
		if(i<fanz_R)
		{
			fR_zu_R0.at(i)=apunkt.at(i)+delta_x/2.;
			const double fehlersumme_plus = berechne_fehlersumme();
			fR_zu_R0.at(i)=apunkt.at(i)-delta_x/2.;
			grad.push_back((fehlersumme_plus-berechne_fehlersumme())/delta_x);
			//geaenderten Wert zuruecksetzen
			fR_zu_R0.at(i) = akt_R_zu_R0.at(i);
		}
		else
		{
			fL_zu_L0.at(i-fanz_R)=apunkt.at(i)+delta_x/2.;
			const double fehlersumme_plus=berechne_fehlersumme();
			fL_zu_L0.at(i-fanz_R)=apunkt.at(i)-delta_x/2.;
			grad.push_back((fehlersumme_plus-berechne_fehlersumme())/delta_x);
			//geaenderten Wert zuruecksetzen
			fL_zu_L0.at(i-fanz_R) = akt_L_zu_L0.at(i-fanz_R);
		}			
	}
	return grad; 
}

Die Richtung, in die nun optimiert werden soll, gibt der Gradient vor, aber nicht die Schrittweite. Also helfe ich mir mit ausprobieren:

const CVektor<double> CErsatznetzwerk::berechne_neuen_optimalpunkt(const CVektor<double> & apunkt, const CVektor<double> & agrad)
{

	CVektor<double> neuer_punkt;
	CVektor<double> bester_punkt=apunkt;

	//Fehlersumme des Ausgangspunktes merken
	const CVektor<double> alte_faktoren = gib_alle_faktoren();
	setze_alle_faktoren(apunkt);

	double akt_fehlersumme = berechne_fehlersumme();
	const double start_fehlersumme = akt_fehlersumme;
	double beste_fehlersumme = akt_fehlersumme;

	//Minimum ist gesucht --> Schrittweite < 0
	double schrittweite = -check(fdata->gib_double("Ersatznetzwerk","Start_Schrittweite"),1E-07,1.,"start_schrittweite");
	double schritt = schrittweite;
	const double abbruch_schrittweite = check(fdata->gib_double("Ersatznetzwerk","Abbruch_Schrittweite"),1E-07,1.,"abbruch_schrittweite");
	while(schrittweite<-abbruch_schrittweite)
	{	
		neuer_punkt = apunkt + schritt*agrad;

		setze_alle_faktoren(neuer_punkt);
		akt_fehlersumme = berechne_fehlersumme();
		if(akt_fehlersumme>beste_fehlersumme) 
		{
			schritt-=schrittweite;
			schrittweite/=2.;
			schritt+=schrittweite;
		}
		else
		{
			bester_punkt=neuer_punkt;
			beste_fehlersumme=akt_fehlersumme;		
			schritt+=schrittweite;
		}		
	}
	setze_alle_faktoren(alte_faktoren);
	return bester_punkt;
}

Das ganze ist nun eingebettet in eine Schleife:

while(iZaehler++<max_iterationen)
{

		//Gradienten berechnen
		CVektor<double> grad = fersatznetzwerk->berechne_gradienten(akt_punkt, delta_x);

		//Optimale Schrittweite in Richtung des Gradienten bestimmen
		akt_punkt = fersatznetzwerk->berechne_neuen_optimalpunkt(akt_punkt, grad);

		if((letzter_punkt-akt_punkt).betrag()<abbruch_distanz) 
			break;

		letzter_punkt = akt_punkt;
}

Anmerkungen:

* Es werden noch keine Nebenbedingungen abgefangen, z.B. dass die Faktoren (d.h. die Werte der Widerstände und Induktivitäten) größer null sein müssen.

* Die Funktion berechne_fehlersumme habe ich ausgiebig getestet. Sie liefert identische Wert zu Excel.

* Der Excel-Solver hat eine Menge Einstellmöglichkeiten. Mir geht es gar nicht um die Perfektion auch nicht in meinem eigenen Solver, aber die Größenordnung der Optimierung sollte schon stimmen:
Quadratischer Fehler vor Optimierung = 0,6528
Nach meiner Optimierung = 0,4
Nach Optimerung mit Excel Solver = 0,0122

Da stimmt doch was nicht oder???

Vielen Dank.

Benutzt Excel denn auch die Gradientenmethode oder eine andere? Hast du die Methode irgendwo mal im Klartext?

KeineAhnung78

Bei Excel kann man die Methode wählen. Entweder Newton oder Gradient. Mit beiden Methoden erhalte ich dasselbe Ergebnis.

Was meinst du mit Klartext? Mein Solver im Klartext oder den Excel-Solver im Klartext? Zu Excel gibt´s leider kaum Infos.

Mein Vorgehen in Worten:
Gegeben: Eingangsgrößen Vektor (8-dim) punkt=(x1,...,x8)
Ziel Unbekannte Funktion f(punkt) = min

Solange noch besserer Punkt gefunden werden kann

Berechne Gradienten:
grad = (df/dx1,...,df/dx8)
numerisch df/dx1 = (f(x1-delta_x/2) - f(x1+delta_x/2))/delta_x usw.
Bestimme einen neuen Punkt in Richtung des Gradienten
neuer_punkt = startpunkt - grad*schrittweite
wobei ich keine Ahnung habe, wie ich die schrittweite bestimmen soll,
also probiere ich aus, z.B. schrittweite = 0.5
Berechne neues Minimum
Loop