%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Simple coordinate search algorithm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Kevin Passino
% Version: 4/6/00
%
% This program simulates the minimization of a simple function with
% the simple coordinate search method.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear % Initialize memory
p=2; % Dimension of the search space
Ncs=200; % Maximum number of iterations to produce
% Specify the pattern:
C=[eye(p,p) -eye(p,p) 0*ones(p,1) ];
[temp,sizeC]=size(C); % Specify the number of columns as the number of elements of C
% Next set the parameters of the algorithm:
gammac=1/2;
lambda=0*ones(1,Ncs+1);
lambda(1,1)=1;
% Set the parameter for the stopping criterion (based on small lambda has become)
epsilon=0.0001;
% Initialize to avoid use of memory manager:
J=0*ones(sizeC,1);
thetasbest=0*ones(p,1);
thetas=0*ones(p,sizeC-1);
P=0*ones(p,sizeC,Ncs+1); % Max possible used - stopping criteria may be invoked
% Next, pick the initial set
% With next initialization it falls into the local minimum at (20,15) and the
% stopping criterion is used (as set above)
P(:,1,1)=[16;
21];
% With next initialization it falls into the global minimum at (15,5) and the
% stopping criterion is used (as set above)
%P(:,1,1)=[20;
% 2];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Start the coordinate search loop
for j=1:Ncs
% Step 2: Compute cost at current estimate (note that second time around loop
% actually know this value, but recompute it - you could make this more efficient):
J(1,1)=optexampfunction([P(1,1,j);P(2,1,j)]);
% Step 3: Exploratory moves step
% Step 3 (a):
thetasbest=0*ones(p,1);
rho=0;
Jmin=J(1,1);
% Step 3 (b):
for i=1:sizeC-1
thetas(:,i)=lambda(1,j)*C(:,i); % Define perturbation for pattern point
P(:,i+1,j)=P(:,1,j)+thetas(:,i); % Compute the pattern point by perturbing from
% the from the center point
J(i+1,1)=optexampfunction([P(1,i+1,j);P(2,i+1,j)]); % Find cost at the pattern point
if J(i+1,1)0
P(:,1,j+1)=P(:,1,j)+thetasbest; % Found a better point so take it
lambda(1,j+1)=lambda(1,j); % and do not contract
else
P(:,1,j+1)=P(:,1,j); % Did not find a better point so use old one and contract
lambda(1,j+1)=gammac*lambda(1,j);
end
% Stopping criterion:
if lambda(1,j)