www.gusucode.com > mbcmodels 工具箱 matlab 源码程序 > mbcmodels/gamincon.m
function [X,FVAL,EXITFLAG,OUTPUT] = gamincon(FUN,X,LB,UB,A,B,NONLCON,options,DISPFUN,varargin)
%GAMINCON Finds the constrained minimum of a function of several variables.
% GAMINCON solves problems of the form:
% min f(x)
% x
% subject to:
% A*x <= b
% C(x) <= 0
% LB <= x <= UB
%
% Copyright 2000-2011 The MathWorks, Inc. and Ford Global Technologies, Inc.
defaultopt = struct('Display','none','LargeScale','on', ...
'TolX',1e-6,'TolFun',1e-6,'TolCon',1e-6,'DerivativeCheck','off',...
'Diagnostics','off','FunValCheck','off',...
'GradObj','off','GradConstr','off',...
'HessMult',[],...
'Hessian','off','HessPattern','sparse(ones(numberOfVariables))',...
'MaxFunEvals','100*numberOfVariables',...
'MaxSQPIter','10*max(numberOfVariables,numberOfInequalities+numberOfBounds)',...
'DiffMaxChange',1e-1,'DiffMinChange',1e-8,...
'PrecondBandWidth',0,'TypicalX','ones(numberOfVariables,1)',...
'MaxPCGIter','max(1,floor(numberOfVariables/2))', ...
'TolPCG',0.1,'MaxIter', 10,'OutputFcn',[],'PlotFcns',[],...
'RelLineSrchBnd',[],'RelLineSrchBndDuration',1,'NoStopIfFlatInfeas','off', ...
'PhaseOneTotalScaling','off', ...
'MaxGAIter',200,...
'MaxNoChange',30,...
'MaxSimilarity',0.95,...
'PopSize',20,...
'Hybridization',struct('name','sorted','option',[10,10]),...
'Fitness',struct('name','scaling','option',[]),...
'Selection',struct('name','normGeometric','option',0.08),...
'Crossover',struct('name',{'simple','arithmetic','heuristic'},...
'option',{1,1,[1 3]}),...
'Mutation',struct('name',{'boundary','multiNonUniform','nonUniform','uniform'},...
'option',{2,[3,3],[2,3],2}));
if nargin == 1 && nargout <= 1 && isequal(FUN,'defaults')
X = defaultopt;
return
end
if nargin < 4
error(message('mbc:gamincon:TooFewInputs'));
end
if nargin < 9
DISPFUN = [];
if nargin < 8
options = [];
if nargin < 7
NONLCON = [];
if nargin < 6
B = [];
if nargin < 5
A = [];
end
end
end
end
end
[msg,x,lb,ub] = checkbounds(X,LB,UB);
if ~isempty(msg)
error('mbc:gamincon:InvalidArgument', msg);
end
popSize = getgaopt(options,'PopSize',defaultopt,'fast');
numVars = length(x);
range = ub-lb;
X0 = [x range(:,ones(1,popSize-1)).*rand(numVars,popSize-1)+lb(:,ones(1,popSize-1))];
gradflag = strcmp(getgaopt(options,'GradObj',defaultopt,'fast'),'on');
gradconstflag = strcmp(getgaopt(options,'GradConstr',defaultopt,'fast'),'on');
if isempty(NONLCON)
constflag = false;
else
constflag = true;
end
if isempty(DISPFUN)
displflag = false;
else
displflag = true;
end
if ~isempty(FUN)
funfcn = optimfcnchk(FUN,'gamincon',length(varargin),gradflag);
else
error(message('mbc:gamincon:InvalidArgument1'));
end
if constflag
confcn = optimfcnchk(NONLCON,'gamincon',length(varargin),gradconstflag,1);
else
confcn{1} = '';
end
if (displflag)
disfcn = optimfcnchk(DISPFUN,'gamincon',length(varargin)+2);
else
disfcn{1} = '';
end
%----- Integer variables not excluded TO DO
x = checkOptimProblem(A,B,x,X0,numVars,funfcn,confcn,varargin{:});
global OPT_STOP OPT_STEP
OPT_STEP = 1;
OPT_STOP = 0;
%-----------------
suphybrid = {'sorted','random'};
supfitness = {'scaling','ranking'};
supselect = {'normGeometric','tournament','roulette'};
supxover = {'simple','arithmetic','heuristic'};
supmutat = {'boundary','multiNonUniform','nonUniform','uniform'};
%--------- hybridization, fitness and select must be scalars
opt = getgaopt(options,'Hybridization',defaultopt,'fast');
[hybrid,opthybrid] = opermatch(opt,suphybrid);
opt = getgaopt(options,'Fitness',defaultopt,'fast');
fitness = opermatch(opt,supfitness);
opt = getgaopt(options,'Selection',defaultopt,'fast');
[select,optselect] = opermatch(opt,supselect);
opt = getgaopt(options,'Crossover',defaultopt,'fast');
[xover,optxover] = opermatch(opt,supxover);
opt = getgaopt(options,'Mutation',defaultopt,'fast');
[mutat,optmutat] = opermatch(opt,supmutat);
nxover = length(xover);
nmutat = length(mutat);
L0 = zeros(1,popSize);
for k = 1:popSize
x(:) = X0(:,k);
L0(k) = feval(funfcn{3},x,varargin{:});
end
iter = 0;
switch disfcn{1}
case 'fun'
try
feval(disfcn{3},iter,X0,L0,varargin{:});
catch ME
error(message('mbc:gamincon:InvalidArgument13', 'GAMINCON cannot continue. User supplied display function', 'failed with the following error:', ME.message));
end
case ''
otherwise
error(message('mbc:gamincon:InvalidState2'));
end
F0 = zeros(size(L0));
L1 = zeros(size(L0));
X1 = zeros(size(X0));
c1 = zeros(numVars,1);
c2 = zeros(numVars,1);
p1 = zeros(numVars,1);
p2 = zeros(numVars,1);
bestx = zeros(numVars,1);
done = false;
tolfun = getgaopt(options,'TolFun',defaultopt,'fast');
maxgen = getgaopt(options,'MaxGAIter',defaultopt,'fast');
maxnoc = getgaopt(options,'MaxNoChange',defaultopt,'fast');
maxsim = getgaopt(options,'MaxSimilarity',defaultopt,'fast');
bufnoc = rand(maxnoc,1);
topnoc = 1;
iter = 1;
hciter = opthybrid{1}(1);
hcsize = opthybrid{1}(2);
hindex = zeros(1,hcsize);
switch hybrid
case 1
hcpind = zeros(1,popSize);
case 2
otherwise
end
switch fitness
case 1
case 2
rankfun = zeros(1,popSize);
rankfun(:) = 2*(0:popSize-1)/(popSize-1);
rankind = zeros(1,popSize);
otherwise
end
switch select
case 1
q = optselect{1};
r = q/(1-(1-q)^popSize);
indices = zeros(1,popSize);
indices(:) = popSize-1:-1:0;
gprob = r*(1-q).^indices;
fit = zeros(1,popSize);
urn = zeros(1,popSize);
case 2
trnSize = optselect{1};
tourind = zeros(popSize,trnSize);
indices = zeros(1,popSize);
case 3
rprob = zeros(1,popSize);
urn = zeros(1,popSize);
otherwise
end
for k = 1:nmutat
job = mutat(k);
mopt = optmutat{k};
switch (job)
case 1
case 2
MD = zeros(numVars,1);
CPOINT = false(numVars,1);
case 3
smshape = mopt(2);
case 4
otherwise
end
end
for k = 1:nxover
job = xover(k);
xopt = optxover{k};
switch (job)
case 1
case 2
case 3
retry = xopt(2);
t1 = zeros(numVars,1);
t2 = zeros(numVars,1);
otherwise
end
end
% Record number of calls to fmincon for output structure
NumberOfFminconCalls = 0;
while (~done)
if (rem(iter,hciter) == 0)
OPT_STEP = 1;
switch (hybrid)
case 1
[~,hcpind(:)] = sort(L0);
hindex(:) = hcpind(1:hcsize);
case 2
hindex(:) = floor(rand(1,hcsize)*popSize)+1;
otherwise
end
for k = hindex
x(:) = X0(:,k);
if isempty(options)
options = optimset('fmincon');
end
% Always want to use medium scale algorithm
fmcoptions = optimset(options,...
'largescale','off',...
'Algorithm','active-set');
[x,f,exitflag] = ...
fmincon(FUN,x,A,B,[],[],lb,ub,NONLCON, fmcoptions, varargin{:});
NumberOfFminconCalls = NumberOfFminconCalls + 1;
if exitflag >= 0
% Only use results if fmincon converges
X0(:,k) = x;
L0(k) = f;
end
if OPT_STOP
break
end
end
end
if displflag
feval(disfcn{3},iter,X0,L0,varargin{:});
end
OPT_STEP = 0;
if OPT_STOP
break
end
switch fitness
case 1
F0(:) = 1+max(L0)-L0;
case 2
[~,rankind(:)] = sort(-L0);
L1(:) = L0(rankind);
i1 = 1;
for i2 = [find(L1(1:popSize-1) ~= L1(2:popSize)),popSize]
F0(i1:i2) = sum(rankfun(i1:i2))/(i2-i1+1);
i1 = i2+1;
end
[~,rankind(:)] = sort(rankind);
F0(:) = F0(rankind);
otherwise
end
[~,ind] = max(F0);
bestx(:) = X0(:,ind);
bestl = L0(ind);
bufnoc(topnoc) = bestl;
topnoc = topnoc+1;
if topnoc > maxnoc
topnoc = 1;
end
if (iter > maxgen) || all(abs(bufnoc-bestl) <= tolfun) || ...
(sum(abs(L0-bestl) <= tolfun) >= maxsim*popSize)
done = true;
else
switch (select)
case 1
[~,indices(:)] = sort(F0);
fit(indices) = gprob;
fit(:) = cumsum(fit);
urn(:) = sort(rand(1,popSize));
i1 = 1;
i2 = 1;
while i2 <= popSize
if (urn(i2) < fit(i1))
X1(:,i2) = X0(:,i1);
L1(i2) = L0(i1);
i2 = i2+1;
else
i1 = i1+1;
end
end
case 2
tourind(:) = floor(rand(popSize,trnSize)*popSize)+1;
[~,indices(:)] = max(reshape(F0(tourind),trnSize,popSize));
indices(:) = diag(tourind(:,indices));
X1(:) = X0(:,indices);
L1(:) = L0(indices);
case 3
rprob(:) = cumsum(F0/sum(F0));
urn(:) = sort(rand(1,popSize));
i1 = 1;
i2 = 1;
while i2 <= popSize
if (urn(i2) < rprob(i1))
X1(:,i2) = X0(:,i1);
L1(i2) = L0(i1);
i2 = i2+1;
else
i1 = i1+1;
end
end
otherwise
end
for k = 1:nxover
job = xover(k);
xopt = optxover{k};
nx = xopt(1);
for l = 1:nx
i1 = round(rand*(popSize-1)+1);
i2 = round(rand*(popSize-1)+1);
p1(:) = X1(:,i1);
p2(:) = X1(:,i2);
switch job
case 1
cpoint = round(rand*(numVars-2)+1);
c1(:) = [p1(1:cpoint);p2(cpoint+1:end)];
c2(:) = [p2(1:cpoint);p1(cpoint+1:end)];
case 2
alfa = rand;
beta = 1-alfa;
c1(:) = alfa*p1+beta*p2;
c2(:) = beta*p1+alfa*p2;
case 3
feasible = false;
l1 = 0;
if (F0(i1) > F0(i2))
t2(:) = p1;
t1(:) = p2;
else
t2(:) = p2;
t1(:) = p1;
end
while l1 < retry
alfa = rand;
beta = 1+alfa;
c1(:) = beta*t2-alfa*t1;
if all(c1 <= ub) && all(c1 >= lb)
l1 = retry;
feasible = true;
else
l1 = l1+1;
end
end
if ~feasible
c1(:) = t1;
end
c2(:) = t2;
otherwise
end
if all(c1 == p1)
l1 = L1(i1);
elseif all(c1 == p2)
l1 = L1(i2);
else
l1 = feval(funfcn{3},c1,varargin{:});
end
if all(c2 == p1)
l2 = L1(i1);
elseif all(c2 == p2)
l2 = L1(i2);
else
l2 = feval(funfcn{3},c2,varargin{:});
end
X1(:,[i1,i2]) = [c1,c2];
L1([i1,i2]) = [l1;l2];
if OPT_STOP
break
end
end
if OPT_STOP
break
end
end
if ~OPT_STOP
for k = 1:nmutat
job = mutat(k);
mopt = optmutat{k};
nm = mopt(1);
for l = 1:nm
i1 = round(rand*(popSize-1)+1);
p1(:) = X1(:,i1);
c1(:) = p1;
switch job
case 1
cpoint = round(rand*(numVars-1)+1);
if (rand >= 0.5)
c1(cpoint) = ub(cpoint);
else
c1(cpoint) = lb(cpoint);
end
case 2
ratio = iter/maxgen;
alfa = (rand*(1-ratio)).^smshape;
beta = 1-alfa;
MD = rand(numVars,1);
CPOINT(:) = MD >= 0.5;
c1(CPOINT) = beta*p1(CPOINT)+alfa*ub(CPOINT);
CPOINT(:) = ~CPOINT;
c1(CPOINT) = beta*p1(CPOINT)+alfa*lb(CPOINT);
case 3
cpoint = round(rand*(numVars-1)+1);
ratio = iter/maxgen;
alfa = (rand*(1-ratio)).^smshape;
beta = 1-alfa;
if (rand >= 0.5)
c1(cpoint) = beta*p1(cpoint)+alfa*ub(cpoint);
else
c1(cpoint) = beta*p1(cpoint)+alfa*lb(cpoint);
end
case 4
cpoint = round(rand*(numVars-1)+1);
c1(cpoint) = lb(cpoint)+rand*range(cpoint);
otherwise
end
if all(c1 == p1)
l1 = L1(i1);
elseif all(c1 == p2)
l1 = L1(i2);
else
l1 = feval(funfcn{3},c1,varargin{:});
end
X1(:,i1) = c1;
L1(i1) = l1;
if OPT_STOP
break
end
end
if OPT_STOP
break
end
end
end
X0(:) = X1;
L0(:) = L1;
[~,ind] = max(L0);
X0(:,ind) = bestx;
L0(ind) = bestl;
iter = iter+1;
if OPT_STOP
break;
end
end
end
if ~OPT_STOP
OPT_STEP = 1;
% Run optimization one more time if the GA has still not converged.
x(:) = bestx;
if isempty(options)
options = optimset('fmincon');
end
% Always want to use medium scale algorithm
fmcoptions = optimset(options, ...
'largescale','off',...
'Algorithm','active-set',...
'MaxIter', 30*numVars);
[bestx,FVAL,EXITFLAG] = ...
fmincon(FUN,x,A,B,[],[],lb,ub,NONLCON, fmcoptions, varargin{:});
X(:) = bestx;
NumberOfFminconCalls = NumberOfFminconCalls + 1;
else
% GA has found a solution.
X(:) = bestx;
FVAL = bestl;
EXITFLAG = 1;
end
OUTPUT = struct('iterations',iter,...
'population',X0,...
'fitness',L0, ...
'numberoffminconcalls',NumberOfFminconCalls);
function [msg,xout,lbout,ubout] = checkbounds(x,lbin,ubin)
%CHECKBOUNDS Move the initial points within the (valid) bounds.
%
msg = '';
xout = x(:);
lbout = lbin(:);
ubout = ubin(:);
lenlb = length(lbout);
lenub = length(ubout);
nvars = length(x);
if (lenlb ~= nvars) || (lenub ~= nvars)
msg = 'Length of lower/upper bounds not equal to the number of variables.';
elseif any(~isfinite(lbout)) || any(~isfinite(ubout))
msg = 'Lower/upper bounds are not finite.';
elseif any(lbout >= ubout)
msg = 'Lower bounds exceed the corresponding upper bounds.';
else
if any(xout < lbout) || any(xout > ubout)
msg = 'Initial solution is not feasible.';
end
end
function [oper,opts] = opermatch(options,FIELDS)
%OPERMATCH Get operators and the associated parameters.
%
if ~isempty(options) && ~isstruct(options)
error(message('mbc:gamincon:InvalidState3'));
end
names = fieldnames(options);
if ~any( strncmp( 'name',names,4 ) ) || ~any( strncmp( 'option',names,6 ) )
error(message('mbc:gamincon:InvalidArgument14'));
end
NAMES = {options.name};
nop = length(NAMES);
oper = zeros(1,nop);
opts = cell(1,nop);
for k = 1:nop
name = NAMES{k};
ind = find( strncmpi( name,FIELDS,length(name) ) );
if isempty(ind)
error(message('mbc:gamincon:InvalidState4', NAMES{ k }));
elseif length(ind) > 1
ind = find( strcmpi( name,FIELDS ) );
if length(ind) ~= 1
error(message('mbc:gamincon:InvalidState5', NAMES{ k }));
end
end
oper(k) = ind;
opts{k} = options(k).option;
end
function x = checkOptimProblem(A,B,x,X0,numVars,funfcn,confcn,varargin)
x(:) = X0(:,1);
% Evaluate function
GRAD = zeros(numVars,1);
switch (funfcn{1})
case 'fun'
try
f = feval(funfcn{3},x,varargin{:});
catch ME
error(message('mbc:gamincon:InvalidArgument2', 'GAMINCON cannot continue. User supplied objective function', 'failed with the following error:', ME.message));
end
case 'fungrad'
try
[f,GRAD(:)] = feval(funfcn{3},x,varargin{:});
catch ME
error(message('mbc:gamincon:InvalidArgument3', 'GAMINCON cannot continue. User supplied objective/gradient function', 'failed with the following error:', ME.message));
end
case 'fun_then_grad'
try
f = feval(funfcn{3},x,varargin{:});
catch ME
error(message('mbc:gamincon:InvalidArgument4', 'GAMINCON cannot continue. User supplied objective function', 'failed with the following error:', ME.message));
end
try
GRAD(:) = feval(funfcn{4},x,varargin{:});
catch ME
error(message('mbc:gamincon:InvalidArgument5', 'GAMINCON cannot continue. User supplied objective gradient function', 'failed with the following error:', ME.message));
end
otherwise
error(message('mbc:gamincon:InvalidArgument6'));
end
% Evaluate constraints
switch (confcn{1})
case 'fun'
try
ctmp = feval(confcn{3},x,varargin{:});
c = ctmp(:);
cGRAD = zeros(numVars,length(c));
catch ME
error(message('mbc:gamincon:InvalidArgument7', 'GAMINCON cannot continue. User supplied nonlinear constraint function', 'failed with the following error:', ME.message));
end
case 'fungrad'
try
[ctmp,cGRAD] = feval(confcn{3},x,varargin{:});
c = ctmp(:);
catch ME
error(message('mbc:gamincon:InvalidArgument8', 'GAMINCON cannot continue. User supplied nonlinear constraint/gradient function', 'failed with the following error:', ME.message));
end
case 'fun_then_grad'
try
ctmp = feval(confcn{3},x,varargin{:});
c = ctmp(:);
catch ME
error(message('mbc:gamincon:InvalidArgument9', 'GAMINCON cannot continue. User supplied nonlinear constraint', 'function failed with the following error:', ME.message));
end
try
cGRAD = feval(confcn{4},x,varargin{:});
catch ME
error(message('mbc:gamincon:InvalidArgument10', 'GAMINCON cannot continue. User supplied nonlinear constraint', 'gradient function failed with the following error:', ME.message));
end
case ''
c = [];
cGRAD = zeros(numVars,length(c));
otherwise
error(message('mbc:gamincon:InvalidState'));
end
nonIneq = length(c);
[linIneq,Acol] = size(A);
[cgrow,cgcol]= size(cGRAD);
if ~isempty(A) && (Acol ~= numVars)
error(message('mbc:gamincon:InvalidArgument11'))
end
if ~isempty(A) && (linIneq ~= length(B))
error(message('mbc:gamincon:InvalidArgument12'));
end
if (cgrow ~= numVars) && (cgcol ~= nonIneq)
error(message('mbc:gamincon:InvalidState1'));
end