MATLAB: How to find minima of a 3d surface using Nelder-Mead

downhill simplexerrorglobal minimanelder-meadobjective-functionresponse surface function

I have function to generate a response surface function of the RMSError of a model. I need to then use the Nelder-Mead (Downhill Simplex) method to find the global minima. I can't figure out how to connect the two pieces (the response surface function and the Nelder-Mead code).

My first problem is calculating the centroid in line 27 of the Nelder-Mead code.

Any help would be appreciated, especially with explicit suggestions. I see what I need to do but not how.

Code is below, data is attached.

Code for RMSE Response Surface Function:

function [err_rmse] = RMSERespSurf(TimeDay, Precipmmday)
SynthK=0.1;
SynthSmax=20;
Baseflow = zeros(size(TimeDay));
Storage = zeros(size(TimeDay));
for nt = 2:length(Precipmmday)
    Baseflow (nt) = (Storage(nt)+Precipmmday(nt))*SynthK;
    Storage(nt) = (Storage(nt-1)+Precipmmday(nt-1))*(1-SynthK);
end
SynthStorage = Storage;
% SynthStorage (:,1) =[];
SynthSpill=SynthStorage - SynthSmax;
SynthSpill(SynthSpill<0)=0;
SynthStorage(SynthStorage>SynthSmax)=SynthSmax;
SynthBaseflow = Baseflow;
SynthOutflow = SynthBaseflow + SynthSpill;
% temp model variables
ModelBaseflow = zeros(size(TimeDay));
ModelStorage = zeros(size(TimeDay));
ModelSpill = zeros(size(TimeDay));
% grid is the resolution of parameter space (grid x grid)
grid = 100;
% preallocate space for error calcs
err_rmse = zeros(grid);
% define parameters
Smax = linspace(10,100,grid);
K = linspace(0.05,0.5,grid);
% loop for Smax
for ns = 1:length(Smax)
    %     loop for K
    for nk = 1:length(K)
        %         loop for timeseries (k1)
        for nt = 2:length(TimeDay)
            % calculate recursive time series
            ModelBaseflow(nt)  = (ModelStorage(nt-1)+Precipmmday(nt-1))*K(nk);
            ModelStorage(nt) = (ModelStorage(nt-1)+Precipmmday(nt-1))*(1-K(nk));
        end
          % calculate secondary variables (not recursive, but depend on
          % baseflow and storage)
          ModelSpill = ModelStorage - Smax(ns);
          ModelSpill(ModelSpill<0)=0;
          ModelStorage(ModelStorage>Smax(ns)) = Smax(ns);
          ModelOutflow  = ModelBaseflow + ModelSpill;
          % calculate errors 
          Difference = ModelOutflow-SynthOutflow;
          err_rmse(ns,nk) = rms(Difference);
      end
end  
surf(Smax, K, err_rmse);
ylabel('k')
xlabel('Smax')
zlabel('err_rmse')

Code for Nelder-Mead:

function [Minima]=Nelder_Mead_MultiDimMatrixForm(RMSERespSurf)
clc;
StdVal=10; %any value for convergence
n=3; %value of N+1
P=1; %reflection coefficient
Chi=2; %expansion coefficient
Gamma=0.5; %contraction coefficient
Sig=0.5; %shrink coefficient
SortVertices = CreateInitialSimplex(n);
tic;
i=1;
minmas=0;
while(StdVal >= 0.00001)
      SortVertices = BubbleSort(SortVertices,n);
      Centroid = CalculateCentroid(SortVertices,n);
      StdVal = CalculateStd(SortVertices,n);
      minmas(i)=SortVertices(1).value; i=i+1;
  %     Simplex_2D(SortVertices); drawnow; 
      Reflect.coord = (1+P).*Centroid.coord - P.*SortVertices(n).coord; %Reflect
      Reflect.value = RMSERespSurf(Reflect);
      if(SortVertices(1).value <= Reflect.value && Reflect.value < SortVertices(n-1).value)
          SortVertices(n)=Reflect;
          continue; %if the above criterion is sattisfied, then terminate the iteration
      end
      if(Reflect.value < SortVertices(1).value) %Expand
          Expand.coord = (1-Chi).*Centroid.coord+Chi.*Reflect.coord;
          Expand.value = RMSERespSurf(Expand);
          if(Expand.value < Reflect.value)
              SortVertices(n) = Expand;
              continue;
          else
              SortVertices(n) = Reflect;
              continue;
          end 
      end
      if(SortVertices(n-1).value <= Reflect.value) %Contract
              ContractOut.coord = (1-Gamma).*Centroid.coord + Gamma.*Reflect.coord; %Contract Outside
              ContractOut.value = RMSERespSurf(ContractOut);
              ContractIn.coord = (1-Gamma).*Centroid.coord + Gamma.*SortVertices(n).coord;  %Contract Inside

              ContractIn.value= RMSERespSurf(ContractIn);
              if(SortVertices(n-1).value <= Reflect.value && Reflect.value < SortVertices(n).value && ContractOut.value <= Reflect.value)
                  SortVertices(n) = ContractOut;
                  continue;
              elseif(SortVertices(n).value <= Reflect.value && ContractIn.value < SortVertices(n).value) %Contract Inside
                  SortVertices(n) = ContractIn;
                  continue;
              else
                  for i=2:n
                      SortVertices(i).coord = (1-Sig).*SortVertices(1).coord + Sig.*SortVertices(i).coord;
                      SortVertices(i).value = RMSERespSurf(SortVertices(i));   
                  end
              end
      end    
  end
  toc
  Minima=SortVertices(1);
  % plot(minmas);
  % a=text(50,250, strcat('No: of iterations =  ', num2str(i)));
  % set(a, 'FontName', 'Arial', 'FontWeight', 'bold', 'FontSize', 18, 'Color', 'red');
  % xlabel('iterations');
  % ylabel('error');
  % title('Error Vs Iterations');
  end
