%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Torczon's multidirectional search algorithm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Kevin Passino
% Version: 3/21/00
%
% This program simulates the minimization of a simple function with
% the multidirectional search method.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear % Initialize memory
p=2; % Dimension of the search space
Nmds=200; % Maximum number of iterations to produce
% Next set the parameters of the algorithm:
gammae=2;
gammac=1/2;
% Set the parameter for the stopping criterion (based on how much the minimum
% cost vertex has changed from one iteration to the next, after one iteration)
epsilon=0.0000000001;
stopflag=0; % Used to signal termination criterion has been met
% Initialize to avoid use of memory manager:
J=0*ones(p+1,1);
Jrot=0*ones(p,1);
Jexp=0*ones(p,1);
Jcont=0*ones(p,1);
thetarot=0*ones(p,p);
thetaexp=0*ones(p,p);
thetacont=0*ones(p,p);
P=0*ones(p,p+1,Nmds+1); % Max possible used - stopping criteria may be invoked
% Next, pick the initial simplex
% With next initialization it falls into the local minimum at (20,15) and the
% stopping criterion is used (as set above)
P(:,:,1)=[15.2 14.9 16;
20.6 20.1 19.3];
% With next initialization it falls into the upper left corner local minimum and the
% stopping criterion is used (as set above)
%P(:,:,1)=[16 14 16;
% 21 22 19];
% With next initialization it falls into the global minimum at (15,5) and the
% stopping criterion is used (as set above)
%P(:,:,1)=[16 15 14;
% 4 7 3];
% With next initialization it falls into the global minimum at (15,5) and the
% stopping criterion is used (as set above)
%P(:,:,1)=[30 0 15;
% 0 0 30];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Start the MDS loop
for j=1:Nmds
if stopflag==1
break
end
% Step 1: find min cost index and place it at P(:,1,j) (first column)
for i=1:p+1 % This is not the simplest way computationally - for j>1
% this approach results in recomputation of costs that
% were already computed below in rotation, expansion, or contraction
% Just do it this way for programming simplicity (laziness on the author's part)
J(i,1)=optexampfunction([P(1,i,j);P(2,i,j)]);
end
[Jmin,istar]=min(J(:,1)'); % Find the index istar of the minimum cost
firstcolumn=P(:,1,j); % Swap the first column (vertex 0) with the best vertex
P(:,1,j)=P(:,istar,j);
P(:,istar,j)=firstcolumn;
minvertex_notreplaced=1; % =1 represents that the minimum vertex has not been
% replaced yet (and steps 2,3,4 continue till it is)
while minvertex_notreplaced==1
% Need stopping criterion here (based on change in size of best vertex)
if j>1
% Compute min vertex change amount
change=(P(:,1,j)-P(:,1,j-1))'*(P(:,1,j)-P(:,1,j-1));
if change