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

    function out = mbcOSmodalopt(action, in)
%MBCOSMODALOPT CAGE Optimization library function
%   MBCOSMODALOPT contains the configuration and evaluation routines for
%   optimal selection of operating modes in CAGE. This function is called
%   by CAGE using the syntaxes below.
%
%   OPTOBJ = MBCOSMODALOPT('OPTIONS', OPTOBJ) is called by CAGE when an
%   Mixed 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 = MBCOSMODALOPT('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 MBCOSFMINCON

%  Copyright 2009-2015 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

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, 'Modal optimization');

% Add a description
opt = setDescription(opt, 'Modal optimization using a single objective optimization subject to constraints');

% Set up the free variables
%%%% Any number are allowed %%%%

% Set up the objective functions
%%%% Any objective is allowed %%%%
opt = setObjectivesMode(opt,'any');

% 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
%%%% No operating point sets are allowed %%%%
opt = setOperatingPointsMode(opt, 'fixed');

% 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, 'ModeVariable',... 
    {'integer', 'positive'},1,...
    'Index to the mode free variable');
opt = addParameter(opt, 'ModeObjective',... 
    {'integer', 'positive'},1,...
    'Index to the objective to determine best mode');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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);

ModeVariable = getParam(optimstore, 'ModeVariable');

try 
    ModeObjective = getParam(optimstore, 'ModeObjective');
catch
    % old versions don't have ModeObjective
    ModeObjective = 1;
end

% Get the user defined options and put them in an optimization toolbox
% options structure for fmincon
options = optimset(optimstore);
doDisplay = ~strcmp(options.Display,'none') && ~getStopState(optimstore);

x0 = getInitFreeVal(optimstore);
if isScalarFreeVariables(optimstore)
    % modal optimization only works for scalar problems
    
    % Set some extra options in fmincon
    options = optimset(options, ...
        'OutputFcn',@i_Output,...
        'GradObj', 'off', ...
        'GradCon', 'off');

    if ModeVariable>numel(x0)
        error(message('mbc:mbcOSmodalopt:InvalidValue'))
    end
    NumModes = ub(ModeVariable);
    MinMode = lb(ModeVariable);

    if NumModes~=fix(NumModes) || NumModes<1 || ~isfinite(NumModes) || ...
            MinMode~=fix(MinMode) || MinMode~=1 || ~isfinite(MinMode)
       error(message('mbc:mbcOSmodalopt:InvalidValue1'))
    end
    
    IndexCont = (1:length(x0))~=ModeVariable;
    % Initialise return quantities
    bestx = repmat(x0,NumModes,1);
    bestfun = Inf(NumModes, 1);
    exitFlag = zeros(1,NumModes);
    
    % make output structure
    OUTPUT(NumModes).iterations = [];
    OUTPUT(NumModes).funcCount = [];
    OUTPUT(NumModes).stepsize = [];
    OUTPUT(NumModes).algorithm = [];
    OUTPUT(NumModes).firstorderopt = [];
    OUTPUT(NumModes).message = [];

    % remove Mode Variable from bounds and linear constraints
    LB = lb(IndexCont);
    UB = ub(IndexCont);
    if ~isempty(A)
        A = A(IndexCont,:);
        B = B(IndexCont);
    end
    
    Xfull = x0;
    runNo = getRunIndex(optimstore);
    if doDisplay
        fprintf('Run %d\n',runNo);
    end
    Objectives = getObjectives(optimstore);
    ObjectiveFuncTypes = getObjectiveType(optimstore);
    if ModeObjective>length(Objectives)
        error(message('mbc:mbcOSmodalopt:InvalidValue2'))
    end

    for Mode = 1:NumModes
        Xfull(ModeVariable) = Mode;
        Xmode0 = x0(:, IndexCont);
        if getStopState(optimstore)
            y0= NaN;
        else
            y0 = i_evalObj(Xmode0,Xfull,IndexCont,optimstore);
        end 
        bestx(Mode, ModeVariable)  = Mode;
        if isfinite(y0)
            % optimize for each mode
            [bestx(Mode, IndexCont), bestfun(Mode), exitFlag(Mode), outstruct] = fmincon(@i_evalObj, ...
                Xmode0, A, B, [],[], LB, UB, @i_evalCon, options,Xfull,IndexCont,optimstore);
            if ModeObjective>1
                % evaluate a different objective as mode objective
                bestfun(Mode) = evaluateObjective(optimstore, bestx(Mode,:),Objectives(ModeObjective));
                if strcmp(ObjectiveFuncTypes{ModeObjective},'max')
                    bestfun(Mode) = -bestfun(Mode);
                end
            end
            
            if doDisplay
                if exitFlag(Mode)<0
                    fprintf('  Mode %d: Not optimal\n', Mode);
                elseif exitFlag(Mode)==0
                    fprintf('  Mode %d: Incomplete %g\n', Mode,bestfun(Mode));
                else
                    fprintf('  Mode %d: Optimal %g\n', Mode,bestfun(Mode));
                end
            end
            OUTPUT(Mode).iterations = outstruct.iterations;
            OUTPUT(Mode).funcCount = outstruct.funcCount;
            OUTPUT(Mode).stepsize = outstruct.stepsize;
            OUTPUT(Mode).algorithm = outstruct.algorithm;
            OUTPUT(Mode).firstorderopt = outstruct.firstorderopt;
            OUTPUT(Mode).message = sprintf('Mode %d: %s',Mode,iGetShortMessage(outstruct.message));
        else
            bestx(Mode, IndexCont) = NaN;
            bestfun(Mode) = Inf;
            if getStopState(optimstore)
                % user stopped optimization
                OUTPUT(Mode).message = sprintf('Mode %d: Modal optimization stopped by the output or plot function.',Mode);
                exitFlag(Mode) = -1;
            else
                % most likely the optimization didn't run due to an invalid
                % operating point for the switch model
                if doDisplay
                    fprintf('  Mode %d: Invalid\n', Mode);
                end
                OUTPUT(Mode).message = sprintf('Mode %d: Invalid\n', Mode);
                exitFlag(Mode) = -7;
            end
            OUTPUT(Mode).iterations = 0;
            OUTPUT(Mode).funcCount = 0;
            OUTPUT(Mode).stepsize = 0;
            OUTPUT(Mode).algorithm = 'fmincon';
            OUTPUT(Mode).firstorderopt = [];
            OUTPUT(Mode).cgiterations = [];
        end

    end

    % Choose the best run and get the results at that run
    if any(exitFlag>0)
        % only use successful optim modes 
        IndOK = find(exitFlag>0);    
        [~, bestInd] = min(bestfun(IndOK));
        SelectedSolution = IndOK(bestInd);
    elseif any(exitFlag>-7)
        % use all modes
        IndOK = find(exitFlag>-7);    
        [~, bestInd] = min(bestfun(IndOK));
        SelectedSolution = IndOK(bestInd);
    else
        [~, SelectedSolution] = min(bestfun);
    end
    if doDisplay
        fprintf('Best Mode %d: %g \n\n', SelectedSolution,bestfun(SelectedSolution));
    end

