谁有拟牛顿法、共轭梯度法、Powell方法之一编写无约束优化程序??
谁有拟牛顿法、共轭梯度法、Powell方法之一编写无约束优化程序??matlab的。。拜谢了~~~~~~~~~~~~~~~~~ 拟牛顿法、共轭梯度法论坛以前的帖子中好象有，自己找一下吧 最优化-无约束共轭梯度法程序（c++）
/////////////////////////////////////////
/////vector.h头文件/////
/////定义向量及其基本运算/////
/////////////////////////////////////////
#include
#defineMAXLENGTH10
//向量定义
typedefstruct{
inttag;//行列向量标志。行向量为0，列向量为1。
intdimension;//向量的维数
doubleelem;//向量的元素
}vector;
vectorvecCreat(inttag,intn){
//建立维数为n的向量
vectorx;
x.tag=tag;
x.dimension=n;
for(inti=0;i
#include"vector.h"
#defineSIZE10
#defineMAX1e300
doublemin(doublea,doubleb){
//求两个数最小
returna}
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;itempPnt1=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;
for(inti=0;ix=vx.elem;
//----------约束目标函数--------------
//returnx*x+x*x;
//----------无约束正定函数--------------
//returnx*x+4*x*x+9*x*x-2*x+18*x;//例一
//----------无约束非二次函数--------------
//return(1-x)*(1-x)+(1-x)*(1-x)+(x*x-x)*(x*x-x
)+(x*x-x)*(x*x-x)+(x*x-x)*(x*x-x);//例二
}
doublevecFun_Si(vectorvx){
//不等式约束函数
doublex;
for(inti=0;ix=vx.elem;
returnx+1;//不等式约束函数
}
doublevecFun_Hi(vectorvx){
//等式约束函数
doublex;
for(inti=0;ix=vx.elem;
returnx+x-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;isum1.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(g2f1=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<<");";//输出本次结果
//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;
};
}
}
自己改写吧 另外《MATLAB 6.0程序设计与实例应用》好像有matlab的 非常感谢老兄们的回复！
to 2楼：我搜过了，怎么没搜到呢？
to 3楼：看不太懂
to 4楼：我找过了，那个共轭梯度法不是用来无约束优化的 看帖子
http://forum.vibunion.com/thread-732-1-1.html
4楼讲的就是无约束非线性规划问题
[此贴子已经被作者于2005-11-19 20:41:54编辑过]
这有一个，看能不能用
function = fminu(FUN,x,OPTIONS,GRADFUN,varargin)
%多元函数极值拟牛顿法，适用于光滑函数优化�
%x=fminu('fun',x0)牛顿法求多元函数y=f(x)在x0出发的局部极小值点�
% 这里 x,x0均为向量。
%x=fminu('fun',x0,options)输入options是优化中算法参数向量设置,
% 用help foptions可看到各分量的含义。
%x=fminu('fun',x0,options,'grad') grad给定f(x)的梯度函数表达式？杉涌旒扑闼俣取�
%x=fminu('fun',x0,options,'grad',p1,p2…)p1,p2,..是表示fun的M函数中的参数。
%例题 求f(x)=100(x2-x1^2)^2+(1-x1)^2在[-1.2,1]附近的极小值。
% 先写M函数optfun1.m
% function f=optfun1(x)
% f=100*(x(2)-x(1).^2).^2+(1-x(1)).^2
% 求解
% clear;
% x=[-1.2,1];
% =fminu('optfun1',x,options);
% x,options(8)
%
%FMINU Finds the minimum of a function of several variables.
% X=FMINU('FUN',X0) starts at the matrix X0 and finds a minimum to the
% function which is described in FUN (usually an M-file: FUN.M).
% The function 'FUN' should return a scalar function value: F=FUN(X).
%
% X=FMINU('FUN',X0,OPTIONS) allows a vector of optional parameters to
% be defined. OPTIONS(1) controls how much display output is given; set
% to 1 for a tabular display of results, (default is no display: 0).
% OPTIONS(2) is a measure of the precision required for the values of
% X at the solution. OPTIONS(3) is a measure of the precision
% required of the objective function at the solution.
% For more information type HELP FOPTIONS.
%
% X=FMINU('FUN',X0,OPTIONS,'GRADFUN') enables a function'GRADFUN'
% to be entered which returns the partial derivatives of the function,
% df/dX, at the point X: gf = GRADFUN(X).
%
% X=FMINU('FUN',X0,OPTIONS,'GRADFUN',P1,P2,...) passes the problem-
% dependent parameters P1,P2,... directly to the functions FUN
% and GRADFUN: FUN(X,P1,P2,...) and GRADFUN(X,P1,P2,...). Pass
% empty matrices for OPTIONS, and 'GRADFUN' to use the default
% values.
%
% =FMINU('FUN',X0,...) returns the parameters used in the
% optimization method. For example, options(10) contains the number
% of function evaluations used.
%
% The default algorithm is the BFGS Quasi-Newton method with a
% mixed quadratic and cubic line search procedure.
% Copyright (c) 1990-98 by The MathWorks, Inc.
% Revision:1.27 Date:1997/11/2901:23:11
% Andy Grace 7-9-90.
% ------------Initialization----------------
XOUT=x(:);
nvars=length(XOUT);
if nargin < 2, error('fminu requires two input arguments');end
if nargin < 3, OPTIONS=[]; end
if nargin < 4, GRADFUN=[]; end
% Convert to inline function as needed.
if ~isempty(FUN)
= fcnchk(FUN,length(varargin));
if ~isempty(msg)
error(msg);
end
else
error('FUN must be a function name or valid expression.')
end
if ~isempty(GRADFUN)
= fcnchk(GRADFUN,length(varargin));
if ~isempty(msg)
error(msg);
end
else
gradfcn = [];
end
f = feval(funfcn,x,varargin{:});
n = length(XOUT);
GRAD=zeros(nvars,1);
OLDX=XOUT;
MATX=zeros(3,1);
MATL=;
OLDF=f;
FIRSTF=f;
=optint(XOUT,f,OPTIONS);
CHG = 1e-7*abs(XOUT)+1e-7*ones(nvars,1);
SD = zeros(nvars,1);
diff = zeros(nvars,1);
PCNT = 0;
how = '';
OPTIONS(10)=2; % Function evaluation count (add 1 for last evaluation)
status =-1;
while status ~= 1
% Work Out Gradients
if isempty(gradfcn) | OPTIONS(9)
OLDF=f;
% Finite difference perturbation levels
% First check perturbation level is not less than search direction.
f = find(10*abs(CHG)>abs(SD));
CHG(f) = -0.1*SD(f);
% Ensure within user-defined limits
CHG = sign(CHG+eps).*min(max(abs(CHG),OPTIONS(16)),OPTIONS(17));
for gcnt=1:nvars
XOUT(gcnt,1)=XOUT(gcnt)+CHG(gcnt);
x(:) = XOUT;
f = feval(funfcn,x,varargin{:});
GRAD(gcnt)=(f-OLDF)/(CHG(gcnt));
if f < OLDF
OLDF=f;
else
XOUT(gcnt)=XOUT(gcnt)-CHG(gcnt);
end
end
% Try to set difference to 1e-8 for next iteration
% Add eps for machines that can't handle divide by zero.
CHG = 1e-8./(GRAD + eps);
f = OLDF;
OPTIONS(10)=OPTIONS(10)+nvars;
% Gradient check
if OPTIONS(9) == 1
GRADFD = GRAD;
x(:)=XOUT;
GRAD(:) = feval(gradfcn,x,varargin{:});
if isa(gradfcn,'inline')
graderr(GRADFD, GRAD, formula(gradfcn));
else
graderr(GRADFD, GRAD, gradfcn);
end
OPTIONS(9) = 0;
end
else
OPTIONS(11)=OPTIONS(11)+1;
x(:)=XOUT;
GRAD(:) = feval(gradfcn,x,varargin{:});
end
%---------------Initialization of Search Direction------------------
if status == -1
SD=-GRAD;
FIRSTF=f;
OLDG=GRAD;
GDOLD=GRAD'*SD;
% For initial step-size guess assume the minimum is at zero.
OPTIONS(18) = max(0.001, min());
if OPTIONS(1)>0
disp();
end
XOUT=XOUT+OPTIONS(18)*SD;
status=4;
if OPTIONS(7)==0; PCNT=1; end
else
%-------------Direction Update------------------
gdnew=GRAD'*SD;
if OPTIONS(1)>0,
num=;
end
if (gdnew>0 & f>FIRSTF)|~finite(f)
% Case 1: New function is bigger than last and gradient w.r.t. SD -ve
% ...interpolate.
how='inter';
=cubici1(f,FIRSTF,gdnew,GDOLD,OPTIONS(18));
if stepsize<0|isnan(stepsize), stepsize=OPTIONS(18)/2; how='C1f '; end
if OPTIONS(18)<0.1&OPTIONS(6)==0
if stepsize*norm(SD)<eps
stepsize=exp(rand(1,1)-1)-0.1;
how='RANDOM STEPLENGTH';
status=0;
else
stepsize=stepsize/2;
end
end
OPTIONS(18)=stepsize;
XOUT=OLDX;
elseif f<FIRSTF
=cubici3(f,FIRSTF,gdnew,GDOLD,OPTIONS(18));
sk=(XOUT-OLDX)'*(GRAD-OLDG);
if sk>1e-20
% Case 2: New function less than old fun. and OK for updating HESS
% .... update and calculate new direction.
how='';
if gdnew<0
how='incstep';
if newstep<OPTIONS(18), newstep=2*OPTIONS(18)+1e-5; how=; end
OPTIONS(18)=min(),1+sk+abs(gdnew)+max(), (1.2+0.3*(~OPTIONS(7)))*abs(newstep)]);
else % gdnew>0
if OPTIONS(18)>0.9
how='int_st';
OPTIONS(18)=min();
end
end %if gdnew
=updhess(XOUT,OLDX,GRAD,OLDG,HESS,OPTIONS);
gdnew=GRAD'*SD;
OLDX=XOUT;
status=4;
% Save Variables for next update
FIRSTF=f;
OLDG=GRAD;
GDOLD=gdnew;
% If mixed interpolation set PCNT
if OPTIONS(7)==0, PCNT=1; MATX=zeros(3,1); MATL(1)=f; end
elseif gdnew>0 %sk<=0
% Case 3: No good for updating HESSIAN .. interpolate or halve step length.
how='inter_st';
if OPTIONS(18)>0.01
OPTIONS(18)=0.9*newstep;
XOUT=OLDX;
end
if OPTIONS(18)>1, OPTIONS(18)=1; end
else
% Increase step, replace starting point
OPTIONS(18)=max(),0.5*OPTIONS(18)]);
how='incst2';
OLDX=XOUT;
FIRSTF=f;
OLDG=GRAD;
GDOLD=GRAD'*SD;
OLDX=XOUT;
end % if sk>
% Case 4: New function bigger than old but gradient in on
% ...reduce step length.
else %gdnew<0 & F>FIRSTF
if gdnew<0&f>FIRSTF
how='red_step';
if norm(GRAD-OLDG)<1e-10; HESS=eye(nvars); end
if abs(OPTIONS(18))<eps
SD=norm(nvars,1)*(rand(nvars,1)-0.5);
OPTIONS(18)=abs(rand(1,1)-0.5)*1e-6;
how='RANDOM SD';
else
OPTIONS(18)=-OPTIONS(18)/2;
end
XOUT=OLDX;
end %gdnew>0
end % if (gdnew>0 & F>FIRSTF)|~finite(F)
XOUT=XOUT+OPTIONS(18)*SD;
if isinf(OPTIONS(1))
disp()
elseif OPTIONS(1)>0
disp(num)
end
end %----------End of Direction Update-------------------
% Check Termination
if max(abs(SD))<2*OPTIONS(2) & (-GRAD'*SD) < 2*OPTIONS(3)
if OPTIONS(1) > 0
disp('Optimization Terminated Successfully')
disp(' Search direction less than 2*options(2)')
disp(' Gradient in the search direction less than 2*options(3)')
disp([' NUMBER OF FUNCTION EVALUATIONS=', int2str(OPTIONS(10))]);
end
status=1;
elseif OPTIONS(10)>OPTIONS(14)
if OPTIONS(1) >= 0
disp('Maximum number of function evaluations exceeded;')
disp(' increase options(14).')
end
status=1;
else
% Line search using mixed polynomial interpolation and extrapolation.
if PCNT~=0
while PCNT > 0 & OPTIONS(10) <= OPTIONS(14)
x(:) = XOUT;
f = feval(funfcn,x,varargin{:});
OPTIONS(10)=OPTIONS(10)+1;
=searchq(PCNT,f,OLDX,MATL,MATX,SD,GDOLD,OPTIONS(18), how);
OPTIONS(18)=steplen;
XOUT=OLDX+steplen*SD;
end
if OPTIONS(10)>OPTIONS(14)
if OPTIONS(1) >= 0
disp('Maximum number of function evaluations exceeded;')
disp(' increase options(14).')
end
status=1;
end
else
x(:)=XOUT;
f = feval(funfcn,x,varargin{:});
OPTIONS(10)=OPTIONS(10)+1;
end
end
end
x(:)=XOUT;
f = feval(funfcn,x,varargin{:});
if f > FIRSTF
OPTIONS(8) = FIRSTF;
x(:)=OLDX;
else
OPTIONS(8) = f;
end 这个我也找到了 ～不还是看不懂啊～ 为什么哪些例题上程序那么简单啊～真的要编怎么又那么多啊 ～ 呵呵，真辛苦各位拉！ 谢谢
页:
[1]