how to find minima of a 3d surface using Nelder-Mead?
12 views (last 30 days)
Show older comments
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
2 Comments
Star Strider
on 26 Oct 2014
You can actually save yourself a significant amount of work, unless you are supposed to program your own Nelder-Mead routine.
Accepted Answer
Star Strider
on 26 Oct 2014
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.
5 Comments
Star Strider
on 27 Oct 2014
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.
More Answers (0)
See Also
Categories
Find more on Systems of Nonlinear Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!