1.txt.bak

这个是网上的最全的关于共轭梯度法的c++程序。

最优化-无约束共轭梯度法程序（c++）2008/01/03 12:49/////////////////////////////////////////
/////          vector.h头文件         /////
/////       定义向量及其基本运算      /////
/////////////////////////////////////////

#include<math.h>
#define MAXLENGTH 10

//向量定义
typedef struct{             
     int tag;                 //行列向量标志。行向量为0，列向量为1。
     int dimension;           //向量的维数
     double elem[MAXLENGTH];     //向量的元素
}vector;

vector vecCreat(int tag, int n){
//建立维数为n的向量
     vector x;
     x.tag=tag;
     x.dimension=n;
     for(int i=0;i<n;++i){
         cout<<"请输入第"<<i+1<<"分量:";
         cin>>x.elem[i];
     }
     return x;
}

double vecMultiply(vector a, vector b){
//行向量与列向量乘法运算
     if((a.tag!=b.tag)&&(a.dimension==b.dimension)){//相乘条件
         double c=0;
         for(int i=0;i<a.dimension;++i)
             c+=a.elem[i]*b.elem[i];
         return c;
     }
}

vector vecAdd(vector a, vector b){
//向量加法运算
     if((a.tag==b.tag)&&(a.dimension==b.dimension)){//相加条件
         vector c;
         c.dimension=a.dimension;
         c.tag=a.tag;
         for(int i=0;i<c.dimension;++i)
             c.elem[i]=a.elem[i]+b.elem[i];
         return c;
     }
}

vector vecConvert(vector a){
//向量转置
     if(a.tag==0)a.tag=1;
     if(a.tag==1)a.tag=0;
     return a;
}

double vecMole(vector a){
//求向量模
     double sum=0;
     for(int i=0;i<a.dimension;++i)
         sum+=a.elem[i]*a.elem[i];
     sum=sqrt(sum);    
     return sum;
}

vector numMultiply(double n,vector a){
//数乘向量
     for(int i=0;i<a.dimension;++i)
         a.elem[i]=n*a.elem[i];
     return a;
}

void showPoint(vector x, double f){
//显示点，解.
     cout<<"--- x=(";
     for(int i=0;i<x.dimension;++i){
         cout<<x.elem[i];
         if(i!=x.dimension-1)cout<<",";
     }
     cout<<") --- f(x)="<<f<<endl;
     cout<<endl;    
}
/////////////////////////////////////////////////////////
////////////// function.h头文件 /////////////////////////
//////// 函数及其导数的设置均在此文件完成////////////////
/////////////////////////////////////////////////////////
// * 无 约 束 问 题 *                                   //
// 目标函数--在vecFun(vector vx)中修改,x1改成x[1]       //
//               x2改成x[2],依此类推...                  //
/////////////////////////////////////////////////////////

#include <iostream.h>
#include "vector.h"
#define SIZE 10
#define MAX 1e300

double min(double a, double b){
//求两个数最小
     return a<b?a:b;
}

double max(double a, double b){
//求两个数最大
     return a>b?a:b;
}

vector vecGrad(double (*pf)(vector x),vector pointx){
//求解向量的梯度
//采用理查德外推算法计算函数pf在pointx点的梯度grad;
     vector grad;
     grad.tag=1; grad.dimension=pointx.dimension;//初始化偏导向量
     vector tempPnt1,tempPnt2;            //临时向量
     double h=1e-3;
     for(int i=0;i<pointx.dimension;++i){
         tempPnt1=tempPnt2=pointx;
         tempPnt1.elem[i]+=0.5*h;
         tempPnt2.elem[i]-=0.5*h;
         grad.elem[i]=(4*(pf(tempPnt1)-pf(tempPnt2)))/(3*h);
         tempPnt1.elem[i]+=0.5*h;
         tempPnt2.elem[i]-=0.5*h;
         grad.elem[i]-=(pf(tempPnt1)-pf(tempPnt2))/(6*h);
     }
     return grad;
}

double vecFun(vector vx){
//最优目标多元函数
     double x[SIZE];
     for(int i=0;i<vx.dimension;++i)
         x[i+1]=vx.elem[i];
     //----------约束目标函数--------------
     //return x[1]*x[1]+x[2]*x[2];
     //----------无约束正定函数--------------
     //return x[1]*x[1]+4*x[2]*x[2]+9*x[3]*x[3]-2*x[1]+18*x[3];//例一
     //----------无约束非二次函数--------------
     //return (1-x[1])*(1-x[1])+(1-x[4])*(1-x[4])+(x[1]*x[1]-x[2])*(x[1]*x[1]-x
[2])+(x[2]*x[2]-x[3])*(x[2]*x[2]-x[3])+(x[3]*x[3]-x[4])*(x[3]*x[3]-x[4]);//例二

}

double vecFun_Si(vector vx){
//不等式约束函数
     double x[SIZE];
     for(int i=0;i<vx.dimension;++i)
         x[i+1]=vx.elem[i];
     return x[1]+1;//不等式约束函数
}

