最优化_三等分法+黄金分割法+牛顿法
一、实验目的
-
掌握一维优化方法的集中算法；
-
编写三分法算法
-
编写黄金分割法算法
-
编写牛顿法算法
二、系统设计
三分法
1.编程思路：
三分法用于求解单峰函数的最值。对于单峰函数，在区间内用两个mid将区间分成三份，这样的查找算法称为三分查找,也就是三分法。在区间[a,b]内部取n=2个内等分点，区间被分为n+1=3等分，区间长度缩短率\(=\frac{1}{3}\) .
各分点的坐标为 \(x_k=a+\frac{b-a}{n+1}\cdot k\ \ (k=1,2)\) ，
然后计算出\(x_1,x_2,\cdots ;y_1,y_2,\cdots;\)
找出\(y_{min}=min\{y_k,k=1,2\}\) ，新区间\((a,b)\Leftarrow (x_{m-1},x_{m+1})\) .
coding中，建立\(left,mid1,mid2,right\) 四个变量用于计算，用新的结果赋值给旧区间即可。
2.算法描述
function [left]=gridpoint(left,right,f)
epsilon=1e-5; %给定误差范围
while((left+epsilon)<right) %检查left,right区间精度
margin=(right-left)/3; %将区间三等分,每小段长度=margin
m1=left+margin; %left-m1-m2-right，三等分需要两个点
m2=m1+margin; %m2=left+margin+margin
if(f(m1)<=f(m2))
right=m2; %离极值点越近，函数值越小(也有可能越大，视函数而定)。
else %当f(m1)>f(m2),m2离极值点更近。缩小区间范围，逼近极值点
left=m1; %所以令left=m1.
end
end %这是matlab的.m文件，不用写return.
黄金分割法
1.编程思路
三分法进化版，区间长度缩短率\(\approx 0.618\).
在区间\([a,b]\)上取两个内试探点，\(p_i,q_i\) 要求满足下面两个条件：
\(1.[a_i,q_i]与[p_i,b_i]的长度相同，即b_i-p_i=q_i-a_i;\)
\(2.区间长度的缩短率相同，即b_{i+1}-a_{i+1}=t(b_i-a_i)]\)
2.算法描述
自己编写的：
function [s,func_s,E]=my_golds(func,left,right,delta)
tic
%输入: func:目标函数，left,right:初始区间两个端点
% delta：自变量的容许误差
%输出: s,func_s:近似极小点和函数极小值
% E=[ds,dfunc] ds,dfunc分别为s和dfunc的误差限
%0.618法的改进形式：每次缩小区间时，同时比较两内点和两端点处的函数值。
%当左边第一个或第二个点是这四个点中函数值最小的点时，丢弃右端点，构成新的搜索区间
%否则，丢弃左端点。
%left---------|--mu lambda，mu为了两个试探点，且lambda<mu
% lambda--|------right
%根据定理2.1.4
%若f(lambda_k)<=f(mu_k)
% 则left_(k+1)=a_k, right_(k+1)=mu_k
%若f(lambda_k)>f(mu_k)
% 则left_(k+1)=lambda_k, right_(k+1)=right_k
t=(sqrt(5)-1)/2;h=right-left; %区间长度缩短率t
lambda=left+(1-t)*h;
mu=left+t*h;
func_left=feval(func,left);
func_right=feval(func,right);
func_lambda=feval(func,lambda);
func_mu=feval(func,mu);
k=1;
G=[func_left,func_lambda,func_mu,func_right];
while(abs(right-left)>delta)
func_t=min(G);
if(func_t==func_left|func_t==func_lambda)
%当func_left,func_lambda较小时，丢弃右端点，即t<3，转步4
right=mu;
func_right=func_mu;
mu=lambda;
func_mu=func_lambda;
h=right-left;%漏行！！！！！！！！！！！！
lambda=left+(1-t)*h;
func_lambda=feval(func,lambda);%漏行！！！！！！！！
else %当func_mu,func_right较小时，丢弃左端点，即t不<3,转步3
left=lambda;
func_left=func_lambda;
lambda=mu;
func_lambda=func_mu;
h=right-left;
mu=left+t*h;
func_mu=feval(func,mu);
end
k=k+1;
G=[func_left,func_lambda,func_mu,func_right];
end
ds=abs(right-left);dfunc=abs(func_right-func_left);
if(func_lambda<=func_mu)
s=lambda;func_s=func_lambda;
else
s=mu;func_s=func_mu;
end
E=[ds,dfunc];toc;
书上例程：
%黄金分割法（0.618）
function [s,phi_s,k,G,E]=golds(phi,left,right,delta,epsilon)
tic;
%输入：phi:目标函数；left,right:搜索区间两个端点
% delta,epsilon:自变量和函数值的容许误差
%输出：s,phi_s:近似极小点和极小值，G为n*4矩阵
% 其第k行分别为left,m1,m2,right的第k次迭代值[ak,m1k,m2k,nk]
% E=[ds,dphi] ,分别为s,phis的误差限
t=(sqrt(5)-1)/2;
h=right-left;
phi_left=feval(phi,left);phi_right=feval(phi,right);
%[y1,y2,...] = feval(fhandle,x1,x2,...,xn)
m1=left+(1-t)*h;m2=left+t*h;
phi_m1=feval(phi,m1);phi_m2=feval(phi,m2);
k=1;
G(k,:)=[left,m1,m2,right];
while((abs(phi_m2-phi_m1)>epsilon)|(h>delta))
if(phi_m1<phi_m2)
right=m2;phi_right=phi_m2;
m2=m1;phi_m2=phi_m1;
h=right-left;
m1=left+(1-t)*h;phi_m1=feval(phi,m1);
else
left=m1;phi_left=phi_m1;
m1=m2;phi_m1=phi_m2;
h=right-left;
m2=left+t*h;phi_m2=feval(phi,m2);
end
k=k+1;G(k,:)=[left,m1,m2,right];
end
ds=abs(right-left);dphi=abs(phi_right-phi_left);
if(phi_m1<=phi_m2)
s=m1;phi_s=phi_m1;
else
s=m2;phi_s=phi_m2;
end
E=[ds,dphi] ;toc;
牛顿法
1.编程思路
牛顿法用于求解无约束优化问题，基本思想是用迭代点\(x_k\)处的一阶导数（梯度）和二阶导数（Hessian阵）对目标函数进行二次函数近似，然后用二次模型的极小点作为新的迭代点，并不断重复，直到求得满足精度要求的近似极小点。
其迭代公式为：\(x_{k+1}=x_k-G^{-1}_kg_k\)
【书上的推导看得脑壳痛，还是PPT简单易懂】
2.算法描述
DIY：
function [min,k]=NR1(f,x0,epsilon,var)
%x0表示起始点
%f为函数
%epsilon表示误差范围
%var表示函数变量名
error=epsilon+1;
k=0;
diff(f)
while((error>epsilon)&(k<1000))
if(subs(diff(f),var,x0)==0)
break
end
x1=x0-subs(diff(f),var,x0)/subs(diff(f,2),var,x0);
x1=double(x1);
error=abs(x1-x0);
x0=x1;k=k+1;
end
min=x0;
例程：
-
修正牛顿法
-
阻尼牛顿法
function [x,val,k]=dampnm(fun,gfun,Hess,x0) %功能：用阻尼牛顿法求解无约束问题：min f(x) %输入：x0是初始点，fun,gfun,Hess分别求目标函数值，梯度，Hessian阵的函数 %输出：x,val分别为近似最优点和最优值，k为迭代次数 maxk=100; %给出最大迭代次数 rho=0.55;sigma=0.4; k=0;epsilon=1e-5; while(k<maxk) gk=feval(gfun,x0); %计算梯度 Gk=feval(Hess,x0); %计算Hessian阵 dk=-Gk/gk; %解方程组Gk*dk=-gk,计算搜索方向 if(norm(gk)<epsilon)%n=norm(A,p),计算矩阵范数 break; end m=0;mk0; while(m<20) %用Armijo搜索求步长 if(feval(fun,x0+rho^m*dk)<feval(fun,x0)+sigma*rho^m*gk'*dk) mk=m; break; end m=m+1; end x0=x0+rho^mk*dk; k=k+1; end x=x0; val=feval(fun,x);举例：求解无约束优化问题：\(\mathop{min}\limits_{x\in\mathbb{R}^2} f(x)=100(x^2_1-x_2)^2+(x_1-1)^2\)
minx∈R2f(x)=100(x12−x2)2+(x1−1)2 \mathop{min}\limits_{x\in\mathbb{R}^2} f(x)=100(x^2_1-x_2)^2+(x_1-1)^2先建立两个分别计算目标函数和梯度的.m文件
function f=fun(x) function g=gfun(x)再建立求\(Hessian\)阵的.m文件
function He=Hess(x) n=length(X); He=zeros(n,n); He=[1200*x(1)^2-400*x(2)+2,-400*x(1); -400*x(1), 200 ];终止准则取为\(\parallel \nabla f(x_k)\parallel\leqslant10^{-5}\)
∥∇f(xk)∥⩽10−5 \parallel \nabla f(x_k)\parallel\leqslant10^{-5}程序调用方式为：
\(x0=[-1.2,1]'\)
\([x,val,k]=dampnm('fun','gfun','Hess',x0)\)
三、测试
三分法
黄金分割法
牛顿法
四、总结
matlab 不用声明变量，科学计算的时候很方便，不过没声明的变量编译默认为0，可能会带来奇怪的bug.