%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 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