how to find minima of a 3d surface using Nelder-Mead?

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

[Star Strider]

You can actually save yourself a significant amount of work, unless you are supposed to program your own Nelder-Mead routine.

See the documentation for the fminsearch Algorithm.

[Vert]

That looks perfect! but I haven't been able to figure out how to connect it to my code for the RMSE response surface.

Maybe I'm missing something obvious but I don't see how to tell it to locate the minima within the RMSE, using three random initial coordinate points of K and Smax, and return the coordinates of the minima. I've read through the info and tried a few ways but I'm not getting it.

Answer 1

[Star Strider]

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.

[Vert]

Thank you for your help! Yes, I am optimizing with respect to just two parameters (K and Smax). I am optimizing using root mean square error as the objective function.

I tried doing what you said (see below) but it says, "Error in fminsearch (line 191) fv(:,1) = funfcn(x,varargin{:}); ". Do you know what that means?

Also why does the sum of squares cost need to be there?

Model = @ RMSERespSurf;                    
  Cost = @(b) sum((Model(b,Precipmmday) - Streamflow).^2);   % Sum-Of-Squares Cost Function
  f = 10;
  g = 100;
  TestSmax = (g-f).*rand + f;
  h = 0.01;
  j = 1;
  TestK = (j-h).*rand + h;
  b0 = [TestK, TestSmax];                      % Initial Parameter Estimates
  b = fminsearch(Cost, b0);

[Star Strider]

My pleasure!

I’m not certain what the fminsearch error refers to. Something in your code gave it indigestion, and it complained.

I would like to see the code for your model. It needs to be of the form my prototype model function is, in that it takes a parameter vector (I call it b out of convention), and your independent variable (in your application, Precipmmday), and returns a vector of estimates for you dependent variable corresponding to every element of your independent variable.

The cost function definitely needs to be there. The fminsearch function is an optimisation function, not a curve-fitting function, so you have to provide the cost function to minimise the sum-squared error. It essentially calculates the RMSE (or a version of it), so don’t calculate it in your objective function (the function you want to fit).

If you have a description of the function you want to fit, post it or attach it here. I may not understand it if it is heavily technical with surface hydrology terms, but so long as you’re here to translate it, this shouldn’t be too difficult.

[Vert]

I think my problem is that I'm missing something basic about how to write a function or model in MatLab and then refer to it. Below is what I'm using as my function and I was trying to call it as "@RMSERespSurf" above.

I want to be able to refer to the surface of error output by the function below. I'm trying to refer to (K, Smax) coordinate points.

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')

[Star Strider]

We need to discuss what you are doing because I don’t understand what you want to optimise. If you want to minimise the RMSE, that’s easy enough to do in the ‘Cost’ function directly if we’re doing curve-fitting.

You might be able to optimise your ‘RMSERespSurf’ function directly without using ‘Cost’. I don’t understand enough just now to be able to determine that.

The point is that in an optimisation problem, you are optimising your model with respect to certain parameters to achieve a desired goal. I assume from our previous discussions that you want to optimise ‘K’ and ‘Smax’ (correct me if I’m wrong) but I have no idea how you are constructing your model. The RMSE has to be calculated by comparing your model to some real-world data. (BTW, I’ve figured out that ‘K’ is a dimensionless decimal fraction on [0,1], but I still have no idea what ‘Smax’ is. Please enlighten me.)

Please describe your model in simple terms (not code) so that I can understand it and what you are doing. What real world data are you comparing it to? Include references to ‘K’ and ‘Smax’ in your description so I know where they fit in.

[Vert]

sorry for not being clear! I think I've been staring at this too long. I am looking at a tank.

Rain falls into the tank (the amount in mm per day is Precipmmday). The tank has a maximum storage capacity of Smax. At the base of the tank is a pipe. Baseflow comes out of the pipe. The amount of baseflow is determined by the coefficient K. If there is too much rain coming in, it exceeds the maximum storage capacity and K is too small to allow it to leave through the pipe, it spills over the side of the tank. Outflow is the sum of baseflow and spillage.

The data I have is for the daily rainfall. Further, I am choosing values for K and Smax and using these values as observed reality. I'm calling these the "synthetic observations". I'm then calculating the amount of outflow given the synthetic observation values for K and Smax. The next step is to compare this to the parameter space defined by the range of possible K and Smax values. Finally, I am calculating the RMSE for the parameter space.

With that done, I am now going back and finding the optimum K and Smax values using Nelder-Mead. The global error minima should get back to the original K and Smax values that I used as synthetic observations.