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