下面是一个uncopt的class 用了4中不同的方法去求解 无约束优化问题
使用方法, 使用时候 构造一个新类继承class UnconstrainedOptim
并重写!.....
virtual void F() = 0; // f=F(x) 目标函数
virtual void G() = 0; // g=Grad F(x) 目标函数的倒数
有相应的选择..
DFP method, BFGS method,(拟牛顿法)
conjugate gradient methods(共饿梯度法)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
uncoptim.h
// =============================================================================
// UNConstrained OPTIMizatoin library.
// keywords: unconstrained optimization, secant methods, DFP method, BFGS method,
// conjugate gradient methods, FR method, PR method, line search algorithms.
// =============================================================================
#ifndef __UNCOPTIM_H
#define __UNCOPTIM_H
#include <math.h>
#define GOLDENSEARCH 0
#define BSSEARCH 1
#define SSEARCH 2
#define RSSEARCH 3
// ========= Preset Parameter for Line Search Algorithms ================
#define DEFAULT_ZERO 1E-8
#define DEFAULT_GOLDEN1 0.61803398875
#define DEFAULT_GOLDEN2 0.38196601125
#define DEFAULT_RHO 0.01
// 0 < RHO < 0.5
#define DEFAULT_SIGMA 0.2
// RHO < SIGMA < 1
#define DEFAULT_TAW1 9
// TAW1 > 1
#define DEFAULT_TAW2 0.1
// TAW2 < SIGMA , Recommeded
#define DEFAULT_TAW3 0.5
// 0 < TAW2 < TAW3 < 0.5
#define DEFAULT_T1 0.05
// ~0.1
#define DEFAULT_T2 0.05
// ~0.5
#define DEFAULT_T3 1.05
// ~1+T1
#define DEFAULT_T4 10
// 4-10
#define DEFAULT_T5 0.0025
// ~T1^2/(1-T1)
#define DEFAULT_T6 0.5
// ~0.5
inline double min(double a, double b) {
return a<b ? a:b;
}
inline double max(double a, double b) {
return a>b ? a:b;
}
inline double arg(double x, double y) {
double z;
if(fabs(x)>fabs(y)) {
z=atan(y/x);
if(x<0)
z+=3.1415926536;
}
else {
z=atan(x/y);
if(y>0)
z=1.5707963268-z;
else if(x>0)
z=-z-1.5707963268;
else
z=4.7123889804-z;
}
return z;
}
class UnconstrainedOptim {
public:
int n;
double *x;
double f;
double *g;
int cnls,cnf,cng;
int defaultLineSearchMethod;
private:
double flb, fopt;
double *xk, *sk, *gk;
int gknown;
const double ZERO;
const double GOLDEN1;
const double GOLDEN2;
const double RHO;
const double SIGMA;
const double TAW1;
const double TAW2;
const double TAW3;
const double T1;
const double T2;
const double T3;
const double T4;
const double T5;
const double T6;
public:
UnconstrainedOptim(int N, double FLB, double FOPT, int lineSearchMethod=BSSEARCH);
UnconstrainedOptim(int N, double FLB, double FOPT, int lineSearchMethod, double _ZERO,
double _GOLDEN1, double _GOLDEN2, double _RHO, double _SIGMA,
double _TAW1, double _TAW2, double _TAW3,
double _T1, double _T2, double _T3, double _T4, double _T5, double _T6);
virtual ~UnconstrainedOptim();
void resetCounter() {
cnls=cnf=cng=0;
}
virtual void F() = 0; // f=F(x)
virtual void G() = 0; // g=Grad F(x)
private:
double fx() {
F();
++cnf;
return f;
}
void gx() {
G();
++cng;
gknown=1;
}
void xh(double h) {
int i;
i=n;
while(--i>=0)
x[i]=xk[i]+h*sk[i];
gknown=0;
}
double gs(void) {
double w=0;
int i=n;
while(--i>=0)
w+=g[i]*sk[i];
return w;
}
// ==========================================================================
// ==========================================================================
// Q(x)=const+ px+ qx^2 ; interpolation 3 point
static void makeQ3p(double a, double fa, double b, double fb,
double c, double fc, double &p, double &q) {
double d1,d2,d3;
double z1,z2,z3;
d1=c-b; d2=c-a; d3=b-a;
z1=fa/d2/d3; z2=fb/d1/d3; z3=fc/d1/d2;
p=-((a+b)*z3-(a+c)*z2+(b+c)*z1);
q=z1-z2+z3;
}
// Q(x)=const+ px+ qx^2 ; interpolation point, slope, point
static void makeQpsp(double a, double fa, double ga, double b, double fb,
double &p, double &q) {
double d,z;
d=b-a;
z=(fb-fa)/d;
p=((a+b)*ga-2*a*z)/d;
q=(z-ga)/d;
}
// C(x)=const+ px+ qx^2+ rx^3 ; interpolation 2 point, 2 slope
static void makeC2ps(double a, double fa, double ga,
double b, double fb, double gb,
double &p, double &q, double &r) {
double d, z0, z1, z2, z3;
d=b-a; z0=a+b; z1=a+z0; z2=b+z0;
z3=(fb-fa)/d;
d=d*d;
p=(a*z2*gb+b*z1*ga-6*a*b*z3)/d;
q=-(z1*gb+z2*ga-3*z0*z3)/d;
r=(gb+ga-2*z3)/d;
}
// C(x)=const+ px+ qx^2+ rx^3 ; interpolation point, slope and 2 point
static void makeCps2p(double a, double fa, double ga,
double b, double fb, double c, double fc,
double &p, double &q, double &r) {
double d1, d2, d3;
double z1, z2, z3;
d1=c-b; d2=c-a; d3=b-a;
z1=fa/d2/d3; z2=fb/d1/d3; z3=fc/d1/d2;
r=(ga+(d2+d3)*z1-d2*z2+d3*z3)/d2/d3;
p=(a*b+b*c+c*a)*r-((a+b)*z3-(a+c)*z2+(b+c)*z1);
q=z1-z2+z3-(a+b+c)*r;
}
// Q(returned value)=Min Q(x)=const+ px+ qx^2 on [a,b]
static double cminQuad(double p, double q, double a, double b) {
double w;
if(fabs(q)<DEFAULT_ZERO) {
if(p<0)
return b;
else
return a;
}
else {
w=-p/2/q;
if(q>0) {
if(w<a)
return a;
else if(w>b)
return b;
else
return w;
}
else {
if(2*w<a+b)
return b;
else
return a;
}
}
}
// C(returned value)=Min C(x)=const+ px+ qx^2+ rx^3 on [a,b]
static double cminCubic(double p, double q, double r, double a, double b) {
double delta,w,z,cw,cz;
if(fabs(r)<DEFAULT_ZERO)
z=cminQuad(p,q,a,b);
else {
delta=q*q-3*p*r;
if( delta>0 ) {
w=(sqrt(delta)-q)/3/r;
z=b;
cz=(p+(q+r*z)*z)*z;
if( a<w && w<b ) {
cw=(p+(q+r*w)*w)*w;
if( cw<cz ) {
z=w;
cz=cw;
}
}
if((p+(q+r*a)*a)*a<cz)
z=a;
}
else { // delta<=0
if( r>0 )
z=a;
else
z=b;
}
}
return z;
}
static double chooseQ(double a, double b, double h1,double fh1, double gh1, double h2, double fh2) {
double p,q;
makeQpsp(h1,fh1,gh1,h2,fh2,p,q);
if( a<b )
return cminQuad(p,q,a,b);
else
return cminQuad(p,q,b,a);
}
static double chooseC(double a, double b, double h1, double fh1, double gh1, double h2, double fh2, double gh2) {
double p,q,r;
makeC2ps(h1,fh1,gh1,h2,fh2,gh2,p,q,r);
if( a<b )
return cminCubic(p,q,r,a,b);
else
return cminCubic(p,q,r,b,a);
}
static double chooseCC(double a, double b, double h1, double fh1, double gh1, double h2, double fh2, double h3, double fh3) {
double p,q,r;
if( fabs(h1-h2)<DEFAULT_ZERO || fabs(h2-h3)<DEFAULT_ZERO ) {
makeQpsp(h1,fh1,gh1,h3,fh3,p,q);
if( a<b )
return cminQuad(p,q,a,b);
else
return cminQuad(p,q,b,a);
}
else {
makeCps2p(h1,fh1,gh1,h2,fh2,h3,fh3,p,q,r);
if( a<b )
return cminCubic(p,q,r,a,b);
else
return cminCubic(p,q,r,b,a);
}
}
public:
void update() {
int i=n;
while(--i>=0) {
xk[i]=x[i];
gk[i]=g[i];
}
}
void goldenSearch(double &hk);
void BSsearch(double &h, double &g0, double &f_1);
void Ssearch(double &h, double &g0, double &f_1);
void RSsearch(double &h, double &g0, double &f_1);
int DFPoptimize() {
return DFPoptimize(defaultLineSearchMethod);
}
int BFGSoptimize() {
return BFGSoptimize(defaultLineSearchMethod);
}
int FRoptimize(int reset=0) {
return FRoptimize(defaultLineSearchMethod, reset);
}
int PRoptimize(int reset=0) {
return PRoptimize(defaultLineSearchMethod, reset);
}
int DFPoptimize(int lineSearchMethod);
int BFGSoptimize(int lineSearchMethod);
int FRoptimize(int lineSearchMethod, int reset);
int PRoptimize(int lineSearchMethod, int reset);
private:
int FRoptimizeWithReset(int lineSearchMethod);
int PRoptimizeWithReset(int lineSearchMethod);
};
#endif
!!!!!!!!!!!!!!!
uncopt.cpp
// =============================================================================
// UNConstrained OPTIMizatoin library.
// programmer: A. Katanforoush,
// last update: 7/2/2003,
// related web page: http://math.ipm.ac.ir/scc/PointsOnSpheres,
// keywords: unconstrained optimization, secant methods, DFP method, BFGS method,
// conjugate gradient methods, FR method, PR method, line search algorithms.
// =============================================================================
#include "uncoptim.h"
#include <math.h>
UnconstrainedOptim::UnconstrainedOptim(int N, double FLB, double FOPT, int lineSearchMethod)
:ZERO(DEFAULT_ZERO),GOLDEN1(DEFAULT_GOLDEN1),GOLDEN2(DEFAULT_GOLDEN2),
RHO(DEFAULT_RHO),SIGMA(DEFAULT_SIGMA),TAW1(DEFAULT_TAW1),TAW2(DEFAULT_TAW2),TAW3(DEFAULT_TAW3),
T1(DEFAULT_T1),T2(DEFAULT_T2),T3(DEFAULT_T3),T4(DEFAULT_T4),T5(DEFAULT_T5),T6(DEFAULT_T6) {
n=N;
x=new double[n];
g=new double[n];
flb=FLB;
fopt=FOPT;
xk=new double[n];
sk=new double[n];
gk=new double[n];
cnls=cnf=cng=0;
defaultLineSearchMethod=lineSearchMethod;
}
UnconstrainedOptim::UnconstrainedOptim(int N, double FLB, double FOPT, int lineSearchMethod, double _ZERO,
double _GOLDEN1, double _GOLDEN2, double _RHO, double _SIGMA,
double _TAW1, double _TAW2, double _TAW3,
double _T1, double _T2, double _T3, double _T4, double _T5, double _T6)
:ZERO(_ZERO),GOLDEN1(_GOLDEN1),GOLDEN2(_GOLDEN2),
RHO(_RHO),SIGMA(_SIGMA),TAW1(_TAW1),TAW2(_TAW2),TAW3(_TAW3),
T1(_T1),T2(_T2),T3(_T3),T4(_T4),T5(_T5),T6(_T6) {
n=N;
x=new double[n];
g=new double[n];
flb=FLB;
fopt=FOPT;
xk=new double[n];
sk=new double[n];
gk=new double[n];
cnls=cnf=cng=0;
defaultLineSearchMethod=lineSearchMethod;
}
UnconstrainedOptim::~UnconstrainedOptim() {
delete x;
delete g;
delete xk;
delete sk;
delete gk;
}
void UnconstrainedOptim::goldenSearch(double &hk) {
double h1,h2,fh1,fh2,u,b;
b=0; u=1.05;
h1=b+GOLDEN2*(u-b);
h2=b+GOLDEN1*(u-b);
xh(h1); fh1=fx();
xh(h2); fh2=fx();
while( fabs(fh2-fh1)>ZERO ) {
if( fh1>fh2 ) {
b=h1;
h1=h2; fh1=fh2;
h2=b+GOLDEN1*(u-b);
xh(h2); fh2=fx();
}
else {
u=h2;
h2=h1; fh2=fh1;
h1=b+GOLDEN2*(u-b);
xh(h1); fh1=fx();
}
}
if( fh1>fh2 )
hk=h2;
else
hk=h1;
}
void UnconstrainedOptim::BSsearch(double &h, double &g0, double &f_1) {
// f0=f(xk), f_1=f(x(k-1))
double f0=f;
double m,gh;
double h0,fh0,gh0;
double a,b,fa,fb,ga,gb;
int gbknown;
int termin;
if( fabs(g0)<1E-10 ) {
f_1=f0;
return;
}
m=(flb-f0)/(RHO*g0);
h0=0; fh0=f0; gh0=g0;
// f_1:=1E30; // Makes h become 1
h=max(f_1-f0, 10*ZERO);
h=min(-2*h/g0, 1);
termin=0;
// ===== Bracketing Phase ============================
while( termin==0 ) {
xh(h); fx();
if( f<=flb )
termin=2;
else if( f>f0+h*RHO*g0 || f>=fh0 ) {
a=h0; fa=fh0; ga=gh0;
b=h; fb=f; gbknown=0;
termin=1;
}
else {
gx();
gh=gs();
if( fabs(gh)<=-SIGMA*g0 )
termin=2;
else if( gh>0 ) {
a=h; fa=f; ga=gh;
b=h0; fb=fh0; gb=gh0;
gbknown=1;
termin=1;
}
else {
a=2*h-h0;
if( m<=a ) {
h0=h; fh0=f; gh0=gh;
h=m;
}
else {
b=min(h+TAW1*(h-h0), m);
a=UnconstrainedOptim::chooseC(a,b,h0,fh0,gh0,h,f,gh);
h0=h; fh0=f; gh0=gh;
h=a;
}
}
}
}
if( termin==1 ) {
termin=0;
// ===== Sectioning Phase ============================
while( termin==0 ) {
if( gbknown )
h=UnconstrainedOptim::chooseC(a+TAW2*(b-a), b-TAW3*(b-a), a, fa, ga, b, fb, gb);
else
h=UnconstrainedOptim::chooseQ(a+TAW2*(b-a), b-TAW3*(b-a), a, fa, ga, b, fb);
xh(h); fx();
// **** Line Search End Condition by Fletcher ****
if( (a-h)*ga<=ZERO )
termin=1;
else {
if( f>f0+h*RHO*g0 || f>=fa ) {
b=h; fb=f; gbknown=0;
}
else {
gx();
gh=gs();
if( fabs(gh)<=-SIGMA*g0 )
termin=1;
else {
if( (b-a)*gh>=0 ) {
b=a; fb=fa; gb=ga; gbknown=1;
}
a=h; fa=f; ga=gh;
}
}
}
}
}
f_1=f0; g0=gh;
}
void UnconstrainedOptim::Ssearch(double &h, double &g0, double &f_1) {
// f0=f(xk), f_1=f(x(k-1))
double f0=f;
double m,a,b,a1;
double gh,fa,ga;
int cont;
if( fabs(g0)<1E-10 ) {
f_1=f0;
return;
}
m=(flb-f0)/(RHO*g0);
a=0; fa=f0; ga=g0;
// f_1:=1E30; // Makes h become 1
h=max(f_1-f0, 10*ZERO);
h=min(-2*h/g0, 1);
b=2*m+1;
cont=1;
while( cont ) {
xh(h); fx();
if( f<=flb )
cont=0;
else if( f>f0+h*RHO*g0 || f>=fa ) {
if( fabs((h-a)*ga)<=ZERO )
cont=0;
else {
b=h;
h=UnconstrainedOptim::chooseQ(a+T1*(h-a),h-T2*(h-a),a,fa,ga,h,f);
}
}
else {
gx();
gh=gs();
if( fabs(gh)<=-SIGMA*g0 )
cont=0;
else {
a1=h;
if( (b-a)*gh<0 ) {
if( b<=m )
h=UnconstrainedOptim::chooseC(h+T5*(b-h),b-T6*(b-h), a,fa,ga, h,f,gh);
else
h=UnconstrainedOptim::chooseC(min(T3*h,m),min(T4*h,m),a,fa,ga,h,f,gh);
}
else {
h=UnconstrainedOptim::chooseC(a+T1*(h-a),h-T2*(h-a),a,fa,ga,h,f,gh);
b=a;
}
a=a1; fa=f; ga=gh;
}
}
}
f_1=f0; g0=gh;
}
void UnconstrainedOptim::RSsearch(double &h, double &g0, double &f_1) {
// f0=f(xk), f_1=f(x(k-1))
double f0=f;
double m,a,b,a1,temp;
double gh,fa,ga,fa1,fb;
int cont;
if( fabs(g0)<1E-10 ) {
f_1=f0;
return;
}
m=(flb-f0)/(RHO*g0);
a=0; fa=f0; ga=g0;
// f_1:=1E30; // Makes h become 1
h=max(f_1-f0, 10*ZERO);
h=min(-2*h/g0, 1);
b=2*m+1;
cont=1;
while( cont==1 ) {
a1=a; fa1=fa;
cont=2;
xh(h); fx();
if( f<=flb )
cont=0;
else if( f>f0+h*RHO*g0 || f>=fa ) {
if( fabs((h-a)*ga)<=ZERO )
cont=0;
else {
temp=h;
if( b<=m )
h=UnconstrainedOptim::chooseCC(a+T1*(h-a),h-T2*(h-a),a,fa,ga,h,f,b,fb);
else
h=UnconstrainedOptim::chooseQ(a+T1*(h-a),h-T2*(h-a),a,fa,ga,h,f);
b=temp; fb=f;
cont=1;
}
}
while( cont==2 && fa1-f>max(10*ZERO,-SIGMA*g0*fabs(a1-h)) ) {
temp=h;
if( b<=m )
h=UnconstrainedOptim::chooseCC(a+T1*(b-a),b-T2*(b-a),a,fa,ga,a1,fa1,h,f);
else
h=UnconstrainedOptim::chooseCC(a+T1*(h-a),h-T2*(h-a),a,fa,ga,a1,fa1,h,f);
a1=temp; fa1=f;
xh(h); fx();
if( f<=flb )
cont=0;
else if( f>f0+h*RHO*g0 || f>=fa ) {
if( fabs((h-a)*ga)<=ZERO )
cont=0;
else {
temp=h;
h=UnconstrainedOptim::chooseCC(a+T1*(h-a),h-T2*(h-a),a,fa,ga,a1,fa1,h,f);
b=temp; fb=f;
cont=1;
}
}
}
if( cont==2 ) {
cont=1;
if( fa1<f ) {
h=a1; f=fa1;
xh(h);
}
gx();
gh=gs();
if( fabs(gh)<=-SIGMA*g0 )
cont=0;
else {
temp=h;
if((b-a)*gh<0) {
if( b<=m )
h=UnconstrainedOptim::chooseC(h+T5*(b-h),b-T6*(b-h), a,fa,ga,h,f,gh);
else
h=UnconstrainedOptim::chooseC(min(T3*h,m),min(T4*h,m),a,fa,ga,h,f,gh);
}
else {
h=UnconstrainedOptim::chooseC(a+T1*(h-a),h-T2*(h-a),a,fa,ga,h,f,gh);
b=a; fb=fa;
}
a=temp; fa=f; ga=gh;
}
}
}
f_1=f0; g0=gh;
}
int UnconstrainedOptim::DFPoptimize(int lineSearchMethod) {
/// DFP //
cnls=0; cnf=0; cng=0;
double *gm=new double[n];
double *hg=new double[n];
int i,j;
double **h;
h=new double *[n];
i=n;
while(--i>=0) {
h[i]=new double[n];
j=n;
while(--j>=0)
h[i][j]=0;
h[i][i]=1;
}
double fk_1=1E30;
double g0;
fx();
gx();
int cont=1;
do {
i=n;
double w;
while(--i>=0) {
w=0;
j=n;
while(--j>=0)
w-=h[i][j]*g[j];
sk[i]=w;
}
g0=gs();
update();
switch(lineSearchMethod) {
case GOLDENSEARCH:
fk_1=f;
goldenSearch(w);
break;
case BSSEARCH:
BSsearch(w,g0,fk_1);
break;
case SSEARCH:
Ssearch(w,g0,fk_1);
break;
case RSSEARCH:
RSsearch(w,g0,fk_1);
}
++cnls;
if( fk_1-f<ZERO )
cont=0;
if(!gknown)
gx();
w=0;
i=n;
while(--i>=0) {
gm[i]=g[i]-gk[i];
w+=(x[i]-xk[i])*gm[i];
}
i=n;
double v;
while(--i>=0) {
v=0;
j=n;
while(--j>=0)
v+=h[i][j]*gm[j];
hg[i]=v;
}
v=0;
i=n;
while(--i>=0)
v+=gm[i]*hg[i];
if( fabs(w)<ZERO || fabs(v)<ZERO )
cont=0;
else {
i=n;
while(--i>=0) {
j=n;
while(--j>=0)
h[i][j]+=(x[i]-xk[i])*(x[j]-xk[j])/w-hg[i]*hg[j]/v;
}
}
}
while(cont);
delete hg;
delete gm;
i=n;
while(--i>=0)
delete h[i];
delete h;
if( fabs(f-fopt)<ZERO )
return 0;
else
return floor((1+(log(fabs(f-fopt))-log(ZERO))/log(10.0))/2.0);
}
int UnconstrainedOptim::BFGSoptimize(int lineSearchMethod) {
/// BFGS //
cnls=0; cnf=0; cng=0;
double *gm=new double[n];
double *hg=new double[n];
int i,j;
double **h;
h=new double *[n];
i=n;
while(--i>=0) {
h[i]=new double[n];
j=n;
while(--j>=0)
h[i][j]=0;
h[i][i]=1;
}
double fk_1=1E30;
double g0;
fx();
gx();
int cont=1;
do {
i=n;
double w;
while(--i>=0) {
w=0;
j=n;
while(--j>=0)
w-=h[i][j]*g[j];
sk[i]=w;
}
g0=gs();
update();
switch(lineSearchMethod) {
case GOLDENSEARCH:
fk_1=f;
goldenSearch(w);
break;
case BSSEARCH:
BSsearch(w,g0,fk_1);
break;
case SSEARCH:
Ssearch(w,g0,fk_1);
break;
case RSSEARCH:
RSsearch(w,g0,fk_1);
}
++cnls;
if( fk_1-f<ZERO )
cont=0;
if(!gknown)
gx();
w=0;
i=n;
while(--i>=0) {
gm[i]=g[i]-gk[i];
w+=(x[i]-xk[i])*gm[i];
}
i=n;
double v;
while(--i>=0) {
v=0;
j=n;
while(--j>=0)
v+=h[i][j]*gm[j];
hg[i]=v;
}
v=0;
i=n;
while(--i>=0)
v+=gm[i]*hg[i];
if( fabs(w)<ZERO )
cont=0;
else {
v=1+v/w;
i=n;
while(--i>=0) {
j=n;
while(--j>=0) {
double di=x[i]-xk[i];
double dj=x[j]-xk[j];
h[i][j]+=v*di*dj/w-(di*hg[j]+dj*hg[i])/w;
}
}
}
}
while(cont);
delete hg;
delete gm;
i=n;
while(--i>=0)
delete h[i];
delete h;
if( fabs(f-fopt)<ZERO )
return 0;
else
return floor((1+(log(fabs(f-fopt))-log(ZERO))/log(10.0))/2.0);
}
int UnconstrainedOptim::FRoptimize(int lineSearchMethod, int reset) {
/// Fletcher_Reeves ///
if(reset)
return FRoptimizeWithReset(lineSearchMethod);
cnls=0; cnf=0; cng=0;
double w,v,g0;
int i;
double fk_1=1E30;
fx();
gx();
v=0;
i=n;
while(--i>=0) {
sk[i]=0;
v+=g[i]*g[i];
}
w=0;
int cont=1;
do {
i=n;
while(--i>=0)
sk[i]=w*sk[i]-g[i];
g0=gs();
if(g0>0) {
i=n;
while(--i>=0)
sk[i]=-g[i];
g0=-v;
}
update();
switch(lineSearchMethod) {
case GOLDENSEARCH:
fk_1=f;
goldenSearch(w);
break;
case BSSEARCH:
BSsearch(w,g0,fk_1);
break;
case SSEARCH:
Ssearch(w,g0,fk_1);
break;
case RSSEARCH:
RSsearch(w,g0,fk_1);
}
++cnls;
if( fk_1-f<ZERO )
cont=0;
if(!gknown)
gx();
w=0;
i=n;
while(--i>=0)
w+=g[i]*g[i];
if( fabs(v)<ZERO )
cont=0;
else {
double z=w;
w=w/v;
v=z;
}
}
while(cont);
if( fabs(f-fopt)<ZERO )
return 0;
else
return floor((1+(log(fabs(f-fopt))-log(ZERO))/log(10.0))/2.0);
}
int UnconstrainedOptim::PRoptimize(int lineSearchMethod, int reset) {
/// Polak_Ribiere ///
if(reset)
return PRoptimizeWithReset(lineSearchMethod);
cnls=0; cnf=0; cng=0;
double w,v,g0;
int i;
double fk_1=1E30;
fx();
gx();
v=0;
i=n;
while(--i>=0) {
sk[i]=0;
v+=g[i]*g[i];
}
w=0;
int cont=1;
do {
i=n;
while(--i>=0)
sk[i]=w*sk[i]-g[i];
g0=gs();
if(g0>0) {
i=n;
while(--i>=0)
sk[i]=-g[i];
g0=-v;
}
update();
switch(lineSearchMethod) {
case GOLDENSEARCH:
fk_1=f;
goldenSearch(w);
break;
case BSSEARCH:
BSsearch(w,g0,fk_1);
break;
case SSEARCH:
Ssearch(w,g0,fk_1);
break;
case RSSEARCH:
RSsearch(w,g0,fk_1);
}
++cnls;
if( fk_1-f<ZERO )
cont=0;
if(!gknown)
gx();
w=0; g0=0;
i=n;
while(--i>=0) {
w+=g[i]*g[i];
g0+=g[i]*gk[i];
}
if( fabs(v)<ZERO )
cont=0;
else {
double z=w;
w=(w-g0)/v;
v=z;
}
}
while(cont);
if( fabs(f-fopt)<ZERO )
return 0;
else
return floor((1+(log(fabs(f-fopt))-log(ZERO))/log(10.0))/2.0);
}
int UnconstrainedOptim::FRoptimizeWithReset(int lineSearchMethod) {
/// Fletcher_Reeves ///
cnls=0; cnf=0; cng=0;
double w,v,g0;
int i;
double fk_1=1E30;
fx();
gx();
v=0;
i=n;
while(--i>=0)
v+=g[i]*g[i];
w=0;
int cont=1;
int r=0;
do {
if(r) {
i=n;
while(--i>=0)
sk[i]=w*sk[i]-g[i];
}
else {
i=n;
while(--i>=0)
sk[i]=-g[i];
g0=-v;
}
++r;
if(r==n)
r=0;
g0=gs();
update();
switch(lineSearchMethod) {
case GOLDENSEARCH:
fk_1=f;
goldenSearch(w);
break;
case BSSEARCH:
BSsearch(w,g0,fk_1);
break;
case SSEARCH:
Ssearch(w,g0,fk_1);
break;
case RSSEARCH:
RSsearch(w,g0,fk_1);
}
++cnls;
if( fk_1-f<ZERO )
cont=0;
if(!gknown)
gx();
w=0;
i=n;
while(--i>=0)
w+=g[i]*g[i];
if( fabs(v)<ZERO )
cont=0;
else {
double z=w;
w=w/v;
v=z;
}
}
while(cont);
if( fabs(f-fopt)<ZERO )
return 0;
else
return floor((1+(log(fabs(f-fopt))-log(ZERO))/log(10.0))/2.0);
}
int UnconstrainedOptim::PRoptimizeWithReset(int lineSearchMethod) {
/// Polak_Ribiere ///
cnls=0; cnf=0; cng=0;
double w,v,g0;
int i;
double fk_1=1E30;
fx();
gx();
v=0;
i=n;
while(--i>=0)
v+=g[i]*g[i];
w=0;
int cont=1;
int r=0;
do {
if(r) {
i=n;
while(--i>=0)
sk[i]=w*sk[i]-g[i];
}
else {
i=n;
while(--i>=0)
sk[i]=-g[i];
g0=-v;
}
++r;
if(r==n)
r=0;
g0=gs();
update();
switch(lineSearchMethod) {
case GOLDENSEARCH:
fk_1=f;
goldenSearch(w);
break;
case BSSEARCH:
BSsearch(w,g0,fk_1);
break;
case SSEARCH:
Ssearch(w,g0,fk_1);
break;
case RSSEARCH:
RSsearch(w,g0,fk_1);
}
++cnls;
if( fk_1-f<ZERO )
cont=0;
if(!gknown)
gx();
w=0; g0=0;
i=n;
while(--i>=0) {
w+=g[i]*g[i];
g0+=g[i]*gk[i];
}
if( fabs(v)<ZERO )
cont=0;
else {
double z=w;
w=(w-g0)/v;
v=z;
}
}
while(cont);
if( fabs(f-fopt)<ZERO )
return 0;
else
return floor((1+(log(fabs(f-fopt))-log(ZERO))/log(10.0))/2.0);
}
!!!!!!!!!!!!!!!!!!
sample!!!!
// =============================================================================
// a Sample application Case of unconstrained optimizatoin library.
// programmer: A. Katanforoush,
// last update: 7/2/2003,
// =============================================================================
//
// Write a suit constructor then override F() and G() in derived class.
//
// I have provided two famous test problem for unconstrained optimization as sample cases.
// Rosenbrock: F(x) = 100 ( x^2 - y )^2 + ( 1 - x )^2
// Powell: F(x) = (x+10y)^2 + 5(z-w)^2 + (y-2*z)^4 + 10(x-w)^4
// =============================================================================
#ifndef SAMPLECASE_H
#define SAMPLECASE_H
#include "uncoptim.h"
class Rosenbrock: public UnconstrainedOptim {
public:
Rosenbrock(int lsmethod):UnconstrainedOptim(2,0,0,lsmethod) {
x[0]=-1.2;
x[1]=1;
}
void F() {
double w1=x[1]-x[0]*x[0];
double w2=1-x[0];
f=100*w1*w1+w2*w2;
}
void G() {
g[1]=200*(x[1]-x[0]*x[0]);
g[0]=-2*(x[0]*g[1]+1-x[0]);
}
};
class Powell: public UnconstrainedOptim {
public:
Powell(int lsmethod):UnconstrainedOptim(4,0,0,lsmethod) {
x[0]=3;
x[1]=-1;
x[2]=0;
x[3]=1;
}
void F() {
double w1=x[0]+10*x[1];
double w2=x[2]-x[3];
double w3=x[1]-2*x[2]; w3*=w3;
double w4=x[0]-x[3]; w4*=w4;
f=w1*w1+5*w2*w2+w3*w3+10*w4*w4;
}
void G() {
double x5=x[0]+10*x[1];
double x6=x[2]-x[3];
double x7=x[0]-x[3]; x7=x7*x7*x7;
double x8=x[1]-2*x[2]; x8=x8*x8*x8;
g[0]=2*x5+40*x7;
g[1]=20*x5+4*x8;
g[2]=10*x6-8*x8;
g[3]=-10*x6-40*x7;
}
};
#endif
// =============================================================================
// a Sample application Case of unconstrained optimizatoin library.
// programmer: A. Katanforoush,
// last update: 7/2/2003,
// =============================================================================
//
// Comparison different optimization methods by the test problem.
//
// =============================================================================
#include <stdio.h>
#include "SampleCase.h"
int main() {
//Rosenbrock r(3);
Powell r(3);
int opt;
opt=r.DFPoptimize();
printf("DFP: x1 = %-0.10lg x2 = %-0.10lg F = %-0.10lg/n", r.x[0], r.x[1], r.f);
printf(" cnls = %i cnf = %i cng = %i optimality = %i/n/n", r.cnls, r.cnf, r.cng, opt);
r.x[0]=3;
r.x[1]=-1;
r.x[2]=0;
r.x[3]=1;
opt=r.BFGSoptimize();
printf("BFGS: x1 = %-0.10lg x2 = %-0.10lg F = %-0.10lg/n", r.x[0], r.x[1], r.f);
printf(" cnls = %i cnf = %i cng = %i optimality = %i/n/n", r.cnls, r.cnf, r.cng, opt);
r.x[0]=3;
r.x[1]=-1;
r.x[2]=0;
r.x[3]=1;
opt=r.FRoptimize(0);
printf("FR: x1 = %-0.10lg x2 = %-0.10lg F = %-0.10lg/n", r.x[0], r.x[1], r.f);
printf(" cnls = %i cnf = %i cng = %i optimality = %i/n/n", r.cnls, r.cnf, r.cng, opt);
r.x[0]=3;
r.x[1]=-1;
r.x[2]=0;
r.x[3]=1;
opt=r.FRoptimize(1);
printf("FR: x1 = %-0.10lg x2 = %-0.10lg F = %-0.10lg/n", r.x[0], r.x[1], r.f);
printf(" cnls = %i cnf = %i cng = %i optimality = %i/n/n", r.cnls, r.cnf, r.cng, opt);
r.x[0]=3;
r.x[1]=-1;
r.x[2]=0;
r.x[3]=1;
opt=r.PRoptimize(0);
printf("PR: x1 = %-0.10lg x2 = %-0.10lg F = %-0.10lg/n", r.x[0], r.x[1], r.f);
printf(" cnls = %i cnf = %i cng = %i optimality = %i/n/n", r.cnls, r.cnf, r.cng, opt);
r.x[0]=3;
r.x[1]=-1;
r.x[2]=0;
r.x[3]=1;
opt=r.PRoptimize(1);
printf("PR: x1 = %-0.10lg x2 = %-0.10lg F = %-0.10lg/n", r.x[0], r.x[1], r.f);
printf(" cnls = %i cnf = %i cng = %i optimality = %i/n/n", r.cnls, r.cnf, r.cng, opt);
getchar();
return 0;
}
8882
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
