compileinterfacedata.m

matlab波形优化算法经常要用到的matlab toolbox工具箱:yalmip

function [interfacedata,recoverdata,solver,diagnostic,F,Fremoved] = compileinterfacedata(F,aux_obsolete,logdetStruct,h,options,findallsolvers,parametric)

persistent CACHED_SOLVERS
persistent EXISTTIME
persistent NCHECKS

%% Initilize default empty outputs
diagnostic = [];
interfacedata = [];
recoverdata = [];
solver = [];
Fremoved  = [];
        
%% Did we make the call from SOLVEMP
if nargin<7
    parametric = 0;
end

%% Clean objective to default empty
if isa(h,'double')
    h = [];
end

% *************************************************************************
%% Exit if LOGDET objective is nonlinear
% *************************************************************************
if ~isempty(logdetStruct)
    for i = 1:length(logdetStruct.P)
        if ~is(logdetStruct.P{i},'linear')
            diagnostic.solvertime = 0;
            diagnostic.problem = -2;
            diagnostic.info = yalmiperror(diagnostic.problem,'');
            return
        end
    end
end

% *************************************************************************
%% EXTRACT LOW-RANK DESCRIPTION
% *************************************************************************
lowrankdetails = getlrdata(F);
if ~isempty(lowrankdetails)
    F = F(~is(F,'lowrank'));
end

% *************************************************************************
%% PERTURB STRICT INEQULAITIES
% *************************************************************************
if isa(options.shift,'sdpvar') | (options.shift~=0)
    F = shift(F,options.shift);
end

% *************************************************************************
%% ADD RADIUS CONSTRAINT
% *************************************************************************
if isa(options.radius,'sdpvar') | ~isinf(options.radius)
    x = recover(unique(union(depends(h),depends(F))));
    if length(x)>1
        F = F + set(cone(x,options.radius));
    else
        F = F + set(-options.radius < x < options.radius);
    end
end

% *************************************************************************
%% CONVERT LOGIC CONSTRAINTS
% *************************************************************************
[F,changed] = convertlogics(F);
if changed
    options.saveduals = 0; % Don't calculate duals since we changed the problem
end

% *************************************************************************
%% Take care of the nonlinear operators by converting expressions such as
% t = max(x,y) to standard conic models and mixed integer models
% This part also adds parts from logical expressions and mpower terms
% *************************************************************************
if options.expand
    % Experimental hack due to support for the PWQ function used for
    % quadratic dynamic programming with MPT.
    % FIX: Clean up and generalize
    try
        variables = uniquestripped([depends(h) getvariables(h)]);
        extendedvariables = yalmip('extvariables');
        index_in_extended = find(ismembc(variables,extendedvariables));
        if ~isempty(index_in_extended)
            extstruct = yalmip('extstruct',variables(index_in_extended));
            if isequal(extstruct.fcn ,'pwq_yalmip')
                [properties,Fz,arguments]=model(extstruct.var,'integer',options,extstruct);
                gain = getbasematrix(h,getvariables(extstruct.var));
                h = replace(h,extstruct.var,0);
                h = h + gain*properties.replacer;
                F = F + Fz;
            end
        end
    catch
    end
    [F,failure,cause,operators] = expandmodel(F,h,options);
    if failure % Convexity propgation failed
        interfacedata = [];
        recoverdata = [];
        solver = '';
        diagnostic.solvertime = 0;
        diagnostic.problem = 14;
        diagnostic.info = yalmiperror(14,cause);
        return
    end
    evalVariables = unique(determineEvaluationBased(operators));%yalmip('evalVariables');
    %evalVariables = yalmip('evalVariables');    
    if isempty(evalVariables)
        evaluation_based = 0;
    else
        evaluation_based = ~isempty(intersect([depends(h) depends(F)],evalVariables));
    end
else
    evalVariables = [];
    evaluation_based = 0;    
end

% *************************************************************************
%% CONVERT CONVEX QUADRATIC CONSTRAINTS
% We do not convert quadratic constraints to SOCPs if we have have
% sigmonial terms (thus indicating a GP problem), if we have relaxed
% nonlinear expressions, or if we have specified a nonlinear solver.
% Why do we convert them already here? Don't remember, should be cleaned up
% *************************************************************************
[monomtable,variabletype] = yalmip('monomtable');
F_vars = getvariables(F);
do_not_convert = any(variabletype(F_vars)==4);
do_not_convert = do_not_convert | strcmpi(options.solver,'snopt');
do_not_convert = do_not_convert | strcmpi(options.solver,'snopt-geometric') | strcmpi(options.solver,'snopt-standard');
do_not_convert = do_not_convert | strcmpi(options.solver,'ipopt');
do_not_convert = do_not_convert | strcmpi(options.solver,'ipopt-geometric') | strcmpi(options.solver,'ipopt-standard');
do_not_convert = do_not_convert | strcmpi(options.solver,'pennon');
do_not_convert = do_not_convert | strcmpi(options.solver,'pennon-geometric') | strcmpi(options.solver,'pennon-standard');
do_not_convert = do_not_convert | strcmpi(options.solver,'pennlp') | strcmpi(options.solver,'penbmi');
do_not_convert = do_not_convert | strcmpi(options.solver,'fmincon') | strcmpi(options.solver,'lindo') | strcmpi(options.solver,'sqplab');
do_not_convert = do_not_convert | strcmpi(options.solver,'fmincon-geometric') | strcmpi(options.solver,'fmincon-standard');
do_not_convert = do_not_convert | strcmpi(options.solver,'bmibnb');
do_not_convert = do_not_convert | (options.convertconvexquad == 0);
do_not_convert = do_not_convert | (options.relax == 1);
if ~do_not_convert & any(variabletype(F_vars))
    [F,socp_changed,infeasible] = convertquadratics(F);
    if infeasible
        diagnostic.solvertime = 0;
        diagnostic.problem = 1;
        diagnostic.info = yalmiperror(diagnostic.problem,'YALMIP');
        return        
    end
    if socp_changed % changed holds the number of QC -> SOCC conversions
        options.saveduals = 0; % We cannot calculate duals since we changed the problem
        F_vars = []; % We have changed model so we cannot use this in categorizemodel
    end
