Statistische Optimierung, Teil 2
Im zweiten Teil dieser Artikelserie wird die Arbeitsweise des verbesserten Random-Walk-Verfahrens sowie dessen C-Source-Listing vorgestellt. Desweiteren werden einige Anregungen zum Experimentieren gegeben, um die Leistungsfähigkeit des C-Programmes evtl, noch weiter zu steigern.
Das verbesserte Random-Walk-Verfahren lehnt sich stark an die Methode des Simulated Annealing (Simulierte Temperung) [2], [3] an, welche Erkenntnisse aus physikalischen Vorgängen, die bei der Züchtung von einkristallinen Substanzen wirksam sind, auf Optimierungsprobleme überträgt. Bei der Temperung eines bis zum flüssigen Zustand erhitzten Stoffes kommt es darauf an, die Temperatur genügend langsam und dazu noch in Schritten zu verringern, wobei in jedem Temperaturschritt lange genug verharrt werden muß, um den Atomen des Materials Gelegenheit zu geben, den thermischen Gleichgewichtszustand einzunehmen. In jedem dieser Temperaturschritte ist die kinetische Energie der Atome statistisch d.h. zufällig verteilt und strebt dem Gleichgewichtszustand und damit dem Energieminimum zu. Aufgrund der statistischen Energieverteilung kann es zwischenzeitlich jedoch durchaus Vorkommen, daß die Atome eine höhere kinetische Energie aufweisen, als der Gleichgewichtszustand angibt. D.h. trotz langfristiger Senkung des Gesamtenergieniveaus kommt es immer wieder zu zwischenzeitlichen Energieerhöhungen. Durch Temperung kann man erreichen, daß das Material langfristig den Zustand des globalen Energieminimums annimmt und ein einkristallines Gitter ausbildet. Die kinetische Energie der Atome bildet hier quasi das Analogon zu unserem Zielfunktionswert, der ja ebenfalls minimiert werden soll.
Die erste prinzipielle Änderung der verwendeten Optimierungsstrategie gegenüber dem einfachen Random-Walk ist nach dem Vorbild der Temperung die Strategie, gelegentlich auch solche Parametervektoren zum Nominalvektor werden zu lassen, die einen größeren Zielfunktionswert liefern als der gegenwärtig gültige Nominalvektor. Damit soll verhindert werden, daß sich das Verfahren in einem lokalen Minimum festsetzt. Die Implementierung dieses „gelegentlichen“ Zulassens von Verschlechterungen ist als Metropolis-Algorithmus [3] bekannt und funktioniert folgendermaßen:
Sei delta = neuer Zielfunktionswert - Zielfunkionswert des alten Nominalvektors
Bei negativem delta, also einer Verringerung des neuen Zielfunktionswertes gegenüber jenem, der zum alten Nominalvektor gehört, wird der neue Parametervektor (wie beim einfachen Random-Walk) stets zum neuen Nominalvektor. Bei positivem delta dagegen wird die Bewertungsfunktion
B = exp(-(delta*ß/Temp)), mit ß und Temp jeweils >0 (1)
gebildet, welche für die angegebenen Wertebereiche von delta, ß und Temp einen Wert zwischen 0 und 1 liefert. Nun wird mittels eines Zufallsgenerators eine im Intervall [0,1] gleich verteilte Zufallszahl zuf erzeugt. Ist nun zuf < B, wird der betrachtete Parametervektor als neuer Nominalvektor akzeptiert, obwohl der alte Nominalvektor einen besseren Zielfunktionswert hatte. Ist zuf dagegen >= B, bleibt der alte Nominalvektor erhalten. In Bild 1 ist ein Beispiel zu sehen, in welchem zuf größer ist als der aktuelle Wert von B, weswegen hier der alte Nominalvektor erhalten bliebe. Man sieht am Verlauf der Exponentialfunktion aus Bild 1, daß die Chance, den alten Nominalvektor durch einen schlechteren Parametervektor zu ersetzen, mit wachsendem delta abnimmt, da B mit wachsendem delta kleiner wird und somit die Wahrscheinlichkeit abnimmt, daß zuf den Wert von B unterschreitet. Damit werden also sehr schlechte Parametervektoren mit besonders geringer Wahrscheinlichkeit zum neuen Nominalvektor.
Diese Wahrscheinlichkeit verringert sich weiter mit einer im Laufe der Optimierung fortwährend vorgenommenen Verringerung von Temp. Der Variablenname Temp soll an die langsam fallende Temperatur bei der Temperung erinnern. Temp wird im Programm immer dann um einen gewissen Faktor reduziert, wenn bei ein und demselben „Temperaturschritt“ eine genügend große Anzahl neuer Parametervektoren erzeugt worden ist. Besitzt der Parametervektor z.B. dim Vektorelemente, so ist dieser Punkt bei entweder 100*dim neuerzeugten Vektoren oder bei 10*dim neuerzeugten Vektoren, die gleichzeitig Nominalvektor wurden, erreicht. Diese Festlegung ist willkürlich, liefert in der Praxis aber recht zufriedenstellende Resultate. Beim Standes Programms werden aufgrund eines relativ großen Wertes von Temp also noch verhältnismäßig viele Verschlechterungen des Zielfunktionswertes akzeptiert, während bei weit fortgeschrittener Optimierung praktisch nur noch Verbesserungen im Zielfunktionswert zu einem neuen Nominalvektor führen. Die Variableß ist eine experimentell festzulegende Steuervariable und wurde in dem verwendeten Programm auf 0.1 gesetzt. Die zweite Verbesserung der neuen Optimierungsstrategie gegenüber dem einfachen Random-Walk-Verfahren besteht in einer automatischen Anpassung des Standardabweichungsvektors sigma[ ], der für jenen Zufallsgenerator benötigt wird, welcher die gaußverteilten Abweichungen (∆x und ∆y im letzten Beispiel aus Teil 1 dieses Artikels) gegenüber dem Nominalvektor erzeugt, der auf diese Weise die neuen Parametervektoren generiert. Die Anpassung wird über das sog. Akzeptanzverhältnis acratio gesteuert, welches sich aus dem Verhältnis der neuakzeptierten Nominalvektoren zu den insgesamt neuerzeugten Parametervektoren pro Temperaturschritt ergibt. Unterschreitet das Akzeptanzverhältnis einen bestimmten Wert level, werden alle Standardabweichungen in sigma[ ] halbiert.
Die Auswirkung der Standardabweichungshalbierung zeigt Bild 2. Man erkennt, daß die Gauß-Verteilungsfunktion bei halbierter Standardabweichung deutlich steiler verläuft, was dazu führt, daß die Zufallszahlen ∆x dichter um den Mittelwert 0 liegen.
Die im folgenden neuerzeugten Parametervektoren weichen damit nicht mehr so stark vom alten Nominalvektor ab. Auf diese Weise ist eine Feinoptimierung möglich, die genau dann benötigt wird, wenn das Programm bei fortgeschrittener Optimierung bereits Parametervektoren in der Nähe des gesuchten globalen Minimums erzeugt. Aufgrund dieser „Halbierungsstrategie“ muß bei der Initialisierung des Programms darauf geachtet werden, daß die Standardabweichungen in sigma[] nicht zu klein gewählt werden.
Das vorliegende C-Programm zur Parameteroptimierung wurde mit dem Public-Domain-C-Compiler „Sozobon“ entwickelt. Die weiterhin benötigte Mathematik-Bibliothek ist oft in den Public-Domain-Pools von Universitäten zu finden. Im Zweifelsfall muß man sich die benötigten Funktionen exp(), fabs(), floor(), sin() und sqrt() selbst schreiben. Hierfür geeigneter Sourcecode ist z.B. in [1] zu finden.
Schachteldimensionierungsproblem
4 -111.4
150 10.0
0.1 10.0
5.1 10.0
7.2 10.0
10.3 Beispiel für INPUT.DAT bei der Sinusapproximation
5 1.0
150 10.0
0.1 10.0
5.1 10.0
7.2 10.0
10.3 10.0
10.4 Beispiel für INPUT.DAT bei der Toleranzschemaeinpassung
Das Programm liest die Werte für dim, ntemp (Anzahl der Temperaturschritte) und beta sowie die Elemente des Startvektors xi[] und des Vektors sigma[] von der Eingabe-Datei INPUT.DAT. Auf diese Weise ist es ohne erneute Compilierung des Programmes möglich, das Programm mit verschiedenen Anfangsbedingungen und Parametereinstellungen starten zu lassen. Zur Protokollierung des Optimierungsverlaufes wird nach jedem Temperaturschritt der gegenwärtig gültige Nominalvektor xi[] sowie einige weitere Programmparameter auf dem Bildschirm ausgegeben. Sind ntemp Temperaturschritte bearbeitet worden, erfolgt eine letzte Werteausgabe, die diesmal jedoch nicht auf den Bildschirm, sondern in die Datei OUTPUT.DAT geschrieben wird.
Möchte man intensiv mit dem Programm experimentieren, empfiehlt es sich, noch mehr Parametereinstellungen, wie z. B. jene für redfctr, tfactr etc., über INPUT.DAT vornehmen zu lassen. Man wird feststellen, daß das Optimierungsziel von vielen ganz und gar unterschiedlichen Startkonfigurationen aus erreicht wird. Wichtig dabei ist, daß die Standardabweichungen in sigma[] nicht zu klein gewählt werden. Am günstigsten ist es selbstverständlich, wenn man den abzudeckenden Wertebereich der Parameter in xi[] bereits kennt, vor allem dann, wenn die Parameter stark unterschiedliche Größenordnungen besitzen. In einem solchen Fall wäre z.B. ein anspruchsvollerer Mechanismus als die „Halbierungsstrategie“ zur Anpassung der Standardabweichungen vorteilhaft, um die einzelnen Standardabweichungen unabhängig voneinander beeinflussen zu können. Ein weiterer entscheidender Punkt für das zufriedenstellende Funktionieren des Optimierungsprogrammes ist ein guter Zufallsgenerator. Der implementierte Zufallsgenerator zuf() liefert zwar recht anständige Resultate, dennoch lohnt sich u.U. ein Experimentieren mit anderen Zufallsgeneratoren.
Wird eine besonders gute Optimierung benötigt, sollte man an Stelle des Datentyps float auf den Datentyp double umsteigen (der vom Sozobon-Compiler allerdings nur unzureichend unterstützt wird). Die ohnehin schon recht langen Optimierungszeiten, die für die Optimierungsbeispiele der Sinusapproximation und der Toleranzschemaeinpassung im Stundenbereich liegen können, werden dadurch natürlich noch länger. Lange Rechenzeiten sind jedoch typisch für statistische Optimierungsverfahren und stellen den zu zahlenden Preis dar, mit welchem die universelle Verwendbarkeit und die einfache Implementierung dieser Verfahren erkauft werden müssen.
Dr. Rainer Storn
Literatur:
[1] Jamsa, K., Bibliothek der C-Routinen, McGraw-Hill, 1986.
[2] Bohachevsky, I.O., Johnson, M.E. und Stein, M.L., „Generalized Simulated Annealing for Function Optimization“, Technometrics, Aug. 1986, pp. 209-217.
[3] Kirkpatrick, S., Gelatt, C.D Jr. und Vecchi, M.P., "Optimization by Simulated Annealing Science, Mai 1983, pp. 671-680.
/************************************************/
/* */
/* Verbessertes Random-Walk-Verfahren zur */
/* statistischen Parameteroptimierung. */
/* */
/* Rainer Storn (c) 1992 MAXON Computer */
/* */
/************************************************/
# include <stdio.h>
# include <math.h>
# define maxdim 30
# define ngauss 200
# define pi 3.141592654
/*------------Globale Variablen-----------------*/
int dim;
float gauss[ngauss], ugauss[ngauss], zz[maxdim];
unsigned long seed1, seed2, seed3;
/*------------Funktionsdeklarationen-----------*/
float zuf();
float phi();
void table();
void ran1();
float zielf();
/*-------------Funktionsdefinitionen----------*/
float zuf()
/************************************************/
/* */
/* zuf() erzeugt eine im Intervall [0,1] */
/* gleichverteilte Zufallszahl. */
/* */
/************************************************/
{
unsigned long m1, m2, m3, a1, a2, a3, c1, c2, c3;
float x;
m1 = 259200;
m2 = 134456;
m3 = 243000;
a1 = 7141;
a2 = 8121;
a3 = 4561;
c1 = 54773;
c2 = 28411;
c3 = 51349;
seed1 = (seed1*a1 + c1)%m1;
seed2 = (seed2*a2 + c2)%m2;
seed3 = (seed3*a3 + c3)%m3;
if ( (float)seed3/(float)m3 > 0.5)
{
return((float)seed1/(float)m1);
}
else
{
return((float)seed2/(float)m2);
}
}/*------Ende von zuf ()-------------------------*/
float phi(x)
/************************************************/
/* */
/* Berechnung von phi(x)=0.5*(1-erf(x)) */
/* mit Hilfe der Simpson-Integration */
/* für positive x. phi(x) ist die Ver- */
/* teilungsfunktion der Gauss'schen */
/* Standard-Nozmalverteilung. */
/* */
/************************************************/
float x;
{
int n, i;
float am, b, g, h, h1, error, e1, s, x1, xn, y;
/*---Untergrenze 0, Obergrenze b von erf(x)------*/
b = fabs(x);
if (b < 1.e-9)
{
y = 0.5;
goto ende;
}
error = 1.e-9;
/*---Berechnung der Streifenzahl n---------------*/
am = b/2.0;
h1 = b/40.0;
e1 = 0.0; /*--Fehlerrechnung--*/
x1 = am - 2.0*h1;
e1 = e1 + exp(-x1*x1);
x1 = am - h1;
e1 = el - 4.0*exp(-x1*x1);
x1 = am;
e1 = e1 + 6.0*exp(-x1*x1);
x1 = am + h1;
e1 = e1 - 4.0*exp(-x1*x1);
x1 = am + 2.0*h1;
e1 = e1 + exp(-x1*x1);
e1 = fabs(e1)/(h1*h1*h1*h1);
xn = sqrt(sqrt(b*b*b*b*b*e1/error/180.0));
xn = 2.0*floor(xn/2.0 + 1.0);
n = (int)xn;
/*-----------Auswertung der Simpson-Formel------*/
s = 1.0 + exp(-b*b);
g = 4.0;
h = b/xn;
for (i=1; i<= (n-1); i=i+1)
{
x1 = (float)i;
x1 = x1*h;
s = s + g*exp(-x1*x1);
g = 6.0 - g;
}
/*-----------Bestimmung von phi(x) = 0.5(1+erf(x))--------*/
y = 0.5*(1.0 + (2.0*s*h/3.0)/1.772453851);
if (x < 0.0)
{
y = 1.0 - y;
}
ende: return(y);
}/*------Ende der Funktion phi(x)------------- ---*/
void table()
/************************************************/
/* */
/* table erstellt e. Tabelle für die Umkehrfunk-*/
/* tion von phi(x). */
/* */
/************************************************/
{
int mgn1, mgn2, mean, k, i;
float xphi, xeig, yinc, error;
xsig =4.0;
/*------Erzeugen einer Tabelle von phi(x) im Feld gauss[]-----*/
for (i=0; i<ngauss; i=i+1)
{
xphi = -xsig + 2.0*xsig*((float)i)/((float)(ngauss-1));
gauss[i] = phi(xphi);
}
/*------Erzeugen der Umkehrfunktionstabelle zu phi (x) .-------*/
/*------Hierzu wird eine Binärsuche und lineare Interpolation--*/
/*------verwendet.----------------------------*/
for (i=0; i<ngauss; i=i+1)
{
yinc = ((float)i)/((float)(ngauss-1));
mgn1 = 0;
mgn2 = ngauss;
mean = (mgn1 + mgn2)/2;
/*--------------------Binärsuche----------------*/
cont: k = mean;
error = gauss[k] - yinc;
if (error <= 0.0)
{
mgn1 = k;
}
else
{
mgn2 = k;
}
mean = (mgn1 + mgn2)/2;
if (k != mean) goto cont;
/*-------------------Interpolation--------------*/
if (error <= 0.0)
{
if (k >= (ngauss-1))
{
ugauss[i] = xsig;
}
else
{
ugauss[i] = 2.0*xsig*((float)k + (yinc-gauss[k])/(gauss[k+1]-gauss[k]))/(float)ngauss-xsig;
}
}
else
{
if (k <= 0)
{
ugauss[i] = - xsig;
}
else
{
ugauss[i] = 2.0*xsig*((float)(k-1) + (yinc-gauss[k-1])/
(gauss[k]-gauss[k-1]))/(float)ngauss-xsig;
}
}
}/*--------Ende der for-Schleife- ------------*/
}/*--------Ende der Funktion table())-------------*/
void ran1()
/************************************************/
/* */
/* ranl() erzeugt e. standard-normalverteilten */
/* Zufallsvektor zz[i] der Dimension dim. */
/* D.h. jed.Vektorelem. ist normalverteilt mit */
/* Mittelwert 0 und Standardabweichung 1. */
/* */
/************************************************/
{
int i, j1, j2;
float x, x1;
/*---Erzeugen eines Zufallsvektors mit gleichverteilten------*/
/*---Vektorelementen aus dem Intervall [0,1].--*/
for (i=0; i<dim; i=i+1)
{
zz[i] = zuf();
}
/*---Umwandlung des gleichverteilten Zufallsvektors in-------*/
/*---einen standard-normalverteilten Zufallsvektor.----------*/
for (i=0; i<dim; i=i+l)
{
x = (float)ngauss*zz[i];
x1 = floor(x);
j1 = (int)x; j2 = jl + 1;
if (j1 >= ngauss)
{
zz[i] = ugauss[ngauss-1];
}
else
{
zz[i] = ugauss[j1] + (ugauss[j2]-ugauss[j1])*(x-x1);
}
}
}/*-------Ende der Funktion ran1()--------*/
float zielf()
/************************************************/
/* */
/* Schachteldimensionierungsproblem */
/* —------------------------------- */
/* */
/* zielf() liefert den Zielfunktionswert */
/* y zurück, den es zu minimieren gilt. */
/* */
/* Benötigte globale Variablen: */
/* zz[maxdim], dim. */
/* */
/************************************************/
{
float area, v, a, b, c, x, y;
area = 4.0; /*---Willkürliche Festlegung der Gesamtfläche auf 4.0---*/
a = zz[0];
x = zz[1];
b = area/a;
c = 8.0;
if ((a<0)||(x<0)||(x>b/2.0)||(x>a/2.0))
{
y = 1000.0;
}
else
{
v = x*(a-2.0*x)*(b-2.0*x);
if (v < 0.0)
{
y = 1000.0;
}
else
{
y = c-v;
}
}
return(y);
}/*---------Ende von zielf()------------*/
/************************************************/
/* */
/* H A U P T P R O G R A M M */
/* */
/************************************************/
main()
{
int i, itemp, itry, idet, iacc, ntemp, maxtry, msucc;
int flag, isucc;
float beta, xi[maxdim], sigma[maxdim], xphi, z, level;
float acratio, tfactr, temp, best, actual, delta, xran;
float reduct;
/*---Initialisierungen-----------------*/
FILE *fpin;
FILE *fpout;
fpin = fopen("INPUT.DAT","r");
fpout = fopen("OUTPUT.DAT","w");
seed1 = 1723; /* Initialisierung von zuf() */
seed2 = 1541; /* Initialisierung von zuf() */
seed3 = 205; /* Initialisierung von zuf() */
tfactr = 0.9; /* Reduktionsfaktor für "Temperatur" */
reduct = 0.5; /* Reduktionsfaktor für sigma[] */
temp = 0.5; /* Anfangstemperatur */
best = 1.0e20; /* bester Anfangs-Zielfunktionswert */
/*---Lese Eingabedaten vom File INPUT.DAT--*/
fscanf (fpin, "%d",&dim);
/*---Anzahl Vektorelemente-------------*/
fscanf (fpin, "%d" , &ntemp) ;
/*---Anzahl Temperaturschritte---------*/
maxtry = 100*dim; /* Maximalzahl an neuen Parametervektoren */
/* pro Temperaturschritt */
msucc = 10*dim; /* Maximalzahl neugenerierter Nominalvektoren pro */
/* Temperaturschritt */
level = 1.0/(float)maxtry; /* Schwelle für Akzeptanzverh. */
fscanf(fpin, "%f",&beta) ;
for (i=0;i<dim;i=i+1)
{
fscanf(fpin,"%f",&xi[i]);
/*---Nominalwert, des Parametervektors */
}
for (i=0;i<dim;i=i+1)
{
fscanf(fpin,"%f",&sigma[i]);
/* Standardabweichungen */
}
/*---Erzeuge Gauss-Umwandlungstabelle---------*/
printf("********************************************\n");
printf("* *\n");
printf("* Verbessertes Random-Walk-Verfahren zur *\n”);
printf("* Parameteroptimierung. *\n");
printf("* *\n");
printf("* Von Rainer Storn *\n");
printf("* *\n");
printf("********************************************\n\n");
printf("Gausstabelle wird erzeugt\n");
printf("Rechner arbeitet\n\n");
table();
printf("Gausstabelle generiert\n");
printf("Optimierung beginnt\n\n");
/*---Starte Optimierungsroutine------------------*/
for (itemp=0; itemp<ntemp; itemp=itemp+1)
{
itry = 0;
isucc = 0;
iacc = 0;
idet = 0;
con1: itry = itry + 1;
/*---Erzeuge neuen Parametervektor---------------*/
ran1();
for (i=0; i<dim; i=i+1)
{
zz[i] = zz[i]*sigma[i];
zz[i] = zz[i] + xi[i];
}
/*---Werte Zielfunktion aus--------------------*/
actual = zielf();
delta = actual - best;
/*----Entscheide, ob d.neue Parametervektor als */
/*----neuer Nominalwert akzeptiert werden soll. */
if (delta < 0.0)
{
flag = 1;
iacc = iacc + 1;
}
else
{
idet = idet + 1;
if (temp >= 0.000001)
/*--Berücksichtigung der einge- */
{ /*--schrankten Genauigkeit des--*/
/*--Datentyps float-------------*/
if (zuf() < exp(-delta*beta/temp))
{
flag = 1;
iacc = iacc + 1;
}
else
{
flag = 0;
}
}
else
{
flag = 0;
}
}
/*----Falls flag=1 wird d.Parametervektor als */
/*----neuer Nominalwert akzeptiert.------------*/
if (flag == 1)
{
for (i=0; i<dim; i=i+1)
{
xi[i] = zz[i];
}
best = actual;
isucc = isucc + 1;
}
/*----Prüfe, ob bereits genug Zufallsänderungen erfolgt sind.--*/
if ((itry < maxtry) && (isucc < msucc))
goto con1;
/*----Falls ja, gehe zur nächstniedrigeren Temperatur.----------*/
temp = temp*tfactr;
acratio = (float)iacc/(float)itry;
if (acratio < level)
{
for (i=0; i<dim; i=i+1)
{
sigma[i] = sigma[i]*reduct;
}
}
/*----Gebe Optimierungsdaten des gegenwärtigen */
/*----Iterationsschrittes aus.-----------------*/
printf("\n");
for (i=0; i<dim; i=i+1)
{
printf("xi(%d) = %f\n",i,xi[i]);
}
for (i=0; i<dim; i=i+1)
{
printf("sigma(%d) = %f\n",i,sigma[i]);
}
printf("Temperaturschritt = %d\n", itemp);
printf("Temperatur = %f\n", temp);
printf("Zahl der Zufallsaenderungen = %d\n", itry);
printf("Akzeptanzverhaeltnis = %f\n", acratio);
printf("Zielfunktionswert = %f\n", best);
if (actual <= 1.0e-6) goto fin;
}/*----Ende for(itemp. . . »-Schleife----*/
/* Schreibe Eingabedaten auf File OUTPUT.DAT */
fin: fprintf(fpout,"\n");
for (i=0; i<dim; i=i+1)
{
fprintf(fpout,"xi(%d) = %f\n",i,xi[i]);
}
for (i=0; i<dim; i=i+1)
{
fprintf(fpout,"sigma(%d) = %f\n",i,sigma[i]);
}
fprintf(fpout,"Temperaturschritt = %d\n",itemp);
fprintf(fpout,"Temperatur = %f\n",temp);
fprintf(fpout,"Zahl der Zufallsaenderungen = %d\n",itry);
fprintf(fpout,"Akzeptanzverhaeltnis = %f\n",acratio);
fprintf(fpout,"Zielfunktionswert = %f\n",best);
fclose(fpin);
fclose(fpout);
} /*---------Ende von main()--------*/
float zielf()
/************************************************/
/* */
/* Schachteldimensionierungsproblem */
/* -------------------------------- */
/* */
/* zielf() liefert den Zielfunktionswert */
/* y zurück, den es zu minimieren gilt. */
/* */
/* Benötigte globale Variablen: */
/* zz[maxdim], dim. */
/* */
/************************************************/
{
float area, v, a, b, c, x, y;
area = 4.0; /*---Willkürliche Festlegung der Gesamtfläche auf 4.0----*/
a = zz[0];
x = zz[1];
b = area/a;
c = 8.0;
if ((a<0)||(x<0)||(x>b/2.0)||(x>a/2.0))
{
y = 1000.0;
}
else
{
v = x*(a-2.0*x)*(b-2.0*x);
if (v < 0.0)
{
y = 1000.0;
}
else
{
y = c-v;
}
}
return(y);
}/*---------Ende von zielf()--------*/
float zielf()
/************************************************/
/* */
/* Sinusapproximation */
/* ------------------ */
/* */
/* zielf() liefert den Zielfunktionswert */
/* y zurück, den es zu minimieren gilt. */
/* */
/* Benötigte globale Variablen: pi, */
/* zz[maxdim], dim. */
/* */
/************************************************/
{
int i, n;
float a, alpha, diff, y;
alpha = pi/100.0;
y = 0.0;
for (i=0; i < 5 0; i=i+1)
{
a = (float)i*alpha;
diff = 0.0;
for (n=dim-1; n>0; n=n-1)
{
diff = (diff + zz[n])*a;
}
diff = diff + zz[0] - sin(a);
y = y + diff*diff;
}
return(y);
}/*--------Ende von zielf()------*/
float zielf()
/************************************************/
/* */
/* Toleranzschemaeinpassung */
/* ------------------------ */
/* zielf() liefert den Zielfunktionswert */
/* y zurück, den es zu minimieren gilt. */
/* */
/* Benötigte globale Variablen: */
/* zz[maxdim], dim. */
/* */
/************************************************/
{
int i, n;
float delta1, delta2, delta3, x, y, maxdev;
/*-----Untersuche Bereich für Argument zwischen -1 und 1-------*/
delta1 = 0.0;
maxdev = 0.0;
for (i=0; i<=100; i=i+1)
{
y = 0.0;
x = -1.0 + (float)i/50;
for (n=dim-1; n>0; n=n-1)
{
y = (y + zz[n])*x;
}
y = y + zz [0];
if (fabs(y) > 1.0) maxdev = (1-fabs(y))*(1-fabs(y));
if (maxdev > delta1) delta1 = maxdev;
}
/*------Untersuche Argumentenwert +1.2------------*/
delta2 = 0.0;
x = 1.2;
y = 0.0;
for (n=dim-1; n>0; n=n-1)
{
y = (y + zz[n])*x;
}
y = y + zz[0];
if (y < 5.9) delta2 = (y-5.9)*(y-5.9);
/*-------Untersuche Argumentenwert -1.2----------*/
delta3 = 0.0;
x = -1.2;
y = 0.0;
for (n=dim-1; n>0; n=n-1)
{
y = (y + zz[n])*x;
}
y = y + zz[0];
if (y < 5.9) delta3 = (y-5.9)*(y-5.9);
return(sqrt(delta1+delta2+delta3));
}/*----------Ende von zielf(--------------------*/
Rainer Storn