最優化_三等分法+黃金分割法+牛頓法
一、實驗目的
-
掌握一維優化方法的集中算法；
-
編寫三分法算法
-
編寫黃金分割法算法
-
編寫牛頓法算法
二、系統設計
三分法
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.