else
    socp_changed = 0;
end

% CHEAT FOR QC
if socp_changed>0 & length(find(is(F,'socc')))==socp_changed
    socp_are_really_qc = 1;
else
    socp_are_really_qc = 0;
end

% *************************************************************************
%% LOOK FOR AVAILABLE SOLVERS
% Finding solvers can be very slow on some systems. To alleviate this
% problem, YALMIP can cache the list of available solvers.
% *************************************************************************
if (options.cachesolvers==0) | isempty(CACHED_SOLVERS)
    getsolvertime = clock;
    solvers = getavailablesolvers(findallsolvers,options);
    getsolvertime = etime(clock,getsolvertime);
    % CODE TO INFORM USERS ABOUT SLOW NETWORKS!
    if isempty(EXISTTIME)
        EXISTTIME = getsolvertime;
        NCHECKS = 1;
    else
        EXISTTIME = [EXISTTIME getsolvertime];
        NCHECKS = NCHECKS + 1;
    end
    if (options.cachesolvers==0)
        if ((NCHECKS >= 3 & (sum(EXISTTIME)/NCHECKS > 1)) | EXISTTIME(end)>2)
            if warningon
                info = 'Warning: YALMIP has detected that your drive or network is unusually slow.\nThis causes a severe delay in SOLVESDP when I try to find available solvers.\nTo avoid this, use the options CACHESOLVERS in SDPSETTINGS.\nSee the FAQ for more information.\n';
                fprintf(info);
            end
        end
    end
    if length(EXISTTIME) > 5
        EXISTTIME = EXISTTIME(end-4:end);
        NCHECKS = 5;
    end
    CACHED_SOLVERS = solvers;
else
    solvers = CACHED_SOLVERS;
end

% *************************************************************************
%% NO SOLVER AVAILABLE
% *************************************************************************
if isempty(solvers)
    diagnostic.solvertime = 0;
    if isempty(options.solver)
        diagnostic.info = yalmiperror(-2,'YALMIP');
        diagnostic.problem = -2;
    else
        diagnostic.info = yalmiperror(-3,'YALMIP');
        diagnostic.problem = -3;
    end
    if warningon & options.warning & isempty(findstr(diagnostic.info,'No problems detected'))
        disp(['Warning: ' diagnostic.info]);
    end
    return
end

% *************************************************************************
%% WHAT KIND OF PROBLEM DO WE HAVE NOW?
% *************************************************************************
[ProblemClass,integer_variables,binary_variables,parametric_variables,uncertain_variables,quad_info] = categorizeproblem(F,logdetStruct,h,options.relax,parametric,evaluation_based,F_vars);

% *************************************************************************
%% SELECT SUITABLE SOLVER
% *************************************************************************
[solver,problem] = selectsolver(options,ProblemClass,solvers,socp_are_really_qc);
if isempty(solver)
    diagnostic.solvertime = 0;
    if problem == -4
        diagnostic.info = yalmiperror(problem,options.solver);
    else
        diagnostic.info = yalmiperror(problem,'YALMIP');
    end
    diagnostic.problem = problem;

    if warningon & options.warning
        disp(['Warning: ' diagnostic.info]);
    end
    return
end
if length(solver.version)>0
    solver.tag = [solver.tag '-' solver.version];
end

% *************************************************************************
%% DID WE SELECT THE INTERNAL BNB SOLVER
% IN THAT CASE, SELECT LOCAL SOLVER
% (UNLESS ALREADY SPECIFIED IN OPTIONS.BNB)
% *************************************************************************
localsolver.qc = 0;
localsolver = solver;
if strcmpi(solver.tag,'bnb')
    temp_options = options;
    temp_options.solver = options.bnb.solver;
    tempProblemClass = ProblemClass;
    tempProblemClass.constraint.binary  = 0;
    tempProblemClass.constraint.integer = 0;
    localsolver = selectsolver(temp_options,tempProblemClass,solvers,socp_are_really_qc);
    if isempty(localsolver) | strcmpi(localsolver.tag,'bnb')
        diagnostic.solvertime = 0;
        diagnostic.info = yalmiperror(-2,'YALMIP');
        diagnostic.problem = -2;
        return
    end
    solver.lower = localsolver;
end

% *************************************************************************
%% DID WE SELECT THE MPCVX SOLVER
% IN THAT CASE, SELECT SOLVER TO SOLVE BOUND COMPUTATIONS
% *************************************************************************
localsolver.qc = 0;
localsolver = solver;
if strcmpi(solver.tag,'mpcvx')