else
    
    % Not required to run the optimization. Return the initial conditions
    % as the "solution".
    bestx = getInitFreeVal(optimstore);
    exitFlag = -8;
    OUTPUT.iterations = 0;
    OUTPUT.funcCount = 0;
    OUTPUT.stepsize = 0;
    OUTPUT.algorithm = 'fmincon';
    OUTPUT.firstorderopt = [];
    OUTPUT.cgiterations = [];
    OUTPUT.message = 'Modal optimization cannot be used for nonscalar problems. Algorithm not run.';
    SelectedSolution = [];
end    

% Save the results
optimstore = setFreeVariables(optimstore, bestx,SelectedSolution);
optimstore = setExitStatus(optimstore, exitFlag, {OUTPUT.message});

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Options plus objective and constraint evaluation %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [y, yg] = i_evalObj(x,Xfull,IndexCont,optimstore)

Xfull(IndexCont) = x(:);
yg = [];
Objectives = getObjectives(optimstore);
y = evaluateObjective(optimstore, Xfull,Objectives(1));
ObjectiveFuncTypes = getObjectiveType(optimstore);
switch ObjectiveFuncTypes{1}
    case 'max'
        % Function and gradient need to be negated
        y = -y;
    case 'min'
    otherwise
        error(message('mbc:mbcOSfmincon:InvalidState'));
end
end

function [c, ceq, cg, cgeq] = i_evalCon(x,Xfull,IndexCont,optimstore)

Xfull(IndexCont) = x(:);
if nargout>2
    % evaluate nonlinear inequality constraints
    [c, cg] = evaluateIneqCon(optimstore, Xfull);
    % evaluate nonlinear equality constraints
    [ceq, cgeq] = evaluateEqCon(optimstore, Xfull);
else
    c = evaluateIneqCon(optimstore, Xfull);
    ceq = evaluateEqCon(optimstore, Xfull);
end

end


function stop= i_Output(varargin)
% get stop state so the optimization can stop in the middle of run
stop= getStopState(varargin{end});
end


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
end