読者です 読者をやめる 読者になる 読者になる

MyEnigma

とあるエンジニアのブログです。#Robotics #Programing #C++ #Python #MATLAB #Vim #Mathematics #Book #Movie #Traveling #Mac #iPhone

Nelder-Mead法(シンプレックス法)による非線形最適化MATLABサンプルプログラム

MATLAB

これなら分かる最適化数学―基礎原理から計算手法まで

これなら分かる最適化数学―基礎原理から計算手法まで

 

Nelder-Mead法

Nelder-Mead法は、非線形最適化法の一種です。

シンプレックス法やアメーバ法とも呼ばれます。


このNelder-Mead法は、多角形の探索領域を

広げたり、縮小したり、移動させることにより、

多次元非線形関数の最小値を探索します。


詳しいロジックに関しては、下記のwikipediaの記事と

末尾の参考文献を参照ください。

Nelder-Mead法の特徴


一般的な非線形最適化手法は、

目的関数の微分情報を使用しますが、

目的関数が不明な時や、

目的関数が複雑すぎて微分出来ない時などには使えません。


このNelder-Mead法はそのような微分情報が無くても

非線形最適化を実施することができます。


一方で、

微分情報を使わないため、

収束計算に時間がかかってしまう時がある

という問題があります。

MATLABサンプルプログラム

% -------------------------------------------------------------------------
%
% File : NelderMead.m
%
% Discription : Non-Linear optimizion with Nelder-Mead Method
%
% Reference: http://en.wikipedia.org/wiki/Nelder%E2%80%93Mead_method
% 
% Environment : Matlab
%
% Author : Atsushi Sakai
%
% Copyright (c): 2014 Atsushi Sakai
%
% License : Modified BSD Software License Agreement
% -------------------------------------------------------------------------

function [] = NelderMead()
close all;
clear all;

disp('NelderMead Simulation start!!')
%シミュレーション空間の最大最小範囲
maxxy=5;
minxy=-5;

%メッシュデータの生成
[x1range,x2range,yrange]=CreateMeshData(minxy,maxxy);

%初期パラメータの生成[x1 x2 y]
param=InitSearchParameter(minxy,maxxy);

iterationMax=50;%イタレーションの最大回数

%Parameter
row=1;%Reflection Parameter
kai=2;%Expansion Parameter
gamma=0.5;%Contraction Paremter 
sigma=0.5;%Shrink Parameter

for i=1:iterationMax
    %xmax change flag
    cf=0;
    
    %(1) Sort
    param=sortrows(param,3);
    
    %(2) Reflection
    xavg=mean(param(1:2,1:2));
    xr=xavg+row*(xavg-param(3,1:2));
    xr=[xr f(xr(1),xr(2))];
    
    %(3) Expand
    if xr(3)<param(1,3)
        xe=kai*xr(1:2)+(1-kai)*xavg;
        xe=[xe f(xe(1),xe(2))];
        if xe(3)<param(1,3)
            param(3,:)=xe;
            cf=1;
        elseif param(1,3)<xe(3)
            param(3,:)=xr;
            cf=1;
        end
    end
    %(4) contract
    if param(2,3)<xr(3)
        %(4-1) contract outside
        if param(2,3)<xr(3)&& xr(3)<param(3,3)
            xc=gamma*xr(1:2)+(1-gamma)*xavg;
            xc=[xc f(xc(1),xc(2))];
            if xc(3)<xr(3)
                param(3,:)=xc;
                cf=1;
            end
            %(4-2) contract inside
        elseif param(3,3)<xr(3)
            xcc=param(3,1:2)+(1-gamma)*xavg;
            xcc=[xcc f(xcc(1),xcc(2))];
            if xcc(3)<xr(3)
                param(3,:)=xcc;
                cf=1;
            end
        end
    end
    
    %(5) shrink
    if cf==0
        param(2,:)=sigma*(param(2,:)+param(1,:));
        param(3,:)=sigma*(param(3,:)+param(1,:));
    end
    
    %収束判定
    norm12=norm(param(1,:)-param(2,:));
    norm23=norm(param(2,:)-param(3,:));
    norm31=norm(param(3,:)-param(1,:));
    normSum=norm12+norm23+norm31;
    if normSum<0.001
        disp('Convergence!! final param [x1, x2, y] is');
        mean(param,1)
        break;
    end
    
    %シミュレーション結果の描画
    figure(1)
    hold off;
    contour(x1range,x2range,yrange,100); hold on;
    ag=[param;param];
    plot3(ag(:,1),ag(:,2),ag(:,3),'-k','markersize',10);
    plot3(ag(:,1),ag(:,2),ag(:,3),'.k','markersize',10);
    xlabel('x1');
    ylabel('x2');
    view(2)
    drawnow;
    pause(0.5);
     
end


function param=InitSearchParameter(minxy,maxxy)
%初期パラメータを作成する関数

x1=maxxy - (maxxy-minxy).*rand(1,3);
x2=maxxy - (maxxy-minxy).*rand(1,3);
y=f(x1,x2);
param=[x1' x2' y'];

function [x1range,x2range,yrange]=CreateMeshData(minxy,maxxy)
%シミュレーション空間のメッシュデータを作成する関数

%メッシュデータの作成
[x1range,x2range]=meshgrid(minxy:0.3:maxxy);
yrange=f(x1range,x2range);

function y = f(x1,x2)
% Himmelblau's function
% see Himmelblau's function - Wikipedia, the free encyclopedia 
% http://en.wikipedia.org/wiki/Himmelblau%27s_function
y=(x1.^2+x2-11).^2+(x1+x2.^2-7).^2;

MATLABのサンプルプログラムは

下記のGitHubリポジトリでも公開されています。

MATLABRobotics/NelderMead.m at master · AtsushiSakai/MATLABRobotics