function [value]= RMSERespSurf(V) %Write your function in matrix form
%in case of any of the three trial functions un comment two lines   
    x=V.coord(1); y=V.coord(2);  %Rosenbrock's Function
    value=(1-x)^2+100*(y-x^2)^2;
%    x1=V.coord(1);x2=V.coord(2);x3=V.coord(3);x4=V.coord(4); %Powell's Quadratic Function
%    value=(x1+10*x2)^2+5*(x3-x4)^2+(x2-2*x3)^4+10*(x1-x4)^4;
%    x1=V.coord(1); x2=V.coord(2); x3=V.coord(3); %Fletcher and Powell's Helical Valley Function
%     theta=atan2(x1,x2);
%     value=100*(x3-10*theta)^2+(sqrt(x1^2+x2^2)-1)^2+x3^2;
%   value=([1,-1;0,1]*V.coord-[4;1])'*V.coord; %f(x,y)=x^2-4*x+y^2-y-x*y;
%   value=([1,0,-1;0,1,0;1,0,1]*V.coord-[3;10;-7])'*V.coord; %f(x,y,z)=x^2+y^2+z^2-3*x-10*y+7
end
function [Vertices]=CreateInitialSimplex(n)
      ExpectMin=rand((n-1),1).*50;
      Vertices(1).coord=ExpectMin; % expected minima
      Vertices(1).value=RMSERespSurf(Vertices(1));
      for i=2:n
          Vertices(i).coord=ExpectMin+50.*rand((n-1),1); %100 is the scale factor
          Vertices(i).value=RMSERespSurf(Vertices(i));
      end
end
function [SortVertices]=BubbleSort(Vertices,n)
      while(n~=0)
          for i=1:n-1
              if(Vertices(i).value<=Vertices(i+1).value)
                  continue;
              else
                  temp=Vertices(i);
                  Vertices(i)=Vertices(i+1);
                  Vertices(i+1)=temp;
              end
          end
          n=n-1;          
      end
      SortVertices=Vertices;
end
function [Centroid]=CalculateCentroid(Vertices,n)
      Sum=zeros((n-1),1);
      for i=1:n-1
          Sum=Sum+Vertices(i).coord;
      end
      Centroid.coord=Sum./(n-1);
      Centroid.value=f(Centroid);
end
function[StdVal]=CalculateStd(Vertices,n) % this is the tolerance value, the standard deviation of the converging values
      for i=1:n
          ValueArray(i)=Vertices(i).value;
      end
      StdVal=std(ValueArray,1);
end

Best Answer

I am not quite sure that I understand what you’re doing, but if you have a model of your function and the data you are fitting, you can do something similar to a nonlinear curve fit. So (remembering this from our earlier conversations) you are fitting a model to data, you might do something like this:

Precipmmday = [(1xN) vector];                              % Input Data (Known)
Streamflow  = [(1xN) vector];                              % Output Data (Guess)
% K = b(1);  Smax = b(2);                                  % Parameter Vector For Model
Model = @(b,x) b(1).*x + log(b(2).*x);                     % Fictitious Model
Cost = @(b) sum((Model(b,Precipmmday) - Streamflow).^2);   % Sum-Of-Squares Cost Function
b0 = realistic estimate of [K, Smax];                      % Initial Parameter Estimates
b = fminsearch(Cost, b0);                                  % Minimise ‘Cost(b)’

If you are optimising with respect to other (different or additional) parameters, I need to know what they are, and the model you are fitting. Unfortunately, this has managed to escape me thus far. I’ve been primarily concerned with helping you get your code running.

I’ll help as much as I can.

Best Answer

Related Solutions

MATLAB: How to loop through a function twice

MATLAB: How to show variation of two variables

Related Question