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