double vecFun_Hi(vector vx){
//等式约束函数
     double x[SIZE];
     for(int i=0;i<vx.dimension;++i)
         x[i+1]=vx.elem[i];
     return x[1]+x[2]-2;//等式约束函数
}

double vecFun_S(vector x,double v,double l,double u){
//约束问题转化为无约束问题的增广目标函数F(x,v,l,u)
     double sum1=0,sum2=0,sum3=0;//分别定义三项的和
     //sum1
     double temp=max(0,v-2*u*vecFun_Si(x));
     sum1=(temp*temp-v*v)/(4*u);
     //sum2
     sum2=l*vecFun_Hi(x);
     //sum3
     sum3=u*vecFun_Hi(x)*vecFun_Hi(x);    
     //F(x,v,l,u)=f(x)+sum1-sum2+sum3
     return vecFun(x)+sum1-sum2+sum3;
}

vector vecFunD_S(vector x,double v,double l,double u){//利用重载函数实现目标函
数F(x,v,l,u)的导数
//约束问题转化为无约束问题的增广目标函数F(x,v,l,u)的导函数
     vector sum1,sum2,sum3;//分别定义三项导数的和
     //sum1
     sum1.dimension=x.dimension;     sum1.tag=1;
     sum2.dimension=x.dimension;     sum2.tag=1;
     sum3.dimension=x.dimension;     sum3.tag=1;
     double temp=max(0,v-2*u*vecFun_Si(x));
     if(temp==0){
         for(int i=0;i<sum1.dimension;++i)
             sum1.elem[i]=0;
     }
     else{
         sum1=numMultiply(-(v-2*u*vecFun_Si(x)),vecGrad(vecFun_Si,x));
     }
     //-sum2
     sum2=numMultiply(-l,vecGrad(vecFun_Hi,x));
     //sum3
     sum3=numMultiply(2*u*vecFun_Hi(x),vecGrad(vecFun_Hi,x));    
     //F=f(x)+sum1-sum2+sum3
     return vecAdd(vecAdd(vecGrad(vecFun,x),sum1),vecAdd(sum2,sum3));
}
///////////////////////////////////////////////////
/////             nonrestrict.h头文件           /////
/////   包含无约束问题的求解函数: 直线搜索      /////
/////          共轭梯度法,H终止准则             /////
///////////////////////////////////////////////////
#include "restrict.h"

vector lineSearch(vector x, vector p ,double t){
//从点x沿直线p方向对目标函数f(x)作直线搜索得到极小点x2,t为初始步长。
     vector x2;
     double l=0.1,n=0.4;          //条件1和2的参数  
     double a=0,b=MAX;            //区间 
     double f1=0,f2=0,g1=0,g2=0;
     int i=0;                     //迭代次数
     do{
         if(f2-f1>l*t*g1){b=t;t=(a+b)/2;}          //改变b,t循环
         do{
             if(g2<n*g1){a=t;t=min((a+b)/2,2*a);} //改变a,t循环
             f1=vecFun(x);                                    //f1=f(x)
             g1=vecMultiply(vecConvert(vecGrad(vecFun,x)),p);//g1=∨f(x)*p
             x2=vecAdd(numMultiply(t,p),x);                   //x2=x+t*p
             if(i++ && vecFun(x2)==f2){return x2;}            //【直线搜索进入无
限跌代，则此次跌代结束。返回当前最优点】
             f2=vecFun(x2);                                   //f2=f(x2)
             g2=vecMultiply(vecConvert(vecGrad(vecFun,x2)),p);       //g2=∨f(x2
)*p    
             //cout<<"("<<x2.elem[0]<<","<<x2.elem[1]<<","<<x2.elem[2]<<","<<x2
.elem[3]<<");"; //输出本次结果
             //cout<<"t="<<t<<";";
             //cout<<"f(x)="<<f2<<endl;
         }
         while(g2<n*g1);       //不满足条件ii,则改变a,t循环
     }
     while(f2-f1>l*t*g1);      //不满足条件i,则改变b,t循环
     return x2;
}

int Himmulblau(vector x0, vector x1, double f0, double f1){
//H终止准则.给定Xk,Xk+1,Fk,Fk+1,判断Xk+1是否是极小点.返回1是极小,否则返回0
     double c1,c2,c3;//定义并初始化终止限
     c1=c2=10e-5;
     c3=10e-4;
     if(vecMole(vecGrad(vecFun,x1))<c3)
         if(vecMole(vecAdd(x1,numMultiply(-1,x0)))/(vecMole(x0)+1)<c1)
             if(fabs(f1-f0)/(fabs(f0)+1)<c2)
                 return 1;
     return 0;
}

