Fletcher-Reeves法サンプルプログラム

無制約最適化問題で使われる（？）Fletcher-Reeves法のプログラム例を載せます．アルゴリズム習得用で組んだものです．
手法自体は適当に調べました．
http://people.equars.com/~marco/poli/phd/node54.html

例：天井からばね定数1のばねによって連結された4つの重量1の質点の平衡点を求める．ポテンシャルエネルギーを目的関数とする最小化問題．

global k mg L
mg=1;%おもりの重量
L=1;%ばね自然長
y=[1 2 3 4]';%初期値
k=1;%ばね定数
　　gF=zeros(4,2);%betaを計算するために2列用意
gF(:,1)=gradF(y);
%↓目的関数Fの定義
F=@(y) (k/2)**1/(gF(:,1)'*gF(:,1));%betaが定義される
　　　　　　s=-gF(:,2)+beta*s;
　　　　end
　　　　%y_{i+1}=y_{i}+alpha*s;
　　　　%ここでalhpaはF(y_{i+1})を最小化する値
　　　　%alphaを少しずつ大きくしながら探索する
　　　　num_roop2=100;%Fを最小化するalphaに到着するのに十分な値
　　　　step=0.1;%探索の歩幅
　　　　F_temp=zeros(2,1);%一時的な目的関数
　　　　F_temp(1)=F(y);%出発点
　　　　for ii=1:num_roop2;%ループするよ
　　　　　　F_temp(2)=F(y+step*ii*s);%1歩先のF
　　　　　　if F_temp(2)>=F_temp(1);%もし1歩先のFの方が大きいなら
　　　　　　　　break%繰り返し計算を終了する
　　　　　　else%まだFが小さくなりそうなら
　　　　　　　　F_temp(1)=F_temp(2);%参照するFを置き換える
　　　　　　end
　　　　end
　　　　y=y+step*(ii-1)*s;%この半直線上でFを最小化するy
　　　　gF(:,1)=gF(:,2);%繰り返しのために勾配を置き換え
end
if jj==num_roop1;%もし収束していないなら
　　fprintf('計算は収束していない')%その旨を表示
end

====================================================
====gradF.m=========================================
====================================================
function out=gradF(y)
　　global k mg L
　　out(1)=k*(y(1)  　　 -L)-k*(y(2)-y(1)-L)-mg;
　　out(2)=k*(y(2)-y(1)-L)-k*(y(3)-y(2)-L)-mg;
　　out(3)=k*(y(3)-y(2)-L)-k*(y(4)-y(3)-L)-mg;
　　out(4)=k*(y(4)-y(3)-L)　　　　　　　　　　　-mg;

これによって，最適解[5 9 12 14]が得られます．もちろん力のつり合い等で求めた解と一致します．

*1:y(1)-L)^2+(y(2)-y(1)-L)^2+(y(3)-y(2)-L)^2+(y(4)-y(3)-L)^2)-mg*sum(y); count=0;%カウント(一周目かどうかが分かればよい) num_roop1=100;%収束するのに十分な値 for jj=1:num_roop1; 　　gF(:,2)=gradF(y);%勾配 　　if max(abs(gF(:,2)))<=1e-3;%収束判定 　　　　break 　　end 　　　　 %探索方向sを決定する 　　　　if count==0;%もし一周目なら 　　　　　　s=-gF(:,2);%betaは0; 　　　　else%2周目以降なら 　　　　　　beta=(gF(:,2)'*gF(:,2