`

最优化-无约束共轭梯度法程序(c++)

阅读更多

copyright @ http://blog.chinaunix.net/uid-253851-id-2140409.html

 

/////////////////////////////////////////

/////         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;
        };
    }
}

分享到:
评论

相关推荐

    最优化-约束/无约束共轭梯度法程序(c++)

    最优化-约束/ 无约束共轭梯度法程序(c++)

    共轭梯度pdf阻尼牛顿算法c++实现

    在优化领域,共轭梯度法(Conjugate Gradient Method)和阻尼牛顿法(Damped Newton Method)是两种广泛使用的数值优化方法,尤其在解决线性系统和非线性最小化问题时。这两种方法在工程计算和科学计算中扮演着重要...

    最优化程序

    在给定的压缩包文件中,包含了多种最优化算法的实现,包括最速下降法(精确步长搜索和Wolfe条件非精确步长搜索)以及共轭梯度法和预处理共轭梯度法,这些都是解决二次线性方程组优化问题的常用技术。 1. **最速下降...

    最优化作业C++源代码.rar

    本压缩包“最优化作业C++源代码.rar”提供了几种常用的最优化方法的C++实现,让我们逐一探讨这些算法及其在实际中的应用。 首先,牛顿法是一种迭代优化算法,主要用于寻找函数的局部极小值。它通过构建函数的泰勒...

    《最优化理论与算法》-陈宝林_深度学习_C++_

    3. **二阶优化方法**:如牛顿法、拟牛顿法(如BFGS和L-BFGS)以及共轭梯度法。 4. **适应性学习率算法**:如Adagrad、RMSprop和Adam,它们能自适应地调整每个参数的学习率。 5. **约束优化**:处理约束条件下的...

    最优化计算机原理与算法程序设计

    共轭梯度法适用于大规模线性方程组的求解,特别是解决多维二次函数的优化问题。它利用共轭方向的性质,避免了矩阵求逆,特别适合于大规模的稀疏系统。 牛顿型方法是一种利用函数二阶导数信息的迭代法,特别适合于...

    powell鲍威尔法(c++).7z

    在实际应用中,鲍威尔法可能不如其他更现代的优化算法(如梯度下降、共轭梯度法或遗传算法)高效,尤其是在高维问题上。但它的优点在于简单且不需要目标函数的导数,这使得它在某些情况下具有优势,特别是当目标函数...

    Wolfe Powell 最优化方法 C++ 程序

    - **改进搜索方向**:探索不同的搜索方向,例如共轭梯度法等,以提高搜索效率。 - **多线程并行化**:利用现代计算机的多核处理器优势,加快计算速度。 综上所述,给定的 C++ 程序实现了 Wolfe-Powell 准则的一个...

    最优化算法C/C++源码

    9. **Conjugate Gradient法**:共轭梯度法是解决对称正定线性系统的有效方法,也可以用于非线性优化问题。它利用梯度的共轭性质,减少迭代次数。 以上每一种算法都有其适用的场景和局限性。在实际应用中,需要根据...

    C++shuzi.rar_数值优化算法

    常见的数值优化方法包括梯度下降法、牛顿法、拟牛顿法、共轭梯度法、线性规划、二次规划等。 2. C++经典数值算法: C++作为一门强类型、静态类型的编程语言,其高效的性能和丰富的库支持使得它非常适合实现数值计算...

    精通MATLAB最优化计算源代码,最优化计算方法及其matlab程序实现,matlab

    MATLAB提供了多种内置的最优化工具箱,如`fminunc`(非线性无约束优化)、`fmincon`(非线性有约束优化)和`intlinprog`(整数线性规划)等,这些工具箱涵盖了各种优化问题的求解。 二、MATLAB优化算法 1. **梯度...

    常用数值算法与程序(C++版)

    5. **最优化**:如梯度下降法、拟牛顿法、共轭梯度法,用于求解无约束或约束优化问题。 6. **数值解微分方程**:如欧拉方法、龙格-库塔方法,用于求解常微分方程初值问题。 7. **随机数生成和统计分析**:包含各种...

    matlab 精通最优化计算

    这涉及理解优化理论,比如梯度法、牛顿法、共轭梯度法等,并用 MATLAB 的矩阵运算和控制结构实现。 7. **优化报告和结果分析**:在完成优化过程后,通常需要分析结果,包括最佳解、函数值、收敛历史等。MATLAB 提供...

    1_JorgeNocedal_数值优化_数值优化必备书籍_

    非线性方程组的求解也是数值优化的一部分,书中介绍了高斯-塞德尔迭代和共轭梯度法等解决这类问题的技术。这些方法在科学计算中有着广泛的应用。 除了基本的优化算法,Nocedal和Wright还讨论了全局优化和随机搜索...

    非线性优化nlopt.7z

    无约束优化算法如梯度下降法、牛顿法、拟牛顿法(如BFGS和L-BFGS)、共轭梯度法等,而有约束优化则包括了拉格朗日乘子法、罚函数法和屏障方法等。这些算法的选择取决于问题的具体性质,例如目标函数是否具有已知的...

    quadprog c++程序

    3. **优化算法**:在每一步迭代中,可能需要用到求解线性系统的算法,如高斯-塞德尔迭代或共轭梯度法,来更新决策变量。 4. **错误处理与调试**:为了确保代码的健壮性,需要处理可能的错误情况,如输入数据无效、...

    数值分析的常用C++代码

    C++中常见的优化算法有梯度下降法、牛顿法、拟牛顿法(如BFGS和L-BFGS)、共轭梯度法等。对于约束优化问题,可以使用拉格朗日乘子法或罚函数法。 5. **概率统计**:在C++中,数值分析也应用于概率统计,如蒙特卡洛...

    C++ 二次规划源码 quadprog++.rar

    这个库可能采用了诸如内点法、梯度下降法或者共轭梯度法等优化算法来求解问题。 C++作为一个强大的系统编程语言,它的标准库并没有内置二次规划的解决方案。因此,像"quadprog++"这样的第三方库就显得尤为重要。...

    最速下降法,最速下降法和牛顿法的优缺点,C,C++源码.zip

    同时,为了处理非凸或非光滑函数,可以采用拟牛顿法或共轭梯度法等改进策略。 在提供的压缩包文件中,包含了C和C++源码,这可能是实现上述优化算法的具体示例。通过分析和理解这些源码,可以深入掌握这两种优化方法...

Global site tag (gtag.js) - Google Analytics