www.gusucode.com > mbcdata 工具箱 matlab 源码程序 > mbcdata/mbcOSfmincon.m
function out = mbcOSfmincon(action, in) %MBCOSFMINCON CAGE Optimization library function % % MBCOSFMINCON contains the configuration and evaluation routines for the % CAGE Optimization, FOPTCON. This function is called by CAGE using the % syntaxes below. % % OPTOBJ = MBCOSFMINCON('OPTIONS', OPTOBJ) is called by CAGE when an % FOPTCON optimization is created in CAGE. This call defines the % structure of the optimization problem, using a CGOPTIMOPTIONS object. % The completed configuration for the Optimization is returned to CAGE % via OPTOBJ. % % OPTIMSTORE = MBCOSFMINCON('EVALUATE', OPTIMSTORE) is called by CAGE % when a user runs an FOPTCON optimization. Optimization parameters are % passed to the evaluation routine through OPTIMSTORE. The results of the % optimization are stored in the OPTIMSTORE object which is returned to % CAGE for further processing. % % See also MBCOSNBI % Copyright 2000-2016 The MathWorks, Inc. and Ford Global Technologies, Inc. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Deal with the two action cases %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% switch lower(action) case 'options' % Return the updated options object out = i_Options(in); case 'evaluate' % Return a handle to the main evaluation routine out = i_Evaluate(in); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % options Function %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function opt = i_Options(opt) % Set the evaluation function interface to always be the latest one opt = setRunInterfaceVersion(opt, 0); % Add a name opt = setName(opt, 'foptcon'); % Add a description opt = setDescription(opt, 'Single objective optimization subject to constraints'); % Set up the free variables %%%% Any number are allowed %%%% % Set up the objective functions %%%% Only one objective is allowed %%%% % Set up the constraints %%%% Any number of constraints are allowed %%%% % Allow objectives and constraints to be renamed opt = setRenameMode(opt, true); % allow equality constraints opt = setEqConMode(opt,'on'); % Set up the operating point sets %%%% Any number of operating point sets are allowed %%%% opt = setOperatingPointsMode(opt, 'any'); % Set up the optimization parameters opt = addParameter(opt, 'Display', {'list', {'none', 'final', 'iter'}}, 'none'); opt = addParameter(opt, 'MaxIter', {'integer', 'positive'}, 500, 'Maximum iterations'); opt = addParameter(opt, 'MaxFunEvals', {'integer', 'positive'}, 2000, 'Maximum function evaluations'); opt = addParameter(opt, 'TolX', {'number', 'positive'}, 1e-6, 'Variable tolerance'); opt = addParameter(opt, 'TolFun', {'number', 'positive'}, 1e-6, 'Function tolerance'); opt = addParameter(opt, 'TolCon', {'number', 'positive'}, 1e-6, 'Constraint tolerance'); opt = addParameter(opt, 'Algorithm', ... {'list', {'active-set','sqp', 'interior-point'}}, 'interior-point',... 'Constrained optimization algorithm'); opt = addParameter(opt, 'DiffMinChange', {'number', 'positive'}, 1e-8, 'Minimum change in variables for gradient'); opt = addParameter(opt, 'DiffMaxChange', {'number', 'positive'}, 0.1, 'Maximum change in variables for gradient'); opt = addParameter(opt, 'NumStart', {'integer', [1 10000]}, 1, 'Number of start points'); opt = addParameter(opt, 'RunFromFeasOnly', 'boolean', false, 'Run from feasible start points only:'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Evaluation Function %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function optimstore = i_Evaluate(optimstore) % Optimization routine must conform to the following API % optimstore = <OPTIMIZATION_FUNCTION>(optimstore) % % Inputs: optimstore: store of all the CAGE data needed % Outputs: optimstore: updated object containing all the CAGE data, including outputs % % Get bounds and linear constraints A = getA(optimstore); B = getB(optimstore); lb = getLB(optimstore); ub = getUB(optimstore); % Get the start point options. We have to support existing sessions with a % different parameter structure. Test for the existence of an old parameter % structure in a try-catch. try NumStart = getParam(optimstore, 'NumStart'); RunFromFeasOnly = getParam(optimstore, 'RunFromFeasOnly'); catch %#ok<*CTCH> % Old structure - does not support NumStart and RunFromFeasOnly NumStart = 1; RunFromFeasOnly = false; end % Generate the initial conditions x0 = i_generateStartPoints(optimstore, NumStart, RunFromFeasOnly); if ~isempty(x0) % Get the number of start positions NumStart = size(x0, 1); % Set some extra options in fmincon GRADREQ = ~isScalarFreeVariables(optimstore); if GRADREQ gradstr = 'on'; else gradstr = 'off'; end options = optimset(optimset('fmincon'),... 'LargeScale', 'off', ... 'Algorithm', 'active-set', ... 'OutputFcn',@i_Stop,... 'GradObj', gradstr, ... 'GradCon', gradstr); % Get the user defined options and put them in an optimization toolbox % options structure for fmincon options = optimset(options,optimstore); % Initialise return quantities bestx = x0; bestfun = Inf*ones(NumStart, 1); exitFlag = zeros(NumStart, 1); OUTPUT(NumStart).iterations = []; OUTPUT(NumStart).funcCount = []; OUTPUT(NumStart).stepsize = []; OUTPUT(NumStart).algorithm = []; OUTPUT(NumStart).firstorderopt = []; OUTPUT(NumStart).message = []; % Display run number if ismember(options.Display, {'final', 'iter'}) runNo = getRunIndex(optimstore); fprintf('\n\nRun %d\n', runNo); end % Run fmincon over each initial condition hObj = @i_evalObjective; hCon = @i_evalConstraint; for i = 1:NumStart y0 = i_evalObjective(x0(i, :), optimstore); c0 = i_evalConstraint(x0(i, :), optimstore); if isfinite(y0) && all(isfinite(c0(:))) [bestx(i, :), bestfun(i), exitFlag(i), outstruct] = fmincon(... hObj, x0(i, :), A, B, [],[], lb, ub, ... hCon, options,optimstore); OUTPUT(i).iterations = outstruct.iterations; OUTPUT(i).funcCount = outstruct.funcCount; OUTPUT(i).stepsize = outstruct.stepsize; OUTPUT(i).algorithm = outstruct.algorithm; OUTPUT(i).firstorderopt = outstruct.firstorderopt; OUTPUT(i).message = iGetShortMessage(outstruct.message); else % nonfinite initial objective bestx(i, :) = NaN; exitFlag(i) = -3; OUTPUT(i).iterations = 0; OUTPUT(i).funcCount = 0; OUTPUT(i).stepsize = 0; OUTPUT(i).algorithm = 'fmincon'; OUTPUT(i).firstorderopt = []; OUTPUT(i).cgiterations = []; OUTPUT(i).message = 'Start positions is not finite. Algorithm not run.'; end end % Choose the best run and get the results at that run [~, bestInd] = min(bestfun); bestx = bestx(bestInd, :); exitFlag = exitFlag(bestInd); OUTPUT = OUTPUT(bestInd); else % Not required to run the optimization. Return the initial conditions % as the "solution". bestx = getInitFreeVal(optimstore); exitFlag = -2; OUTPUT.iterations = 0; OUTPUT.funcCount = 0; OUTPUT.stepsize = 0; OUTPUT.algorithm = 'fmincon'; OUTPUT.firstorderopt = []; OUTPUT.cgiterations = []; OUTPUT.message = 'All start positions infeasible. Algorithm not run.'; end % Save the results optimstore = setFreeVariables(optimstore, bestx); optimstore = setExitStatus(optimstore, exitFlag, OUTPUT.message); OUTPUT = rmfield(OUTPUT, 'message'); optimstore = setOutput(optimstore, OUTPUT); function [y, yg] = i_evalObjective(x, optimstore) %I_EVALOBJECTIVE evaluate objective for fmincon if nargout>1 [y, yg] = evaluateObjective(optimstore, x); else y = evaluateObjective(optimstore, x); yg = []; end ObjectiveFuncTypes = getObjectiveType(optimstore); switch ObjectiveFuncTypes{1} case 'max' % Function and gradient need to be negated y = -y; yg = -yg; case 'min' otherwise error(message('mbc:mbcOSfmincon:InvalidState')); end function [c, ceq, cg, cgeq] = i_evalConstraint(x, optimstore) %I_EVALCONSTRAINT evaluate constraints for fmincon if nargout>2 % evaluate nonlinear inequality constraints [c, cg] = evaluateIneqCon(optimstore, x); % evaluate nonlinear equality constraints [ceq, cgeq] = evaluateEqCon(optimstore, x); else c = evaluateIneqCon(optimstore, x); ceq = evaluateEqCon(optimstore, x); end function stop = i_Stop(~, ~, ~,optimstore) %i_STOP output function for fmincon call %get stop state from optimstore. This is true after the stop or cancel %button has been pressed. stop = getStopState(optimstore); function x0 = i_generateStartPoints(optimstore, NumStart, RunFromFeasOnly) % Pre R2008a sessions may have set NumStart higher than 10000. Limit % NumStart to 10000. NumStart(NumStart > 10000) = 10000; % Tolerance for feasibility TOL = 1e-4; % Check whether the specified initial condition is feasible specifiedx0 = getInitFreeVal(optimstore); if RunFromFeasOnly sconval = i_evaluateConstraints(optimstore, specifiedx0); if ~all(sconval < TOL) specifiedx0 = []; end end if NumStart > 1 % Parameters for candidate set nDesign = 10000; nFree = sum(getNumFreeVariables(optimstore)); % Generate candidate set hSet = haltonset(nFree); hSet = scramble(hSet, 'RR2'); candx0 = net(hSet, nDesign); % Each variable is scaled on [0, 1]. Transform the variables back % to natural units lb = getLB(optimstore); ub = getUB(optimstore); shift = bsxfun(@times, ub-lb, candx0); candx0 = bsxfun(@plus, lb, shift); % Evaluate the constraints at each point in the candidate set. Filter out % those points that are constrained. conval = i_evaluateConstraints(optimstore, candx0); if isempty(conval) idxFeas = (1:size(candx0, 1))'; else idxFeas = find(all(conval < TOL, 2)); end % Evaluate the objectives over each feasible point in the candidate set. objval = evaluateObjective(optimstore, candx0(idxFeas, :)); ObjectiveType = getObjectiveType(optimstore); if strcmpi(ObjectiveType{1}, 'max') objval = -objval; end % Return the NumStart-1 feasible points with the best objective values if % possible nBestFeasStart = min(NumStart-1, size(objval, 1)); [~, idx] = sort(objval); thisx0 = candx0(idxFeas(idx(1:nBestFeasStart)), :); % If there are less than NumStart-1 feasible points, fill the remaining start % points with infeasible points with the lowest maximum constraint % violation. Only do this if the user wants to run from infeasible % solutions. nAddPoints = NumStart - 1 - nBestFeasStart; if ~RunFromFeasOnly && nAddPoints > 0 idxInfeas = setdiff(1:size(candx0, 1), idxFeas); maxConVal = max(conval(idxInfeas, :), [], 2); [~, idx] = sort(maxConVal); thisx0 = [thisx0; candx0(idxInfeas(idx(1:nAddPoints)), :)]; end else thisx0 = []; end % Return the initial conditions x0 = [specifiedx0; thisx0]; function conval = i_evaluateConstraints(optimstore, x) A = getA(optimstore); B = getB(optimstore); conval = evaluateNonlcon(optimstore, x); if ~isempty(A) lconval = A*x' - B(:, ones(1, size(x, 1))); else lconval = []; end conval = [conval, lconval']; function shortstr = iGetShortMessage(str) % Process output message from fmincon and extract the first short sentence. shortstr = regexp(str, '\n*(.*?)\n\n', 'once', 'tokens'); if ~isempty(shortstr) shortstr = shortstr{1}; else % Regexp did not match anything. Forward the string as-is. shortstr = str; end