最优化无约束共轭梯度法程序c.docx
- 文档编号:14953696
- 上传时间:2023-06-28
- 格式:DOCX
- 页数:16
- 大小:18.30KB
最优化无约束共轭梯度法程序c.docx
《最优化无约束共轭梯度法程序c.docx》由会员分享,可在线阅读,更多相关《最优化无约束共轭梯度法程序c.docx(16页珍藏版)》请在冰点文库上搜索。
最优化无约束共轭梯度法程序c
最优化-无约束共轭梯度法程序(c++)
////////////////////////////////////////
/////vector.h头文件/////
/////定义向量及其基本运算/////
/////////////////////////////////////////
#include
#defineMAXLENGTH10
//向量定义
typedefstructvector
{
inttag;//行列向量标志。
行向量为0,列向量为1。
intdimension;//向量的维数
doubleelem[MAXLENGTH];//向量的元素
};
vectorvecCreat(inttag,intn)
{//建立维数为n的向量
vectorx;
x.tag=tag;
x.dimension=n;
for(inti=0;i { cout<<"请输入第"< "; cin>>x.elem; } returnx; } doublevecMultiply(vectora,vectorb) {//行向量与列向量乘法运算 if((a.tag! =b.tag)&&(a.dimension==b.dimension)) {//相乘条件 doublec=0; for(inti=0;i c+=a.elem*b.elem; returnc; } } vectorvecAdd(vectora,vectorb) { //向量加法运算 if((a.tag==b.tag)&&(a.dimension==b.dimension)) {//相加条件 vectorc; c.dimension=a.dimension; c.tag=a.tag; for(inti=0;i c.elem=a.elem+b.elem; returnc; } } vectorvecConvert(vectora) { //向量转置 if(a.tag==0)a.tag=1; if(a.tag==1)a.tag=0; returna; } doublevecMole(vectora) { //求向量模 doublesum=0; for(inti=0;i sum+=a.elem*a.elem; sum=sqrt(sum); returnsum; } vectornumMultiply(doublen,vectora) { //数乘向量 for(inti=0;i a.elem=n*a.elem; returna; } voidshowPoint(vectorx,doublef) { //显示点,解. cout<<"---x=("; for(inti=0;i { cout< if(i! =x.dimension-1) cout<<","; } cout<<")---f(x)="< cout< } ///////////////////////////////////////////////////////// //////////////function.h头文件///////////////////////// ////////函数及其导数的设置均在此文件完成//////////////// ///////////////////////////////////////////////////////// //*无约束问题*// //目标函数--在vecFun(vectorvx)中修改,x1改成x[1]// //x2改成x[2],依此类推...// ///////////////////////////////////////////////////////// #include #include"vector.h" #defineSIZE10 #defineMAX1e300 doublemin(doublea,doubleb) { //求两个数最小 returna a: b; } doublemax(doublea,doubleb) { //求两个数最大 returna>b? a: b; } vectorvecGrad(double(*pf)(vectorx),vectorpointx) { //求解向量的梯度 //采用理查德外推算法计算函数pf在pointx点的梯度grad; vectorgrad; grad.tag=1;grad.dimension=pointx.dimension;//初始化偏导向量 vectortempPnt1,tempPnt2;//临时向量 doubleh=1e-3; for(inti=0;i tempPnt1=tempPnt2=pointx; tempPnt1.elem+=0.5*h; tempPnt2.elem-=0.5*h; grad.elem=(4*(pf(tempPnt1)-pf(tempPnt2)))/(3*h); tempPnt1.elem+=0.5*h; tempPnt2.elem-=0.5*h; grad.elem-=(pf(tempPnt1)-pf(tempPnt2))/(6*h); } returngrad; } doublevecFun(vectorvx) { //最优目标多元函数 doublex[SIZE]; for(inti=0;i x[i+1]=vx.elem; //----------约束目标函数-------------- //returnx[1]*x[1]+x[2]*x[2]; //----------无约束正定函数-------------- //returnx[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]);//例二 } doublevecFun_Si(vectorvx) { //不等式约束函数 doublex[SIZE]; for(inti=0;i x[i+1]=vx.elem; returnx[1]+1;//不等式约束函数 } doublevecFun_Hi(vectorvx) { //等式约束函数 doublex[SIZE]; for(inti=0;i x[i+1]=vx.elem; returnx[1]+x[2]-2;//等式约束函数 } doublevecFun_S(vectorx,doublev,doublel,doubleu) { //约束问题转化为无约束问题的增广目标函数F(x,v,l,u) doublesum1=0,sum2=0,sum3=0;//分别定义三项的和 //sum1 doubletemp=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 returnvecFun(x)+sum1-sum2+sum3; } vectorvecFunD_S(vectorx,doublev,doublel,doubleu) { //利用重载函数实现目标函数F(x,v,l,u)的导数 //约束问题转化为无约束问题的增广目标函数F(x,v,l,u)的导函数 vectorsum1,sum2,sum3;//分别定义三项导数的和 //sum1 sum1.dimension=x.dimension;sum1.tag=1; sum2.dimension=x.dimension;sum2.tag=1; sum3.dimension=x.dimension;sum3.tag=1; doubletemp=max(0,v-2*u*vecFun_Si(x)); if(temp==0){ for(inti=0;i sum1.elem=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 returnvecAdd(vecAdd(vecGrad(vecFun,x),sum1),vecAdd(sum2,sum3)); } /////////////////////////////////////////////////// /////nonrestrict.h头文件///// /////包含无约束问题的求解函数: 直线搜索///// /////共轭梯度法,H终止准则///// /////////////////////////////////////////////////// #include"restrict.h" vectorlineSearch(vectorx,vectorp,doublet) { //从点x沿直线p方向对目标函数f(x)作直线搜索得到极小点x2,t为初始步长。 vectorx2; doublel=0.1,n=0.4;//条件1和2的参数 doublea=0,b=MAX;//区间 doublef1=0,f2=0,g1=0,g2=0; inti=0;//迭代次数 do{ if(f2-f1>l*t*g1){b=t;t=(a+b)/2;}//改变b,t循环 do{ if(g2 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){returnx2;}//【直线搜索进入无 限跌代,则此次跌代结束。 返回当前最优点】 f2=vecFun(x2);//f2=f(x2) g2=vecMultiply(vecConvert(vecGrad(vecFun,x2)),p);//g2=∨f(x2 )*p //cout<<"("< .elem[3]<<");";//输出本次结果 //cout<<"t="< //cout<<"f(x)="< } while(g2 } while(f2-f1>l*t*g1);//不满足条件i,则改变b,t循环 returnx2; } intHimmulblau(vectorx0,vectorx1,doublef0,doublef1) { //H终止准则.给定Xk,Xk+1,Fk,Fk+1,判断Xk+1是否是极小点.返回1是极小,否则返回0 doublec1,c2,c3;//定义并初始化终止限 c1=c2=10e-5; c3=10e-4; if(vecMole(vecGrad(vecFun,x1)) if(vecMole(vecAdd(x1,numMultiply(-1,x0)))/(vecMole(x0)+1) if(fabs(f1-f0)/(fabs(f0)+1) return1; return0; } voidFletcher_Reeves(vectorxx,vector&minx,double&minf) { //Fletcher-Reeves共轭梯度法. //给定初始点xx.对vecFun函数求得最优点x和最优解f,分别用minx和minf返回 intk=0,j=0;//迭代次数,j为总迭代次数 doublec=10e-1;//终止限c intn=xx.dimension;//问题的维数 doublef0,f,a;//函数值f(x),a为使p方向向量共轭的系数 vectorg0,g;//梯度g0,g vectorp0,p;//搜索方向P0,p vectorx,x0;//最优点和初始点 doublet=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<<"*总迭代次数: "< minx=x;minf=f; break; } else{//不满足H准则 cout<<"*第"< 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; } elseif(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" voidprintSelect() { cout<<"**************最优化方法***************"< cout<<"*****************************************"< cout<<"***请选择: ***"< cout<<"***1.无约束最小化问题(共轭梯度法)***"< cout<<"***2.约束最小化问题(H乘子法)***"< cout<<"***3.退出(exit)***"< cout<<"*****************************************"< } voidsub1(){ charkey; cout<<"--------------------共轭梯度法求无约束最小化问题----------------"<< endl; cout<<"步骤: "< cout<<"◆1.输入f(x)及其导数.(参考function.h文件说明)"< cout<<"-----确认是否已输入<目标函数>及其<导数>? (Y/N)"; cin>>key; if(key=='N'||key=='n')return; elseif(key=='Y'||key=='y'){ vectorx0;//起始点 intm; cout<<"◆2.起始点X0初始化"< "; cin>>m; x0=vecCreat(1,m); vectorminx;//求得的极小点 doubleminf;//求得的极小值 Fletcher_Reeves(x0,minx,minf); cout<<"◆最优解及最优值: "; showPoint(minx,minf); } } voidsub2(){ charkey; cout<<"------------------H乘子法求约束最小化问题-----------------"< ; cout<<"步骤: "< cout<<"◆1.输入f(x)及其导数.(参考function.h文件说明)"< cout<<"-----确认是否已输入<目标函数及导数>和<约束条件>? (Y/N)"; cin>>key; if(key=='N'||key=='n')return; elseif(key=='Y'||key=='y'){ vectorx0;//起始点 intm; cout<<"◆2.起始点X0初始化"< "; cin>>m; x0=vecCreat(1,m); vectorminx;//求得的极小点 doubleminf;//求得的极小值 Hesternes(x0,minx,minf); showPoint(minx,minf); } } voidmain(){ intmark=1; while(mark){ printSelect(); intsel; cin>>sel; switch(sel){ case1: sub1(); break; case2: sub2(); break; case3: mark=0; break; }; } }
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 优化 无约束 共轭 梯度 程序