void Fletcher_Reeves(vector xx,vector &minx, double &minf){
//Fletcher-Reeves共轭梯度法.
//给定初始点xx.对vecFun函数求得最优点x和最优解f,分别用minx和minf返回
     int k=0,j=0;        //迭代次数,j为总迭代次数
     double c=10e-1;     //终止限c
     int n=xx.dimension;//问题的维数
     double f0,f,a;      //函数值f(x),a为使p方向向量共轭的系数
     vector g0,g;        //梯度g0,g
     vector p0,p;        //搜索方向P0,p
     vector x,x0;        //最优点和初始点
     double t=1;         //直线搜索的初始步长=1
     x=xx;                          //x0
     f=vecFun(x);                   //f0=f(x0)
     g=vecGrad(vecFun,x);           //g0
     p0=numMultiply(-1,g);          //p0=-g0,初始搜索方向为负梯度方向
     do{
         x0=x;f0=f;g0=g;              
         x=lineSearch(x0,p0,t);        //从点x出发,沿p方向作直线搜索
         f=vecFun(x);                  //f=f(x)
         g=vecGrad(vecFun,x);                 //g=g(x)
         if(Himmulblau(x0,x,f0,f)==1){//满足H终止准则,返回最优点及最优解。
             cout<<"* 总迭代次数："<<j<<" ----"<<endl;
             minx=x; minf=f;
             break;
         }
         else{                         //不满足H准则
             cout<<"* 第"<<j<<"次迭代";
             showPoint(x,f);           //显示中间结果x,f    

             if(k==n){                 //若进行了n+1次迭代
                 k=0;                  //重置k=0,p0=-g
                 p0=numMultiply(-1,g);
                 t=vecMole(vecAdd(x,numMultiply(-1,x0)))/vecMole(p0);//下次直线
搜索步长
                 continue;             //进入下一次循环,重新直线搜索
             }                        
             else{                      //否则,计算共轭向量p
                 a=vecMole(g)*vecMole(g)/(vecMole(g0)*vecMole(g0));
                 p=vecAdd(numMultiply(-1,g),numMultiply(a,p0));
             }
         }
         if(fabs(vecMultiply(vecConvert(p),g))<=c){//共轭梯度方向下降或上升量很
小
             p0=numMultiply(-1,g);//p0=-g,k=0
             k=0;
         }
         else if(vecMultiply(vecConvert(p),g)<-c){//共轭梯度方向下降并且下降量保
证
             p0=p;                 //p0=p,k=k+1
             ++k;
         }
         else{                                     //共轭梯度方向上升并且下降量保
证
             p0=numMultiply(-1,p);//p0=-p,k=k+1
             ++k;
         }
         t=vecMole(vecAdd(x,numMultiply(-1,x0)))/vecMole(p0);//下次直线搜索步长

     }
     while(++j);
}

///////////////////////////////////////////////////
/////              main.h头文件                 /////
/////   主程序文件，控制整个程序的流程          /////
///////////////////////////////////////////////////
#include "nonrestrict.h"

void printSelect(){
     cout<<"************** 最优化方法 ***************"<<endl;
     cout<<"*****************************************"<<endl;
     cout<<"***   请选择：                          ***"<<endl;
     cout<<"***   1. 无约束最小化问题(共轭梯度法)   ***"<<endl;
     cout<<"***   2. 约束最小化问题(H乘子法)        ***"<<endl;
     cout<<"***   3. 退出(exit)                     ***"<<endl;
     cout<<"*****************************************"<<endl;
}

void sub1(){
     char key;
     cout<<"--------------------共轭梯度法求无约束最小化问题----------------"<<
endl;
     cout<<"步骤:"<<endl;
     cout<<"◆ 1. 输入f(x)及其导数.(参考function.h文件说明)"<<endl;
     cout<<"-----确认是否已输入<目标函数>及其<导数>? (Y/N)";
     cin>>key;
     if(key=='N' || key=='n')return; 
     else if(key=='Y' || key=='y'){
         vector x0; //起始点
         int m;
         cout<<"◆ 2. 起始点X0初始化"<<endl<<"-----请输入X0的维数:";
         cin>>m;
         x0=vecCreat(1,m);
         vector minx;//求得的极小点
         double minf;//求得的极小值
         Fletcher_Reeves(x0,minx,minf);
         cout<<"◆ 最优解及最优值:";
         showPoint(minx,minf);
     }
}

void sub2(){
     char key;
     cout<<"------------------ H乘子法求约束最小化问题 -----------------"<<endl
;
     cout<<"步骤:"<<endl;
     cout<<"◆ 1. 输入f(x)及其导数.(参考function.h文件说明)"<<endl;
     cout<<"-----确认是否已输入<目标函数及导数>和<约束条件>? (Y/N)";
     cin>>key;
     if(key=='N' || key=='n')return; 
     else if(key=='Y' || key=='y'){
         vector x0; //起始点
         int m;
         cout<<"◆ 2. 起始点X0初始化"<<endl<<"-----请输入X0的维数:";
         cin>>m;
         x0=vecCreat(1,m);
         vector minx;//求得的极小点
         double minf;//求得的极小值
         Hesternes(x0,minx,minf);
         showPoint(minx,minf);
     }
}

void main(){
     int mark=1;
     while(mark){
         printSelect();
         int sel;
         cin>>sel;
         switch(sel){
             case 1:
                 sub1();
                 break;
             case 2:
                 sub2();
                 break;
             case 3:
                 mark=0;
                 break;
         };
     }