/* polytope.c ポリトープ法：多変数関数の最小化 */
#include <stdio.h>
#include <math.h>
/*
#define MONITOR
#define MONITOR1
*/
#define M 5 /* パラメータ数（次元数） */
#define MAXITER 10000 /* 最大繰り返し数 */
#define ALPHA 1.0 /* 反射係数 */
#define BETA 2.0 /* 拡大係数 */
#define GAMMA 0.5 /* 縮小係数 */
#define RATIO 1000. /* scalefactor / tolerance */
static double vertex[M + 1][M] = { 0.0, 1.0, 6.0, 3.5, 0.2 };
static double scalefactor[M] = { 1.0, 1.0, 1.0, 1.0, 0.1 };
static double tolerance[M] = { 0.001, 0.001, 0.001, 0.001, 0.0001 };
/*
static double vertex[M + 1][M];
static double scalefactor[M] = { 1., 1., 1., 1. };
static double tolerance[M];
*/
double func(double parameter[]);
int polytope(double xbest[], double *fbest, int *iter);
/* 最小化する関数の定義 */
/* f(Ｍ,γ,ｈ,ｂ,ａ)=Σ[ｈγ^2 /{(Xi - Ｍ)^2 + γ^2} + ｂ + ａXi - Yi] ^ 2 */
#define N 11 /* データの数 */
static double x[N] = { -10., -8., -6., -4., -2., 0., 2., 4., 6., 8., 10. };
static double y[N] = { 2.127, 2.520, 2.629, 2.938, 3.414,
4.669, 8.014, 6.372, 4.596, 4.296, 4.291 };
double func(double parameter[])
{
int i;
double f, g2, dx, d;
double p0, p1, p2, p3, p4;
p0 = parameter[0];
p1 = parameter[1];
p2 = parameter[2];
p3 = parameter[3];
p4 = parameter[4];
#ifdef MONITOR1
printf("%15.9e%15.9e%15.9e%15.9e%15.9e\n", p0, p1, p2, p3, p4);
#endif
f = 0.0;
g2 = p1 * p1;
for(i = 0; i < N; i++)
{
dx = x[i] - p0;
d = p2 * g2 / (dx * dx + g2) + p3 + p4 * x[i] - y[i];
f += (d * d);
}
return f;
}
/* f(x1,x2,x3,x4) = (x1 - 25.)^4 + |x2^3| + (x2 - x3 + 50.)^2 + (x4 - 15)^2 */
/*
double func(double x[])
{
double x1, x2, x3, x4;
#ifdef MONITOR1
printf("%15.9e%15.9e%15.9e%15.9e%15.9e\n", x[0], x[1], x[2], x[3]);
#endif
x1 = x[0] - 25.;
x2 = x[1];
x3 = x[1] - x[2] + 50.;
x4 = x[3] - 15.;
return x1 * x1 * x1 * x1 + fabs(x2 * x2 * x2) + x3 * x3 + x4 * x4;
}
*/
int polytope(double xbest[], double *fbest, int *iter)
{
int flag, j, k, kbest, kworst1, kworst2;
double c, t, x, fnew, *p, *vbest, *vworst;
static double f[M + 1]; /* 各頂点での関数値 */
static double vcenter[M], vnew[M], vwork[M];
#ifdef MONITOR
for(j = 0; j <= M; j++)
{
for(k = 0; k < M; k++) printf("%f ", vertex[j][k]);
printf("\n");
}
#endif
c = 1. / (sqrt((double)M + 1.) + 2.);
for(j = 0, p = xbest; j < M; j++)
{
t = (vertex[0][j] = *p++) + c * scalefactor[j];
for(k = 1; k <= M; k++) vertex[k][j] = t;
vertex[j + 1][j] = vertex[0][j] + scalefactor[j];
}
kbest = kworst1 = kworst2 = 0;
f[0] = func(vertex[0]);
for(k = 1; k <= M; k++)
{
f[k] = func(vertex[k]);
if(f[k] < f[kbest]) kbest = k;
else if(f[k] >= f[kworst1])
{
kworst2 = kworst1;
kworst1 = k;
}
else if(f[k] > f[kworst2]) kworst2 = k;
}
vbest = vertex[kbest];
vworst = vertex[kworst1];
#ifdef MONITOR
printf("iter value\n");
printf(" 1 %-14g", f[kbest]);
#endif
*iter = 1;
while(f[kbest] != f[kworst1])
{
for(j = 0; j < M; j++)
if(fabs(vbest[j] - vworst[j]) > tolerance[j]) break;
if(j == M) break;
if(++(*iter) > MAXITER)
{
fprintf(stderr, "\n未収束！（MAXITER = %d)\n", MAXITER);
break;
}
for(j = 0; j < M; j++)
{
x = 0.0;
for(k = 0; k <= M; k++)
if(k != kworst1) x += vertex[k][j];
vcenter[j] = x / M;
vnew[j] = vcenter[j] + ALPHA * (vcenter[j] - vworst[j]);
}
fnew = func(vnew);
if(fnew < f[kbest]) /* vnew が新しい最良の点 */
{
for(j = 0; j < M; j++)
vworst[j] = vcenter[j] + BETA * (vnew[j] - vcenter[j]);
f[kworst1] = func(vworst);
if(f[kworst1] >= fnew) /* 反射 */
{
#ifdef MONITOR
printf("R");
#endif
for(j = 0; j < M; j++) vworst[j] = vnew[j];
f[kworst1] = fnew;
}
#ifdef MONITOR
else printf("E"); /* 拡大 */
#endif
kbest = kworst1;
vbest = vworst;
#ifdef MONITOR
printf("\n%5d %-14g", *iter, f[kbest]);
#endif
kworst1 = kworst2;
vworst = vertex[kworst1];
kworst2 = kbest;
for(k = 0; k <= M; k++)
if(k != kworst1 && f[k] > f[kworst2]) kworst2 = k;
}
else if(fnew <= f[kworst2]) /* vnew は新しい最悪の点ではない */
{
#ifdef MONITOR
printf("R");
#endif
/* 反射 */
for(j = 0; j < M; j++) vworst[j] = vnew[j];
f[kworst1] = fnew;
kworst1 = kworst2;
vworst = vertex[kworst1];
kworst2 = kbest;
for(k = 0; k <= M; k++)
if(k != kworst1 && f[k] > f[kworst2]) kworst2 = k;
}
else
{
if(fnew < f[kworst1])
{
for(j = 0; j < M; j++) vworst[j] = vnew[j];
f[kworst1] = fnew;
}
do
{
#ifdef MONITOR
printf("C");
#endif
/* 縮小 */
for(j = 0; j < M; j++)
{
vwork[j] = vworst[j];
vworst[j] += GAMMA * (vbest[j] - vworst[j]);
}
fnew = func(vworst);
flag = 1;
for(j = 0; j < M; j++)
{
if(vwork[j] != vworst[j])
{
flag = 0;
break;
}
}
if(flag) break;
} while(fnew >= f[kworst1]);
f[kworst1] = fnew;
if(fnew < f[kbest])
{
kbest = kworst1;
f[kbest] = fnew;
vbest = vworst;
#ifdef MONITOR
printf("\n%5d %-14g", *iter, fnew);
#endif
}
if(fnew < f[kworst2])
{
kworst1 = kworst2;
vworst = vertex[kworst1];
kworst2 = kbest;
for(k = 0; k <= M; k++)
if(k != kworst1 && f[k] > f[kworst2]) kworst2 = k;
}
}
}
for(j = 0; j < M; j++) xbest[j] = vbest[j];
*fbest = f[kbest];
return 1;
}
int main(void)
{
double xbest[M];
double fbest, fold = 1.e+75;
int i, j, iter;
for(i = 0; i < M; i++) xbest[i] = vertex[0][i];
for(i = 0; i < M; i++) tolerance[i] = scalefactor[i] / RATIO;
j = 0;
while(1)
{
polytope(xbest, &fbest, &iter);
if(fbest >= fold) break;
printf("\n*** Polytope-Method Step %d\n", ++j);
printf("繰返し： %4d\n", iter);
printf("最小値： %g\n", (fold = fbest));
printf("パラメータ：\n");
for(i = 0; i < M; i++) printf("%5d %g\n", i, xbest[i]);
for(i = 0; i < M; i++) scalefactor[i] = tolerance[i];
for(i = 0; i < M; i++) tolerance[i] /= RATIO;
}
return 1;
}
|
（トップに戻る）
戻る（メインメニュー）
戻る（その他メニュー）
注意事項