MATLAB模拟退火算法求解器

%% Objective Function
clear, close all, clc
peaks

%% Nonlinear Constraint Function
type circularConstraint

%% Define Optimization Problem
problem = createOptimProblem('fmincon',...
                             'objective',@(x) peaks(x(1),x(2)), ...
                             'nonlcon',@circularConstraint,...
                             'x0',[-1 -1],...
                             'lb',[-3 -3],...
                             'ub',[3 3],...
                             'options',optimset('OutputFcn',...
                                                @peaksPlotIterates))
                             
%% Run the solver |fmincon| from the inital point
[x,f] = fmincon(problem)                                                    

%% Use Simmulated Annealing to Find the Global Minimum

problem.solver  = 'simulannealbnd';
problem.objective = @(x) peaks(x(1),x(2)) + (x(1)^2 + x(2)^2 - 9);
problem.options = saoptimset('OutputFcn',@peaksPlotIterates,...
                             'Display','iter',...
                             'InitialTemperature',10,...
                             'MaxIter',300)

[x,f] = simulannealbnd(problem) 

function varargout = peaksPlotIterates(varargin)
% Output function that plots the iterates of the optimization algorithm.

%  Copyright (c) 2010, The MathWorks, Inc.
%  All rights reserved.

% Check if caller is from global or optimization toolbox
optimValues = varargin{
   2};
state = varargin{
   3};
if nargout > 1
    if isfield(optimValues,'x') % simulated annealing options,optimvalues,flag
        x = optimValues.x;
        varargout{
   1} = false;
        varargout{
   2} = x; % options field
        varargout{
   3} = false;
    else % pattern search optimvalues,options,flag,interval
        optimValues = varargin{
   1};
        state = varargin{
   3};
        if isfield(optimValues,'x')
            x = optimValues.x;
            varargout{
   1} = false;
            varargout{
   2} = x; % options field
            varargout{
   3} = false;
        else % gentic algorithm options,state,flag,interval
            x = varargin{
   2}.Population;
            optimValues.iteration = -1;
            varargout{
   1} = varargin{
   2};
            varargout{
   2} = varargin{
   1};
            varargout{
   3} = false;
        end
    end
else
    x = varargin{
   1};
    varargout{
   1} = false;
end

% Check for state
switch state
    case 'init'
        % Plot objective function surface
        PlotSurface(x,peaks(x(:,1),x(:,2)));
    case 'iter'
        if ~(optimValues.iteration == 0)
            % Update surface plot to show current solution
            PlotUpdate(x,peaks(x(:,1),x(:,2)));
        end
    case 'done'
        if ~(optimValues.iteration == 0)
            % After optimization, display solution in plot title
            DisplayTitle(x,peaks(x(:,1),x(:,2)))
        end
end

% -------------------------------------------------------------------------
% helper function PlotSurface
% -------------------------------------------------------------------------
function PlotSurface(x,z,varargin)

% Check to see if figure exists, if not create it
h = findobj('Tag','PlotIterates');
if isempty(h)
    h = figure('Tag','PlotIterates','Name','Plot of Iterates', ...
        'NumberTitle','off');

% Plot the objective function
[X,Y,Z] = peaks(100);
zlower = -15;
axis([-3 3 -3 3 zlower 10]);
hold on
surfc(gca,X,Y,Z,'EdgeColor','None','FaceColor','interp')
xlabel('X'), ylabel('Y'), zlabel('Z')
view([-45 30])
shading interp
lightangle(-45,30)
set(findobj(gca,'type','surface'),...
    'FaceLighting','phong',...
    'AmbientStrength',.3,'DiffuseStrength',.8,...
    'SpecularStrength',.9,'SpecularExponent',25,...
    'BackFaceLighting','unlit');

% Plot constraint on lower contour plot
hc=0; k=0; r=3; N=256; % circle parameters
t=(0:N)*2*pi/N;
xc = r*cos(t)+hc;
yc = r*sin(t)+k;

% bounds
ax = axis;%.*[1.1 1.1 1.1 1.1 1 1]; 
xbound = ( ax(1):(ax(2)-ax(1))/N*4:ax(2) )';
ybound = ( ax(3):(ax(4)-ax(3))/N*4:ax(4) )';
len = length(xbound);
xbox = [xbound;  xbound(end)*ones(len-1,1);
        xbound(end-1:-1:1); xbound(1)*ones(len-2,1)]; 
ybox = [ybound(1)*ones(len,1); ybound(2:end);
        ybound(end)*ones(len-1,1); ybound(end-1:-1:2)];

boxCon = [(1:length(xbox)-1)' (2:length(ybox))'; length(xbox) 1];
cirCon = [(1:length(xc)-1)' (2:length(yc))'; length(x) 1] + length(xbox);

warning off
DT = DelaunayTri([[xbox(:); xc(:)] [ybox(:); yc(:)]], [boxCon; cirCon]);
warning on
inside = inOutStatus(DT);
cx = caxis;
trisurf(DT(inside,:),DT.X(:,1),DT.X(:,2),...
       zlower*ones(size(DT.X(:,1))),'EdgeColor','none',...
       'FaceColor',[0.9 0.9 0.9]);
caxis(cx)
hold off

% colors to use for multiple staring points
ms.index = 1;
ms.Colors = ['rgbcmyk'];
set(h,'UserData',ms);

end

PlotUpdate(x,z)
if nargin > 2
    DisplayTitle(x,z,varargin{
   1})
else
    DisplayTitle(x,z,'Initial')
end

% -------------------------------------------------------------------------
% helper function PlotUpdate
% -------------------------------------------------------------------------
function PlotUpdate(x,z)

% Check to see if figure exists, if not, create
h = findobj('Tag','PlotIterates');
if isempty(h)
    PlotSurface(x,z,'Current')
    h = gcf;
end
% Update Plot with New Point
figure(h)
ms = get(h