compileinterfacedata.m

optimization toolbox

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

persistent CACHED_SOLVERS
persistent EXISTTIME
persistent NCHECKS

diagnostic = [];
interfacedata = [];
recoverdata = [];
solver = [];
Fremoved  = [];

if nargin<7
    parametric = 0;
end

if isequal(h,0)
    h = [];
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')
                [Fz,properties,arguments]=model(extstruct.var,'milp',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] = 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 = yalmip('evalVariables');
    if isempty(evalVariables)
        evaluation_based = 0;
    else
        evaluation_based = ~isempty(intersect([getvariables(h) getvariables(F)],evalVariables));
    end
else
    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');
% if any(variabletype(getvariables(F))==4)
%     do_not_convert = 1;
% else
%     do_not_convert = 0;
% end
do_not_convert = any(variabletype(getvariables(F))==4);
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');
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
    [F,socp_changed] = convertquadratics(F);
    if socp_changed % changed holds the number of QC -> SOCC conversions
        options.saveduals = 0; % We cannot calculate duals since we changed the problem
    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 avaialble 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,P,h,options.relax,parametric,evaluation_based);

% *************************************************************************
% 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')
    temp_options = options;
    temp_options.solver = options.mpcvx.solver;
    tempProblemClass = ProblemClass;    
    tempProblemClass.objective.quadratic.convex = tempProblemClass.objective.quadratic.convex | tempProblemClass.objective.quadratic.nonconvex;
    tempProblemClass.objective.quadratic.nonconvex = 0;
    tempProblemClass.parametric = 0;
    localsolver = selectsolver(temp_options,tempProblemClass,solvers,socp_are_really_qc);
    if isempty(localsolver) | strcmpi(localsolver.tag,'bnb') | strcmpi(localsolver.tag,'kktqp')
        diagnostic.solvertime = 0;
        diagnostic.info = yalmiperror(-2,'YALMIP');
        diagnostic.problem = -2;
        return
    end
    solver.lower = localsolver;
end

% *************************************************************************
% DID WE SELECT THE INTERNAL EXPERIMENTAL KKT SOLVER
% IN THAT CASE, SELECT SOLVER TO SOLVE THE MILP PROBLEM
% *************************************************************************
localsolver.qc = 0;
localsolver = solver;
if strcmpi(solver.tag,'kktqp')
    temp_options = options;
    temp_options.solver = '';
    tempProblemClass = ProblemClass;
    tempProblemClass.constraint.binary = 1;
    tempProblemClass.objective.quadratic.convex = 0;
    tempProblemClass.objective.quadratic.nonconvex = 0;
    localsolver = selectsolver(temp_options,tempProblemClass,solvers,socp_are_really_qc);
    if isempty(localsolver) | strcmpi(localsolver.tag,'bnb') | strcmpi(localsolver.tag,'kktqp')
        diagnostic.solvertime = 0;
        diagnostic.info = yalmiperror(-2,'YALMIP');
        diagnostic.problem = -2;
        return
    end
    solver.lower = localsolver;
end

% *************************************************************************
% DID WE SELECT THE LMIRANK?
% FIND SDP SOLVER FOR INITIAL SOLUTION
% *************************************************************************
if strcmpi(solver.tag,'lmirank')
    temp_options = options;
    temp_options.solver = options.lmirank.solver;
    tempProblemClass = ProblemClass;
    tempProblemClass.constraint.inequalities.rank = 0;
    tempProblemClass.constraint.inequalities.semidefinite.linear = 1;
    tempProblemClass.objective.linear = 1;
    initialsolver = selectsolver(temp_options,tempProblemClass,solvers,socp_are_really_qc);
    if isempty(initialsolver) | strcmpi(initialsolver.tag,'lmirank')
        diagnostic.solvertime = 0;
        diagnostic.info = yalmiperror(-2,'YALMIP');
        diagnostic.problem = -2;
        return
    end
    solver.initialsolver = initialsolver;
end

% *************************************************************************
% DID WE SELECT THE INTERNAL BMIBNB SOLVER? SELECT UPPER/LOWER SOLVERs
% (UNLESS ALREADY SPECIFIED IN OPTIONS)
% *************************************************************************
if strcmpi(solver.tag,'bmibnb')

    % Relax problem for lower solver
    tempProblemClass = ProblemClass;