www.gusucode.com > mbcdata 工具箱 matlab 源码程序 > mbcdata/mbcgacreationnonlinearfeasible.m

    function Population = mbcgacreationnonlinearfeasible(usegrad,GenomeLength,~,options,ConstrFcn)
%MBCGACREATIONNONLINEARFEASIBLE Creates the initial population for GA.
%   POP = GACREATIONNONLINEARFEASIBLE(NVARS,FITNESSFCN,OPTIONS,CONSTRFCN) 
%   creates an initial population that GA will then evolve into a solution.
%
%   Example:
%     options = gaoptimset('CreationFcn',@gacreationnonlinearfeasible);

%   Copyright 2015 The MathWorks, Inc.

totalPopulation = sum(options.PopulationSize);
initPopProvided = size(options.InitialPopulation,1);
individualsToCreate = totalPopulation - initPopProvided;

A = []; b = []; Aeq = []; beq = []; lb = []; ub = []; 
if isfield(options,'LinearConstr')
    % Extract linear constraint data and append columns for the slack
    % variable of the minimum infeasibility formulation.
    b = options.LinearConstr.bineq;
    beq = options.LinearConstr.beq;
    A = [options.LinearConstr.Aineq  zeros(numel(b),1)];
    Aeq = [options.LinearConstr.Aeq  zeros(numel(beq),1)];
    lb = options.LinearConstr.lb(:);
    ub = options.LinearConstr.ub(:);
end

% Initialize Population to be created
initPts = [];
if initPopProvided > 0
    % Use initial population provided already
    initPts = options.InitialPopulation;
end

if individualsToCreate > 0
   % Create the rest of the initial points from gacreationuniform
   % NOTE: overwrite (or create) the "type" field as 'boundconstraints' to
   % ensure points will be generated.
   
   % TODO: should we call gacreationlinearfeasible if there are linear
   % constraints?
   options.LinearConstr.type = 'boundconstraints';
   initPts = [initPts;
              gacreationuniform(GenomeLength,[],options)];
end

if usegrad
    opts = optimoptions('fmincon','Display','off','Algorithm','interior-point', ...
                    'TolCon',options.TolCon,'GradObj','on','GradConstr','on');
else
    opts = optimoptions('fmincon','Display','off','Algorithm','interior-point', ...
                    'TolCon',options.TolCon,'GradObj','on');
end

% TODO: set MaxTime to a fraction of GA's option MaxTime or some
% pre-determined limit?
if isfinite(options.TimeLimit)    
    opts.OutputFcn = @(a,b,c) maxTimeManager(a,b,c,state.StartTime,options.TimeLimit);
end

% Solve a set of feasibility problems using fmincon
Population = Inf(totalPopulation,GenomeLength);

% choose whether to build parallel or serial
parallelAvailable = mbcfoundation.DCTManager.isDCTAvailable();


if parallelAvailable
    pp = gcp;
    % store ConstrFcn on workers
    fAll = parfevalOnAll(pp,@storeConstrFcn,0,ConstrFcn);
    wait(fAll);
    % distribute calls to callSolver
    for k = 1:totalPopulation
        f(k) = parfeval(pp, @callSolver, 1, GenomeLength,...
                initPts(k,:), A,b,Aeq,beq,lb,ub,[],opts); %#ok<AGROW>
    end
    % collect results from distributed runs
    for k=1:totalPopulation
        [~,value] = fetchNext(f);
        Population(k,:) = value;
        % check if stop button had been pressed
        stop = getStopState(cgoptimstore);
        if stop
            cancel(f);
            break
        end
    end
    % clear stored ConstrFcn on workers
    fAll = parfevalOnAll(pp,@storeConstrFcn,0,[]);
    wait(fAll);

else
    % call callSolver in serial
    for k = 1:totalPopulation
        Population(k,:) = callSolver(GenomeLength,initPts(k,:), ...
                            A,b,Aeq,beq,lb,ub,ConstrFcn,opts);
    end
end
   
%--------------------------------------------------------------------------
function x = callSolver(GenomeLength,x0,A,b,Aeq,beq,lb,ub,nonlcon,options)
% NOTE: It's expected that this function always returns row vectors.


maxInfeasObj = @(x)0;
zeroGrad = zeros(1,GenomeLength); 
maxInfeasGrad = @(x)zeroGrad;
obj = {maxInfeasObj,maxInfeasGrad};
if isempty(nonlcon)
    % retrieve constraint function handle from persistent store
    nonlcon = storeConstrFcn;
end

if isFeasible(nonlcon,x0) && iSatisfyLinearConstraint(x0, A, b)
    x = x0';
else
    try
        [x,~,flag] = fmincon(obj,x0(:),A,b,Aeq,beq,lb,ub,nonlcon,options);
        x = x'; % Make column vector
    catch
        % fmincon errored, we can only return the initial point
        x = x0';
        flag = -1;
    end
    if flag < 0
        % fmincon didn't achieve feasibility. We should return the result if
        % possible. In all likelihood, fmincon should at least maintain the
        % same level of infeasibility as x0, if not improve.
        if isempty(x)
            x = x0';
        end
    end
end

function out = iSatisfyLinearConstraint(x, A, b)
if isempty(A)
    out = true;
else
    out = all(A*x<b);
end

%--------------------------------------------------------------------------
function stop = maxTimeManager(~,~,~,startTime,timeLimit)
stop = cputime - startTime >= timeLimit;

function isfeasable = isFeasible(ConstrFcn,x)
evalInit = ConstrFcn(x);
isfeasable = all(evalInit<eps);


function ConstrFcn = storeConstrFcn(ConstrFcn)
%storeConstrFcn store constraint function handle in a persistent variable on workers
persistent STORE_CONSTRFCN

if nargin==1
    STORE_CONSTRFCN = ConstrFcn;
else
    ConstrFcn = STORE_CONSTRFCN;
end