compileinterfacedata.m

optimization toolbox

    sdp = tempProblemClass.constraint.inequalities.semidefinite;
    tempProblemClass.constraint.inequalities.semidefinite.linear = sdp.linear | sdp.quadratic | sdp.polynomial;
    tempProblemClass.constraint.inequalities.semidefinite.quadratic = 0;
    tempProblemClass.constraint.inequalities.semidefinite.polynomial = 0;

    lp = tempProblemClass.constraint.inequalities.elementwise;
    tempProblemClass.constraint.inequalities.elementwise.linear = lp.linear | lp.quadratic.convex | lp.quadratic.nonconvex | sdp.polynomial;
    tempProblemClass.constraint.inequalities.elementwise.quadratic.convex = 0;
    tempProblemClass.constraint.inequalities.elementwise.quadratic.nonconvex = 0;
    tempProblemClass.constraint.inequalities.elementwise.polynomial = 0;

    equ = tempProblemClass.constraint.equalities;
    tempProblemClass.constraint.equalities.linear = equ.linear | equ.quadratic | equ.polynomial;
    tempProblemClass.constraint.equalities.quadratic = 0;
    tempProblemClass.constraint.equalities.polynomial = 0;

    %     if nnz(quad_info.Q) > 0
    %     if all(eig(quad_info.Q) > -1e-12)
    %     tempProblemClass.objective.quadratic.convex = 0;
    %     tempProblemClass.objective.quadratic.nonconvex = 0;
    %     else
    %    tempProblemClass.objective.quadratic.convex = 0;
    tempProblemClass.objective.quadratic.nonconvex = 0;
    tempProblemClass.objective.polynomial = 0;
    %    end
    %    end

    tempProblemClass.constraint.inequalities.rank  = 0;
    tempProblemClass.evaluation  = 0;

    temp_options = options;
    temp_options.solver = options.bmibnb.lowersolver;

    % If the problem actually is quadratic, try to get a convex problem
    % this will typically allow us to solver better lower bounding problems
    % (we don't have to linearize the cost)
    [lowersolver,problem] = selectsolver(temp_options,tempProblemClass,solvers,socp_are_really_qc);
    if isempty(lowersolver)| strcmpi(lowersolver.tag,'bmibnb') | strcmpi(lowersolver.tag,'bnb')
        % No, probably non-convex cost. Pick a linear solver instead and go
        % for lower bound based on a complete "linearization"
        tempProblemClass.objective.quadratic.convex = 0;
        [lowersolver,problem] = selectsolver(temp_options,tempProblemClass,solvers,socp_are_really_qc);
    end

    if isempty(lowersolver)| strcmpi(lowersolver.tag,'bmibnb') | strcmpi(lowersolver.tag,'bnb')
        tempbinary = tempProblemClass.constraint.binary;
        tempinteger = tempProblemClass.constraint.integer;
        tempProblemClass.constraint.binary = 0;
        tempProblemClass.constraint.integer = 0;
        [lowersolver,problem] = selectsolver(temp_options,tempProblemClass,solvers,socp_are_really_qc);
        tempProblemClass.constraint.binary = tempbinary;
        tempProblemClass.constraint.integer = tempinteger;
    end

    if isempty(lowersolver) | strcmpi(lowersolver.tag,'bmibnb') | strcmpi(lowersolver.tag,'bnb')
        diagnostic.solvertime = 0;
        diagnostic.info = yalmiperror(-2,'YALMIP');
        diagnostic.problem = -2;
        return
    end
    solver.lowercall = lowersolver.call;
    solver.lowersolver = lowersolver;

    temp_options = options;
    temp_options.solver = options.bmibnb.uppersolver;
    temp_ProblemClass = ProblemClass;
    temp_ProblemClass.constraint.binary = 0;
    temp_ProblemClass.constraint.integer = 0;
   % temp_ProblemClass.objective.quadratic.convex = 0;
   % temp_ProblemClass.objective.quadratic.nonconvex = 0;
    [uppersolver,problem] = selectsolver(temp_options,temp_ProblemClass,solvers,socp_are_really_qc);
    if strcmpi(uppersolver.tag,'bnb')
        temp_options.solver = 'none';
        [uppersolver,problem] = selectsolver(temp_options,temp_ProblemClass,solvers,socp_are_really_qc);
    end
    if isempty(uppersolver) | strcmpi(uppersolver.tag,'bmibnb')
        diagnostic.solvertime = 0;
        diagnostic.info = yalmiperror(-2,'YALMIP');
        diagnostic.problem = -2;
        return
    end
    solver.uppercall = uppersolver.call;
    solver.uppersolver = uppersolver;


    temp_options = options;
    temp_options.solver = options.bmibnb.lpsolver;
    tempProblemClass.constraint.inequalities.semidefinite.linear = 0;
    tempProblemClass.constraint.inequalities.semidefinite.quadratic = 0;
    tempProblemClass.constraint.inequalities.semidefinite.polynomial = 0;
    tempProblemClass.constraint.inequalities.secondordercone = 0;
    tempProblemClass.objective.quadratic.convex = 0;
    tempProblemClass.objective.quadratic.nonconvex = 0;
    tempProblemClass.objective.quadratic.nonconvex = 0;
    tempProblemClass.objective.polynomial = 0;

    [lpsolver,problem] = selectsolver(temp_options,tempProblemClass,solvers,socp_are_really_qc);

    if isempty(lowersolver)| strcmpi(lowersolver.tag,'bmibnb')
        tempbinary = tempProblemClass.constraint.binary;
        tempProblemClass.constraint.binary = 0;
        [lpsolver,problem] = selectsolver(temp_options,tempProblemClass,solvers,socp_are_really_qc);
        tempProblemClass.constraint.binary = tempbinary;
    end

    if isempty(lpsolver) | strcmpi(lpsolver.tag,'bmibnb')
        diagnostic.solvertime = 0;
        diagnostic.info = yalmiperror(-2,'YALMIP');
        diagnostic.problem = -2;
        return
    end
    solver.lpsolver = lpsolver;
    solver.lpcall = lpsolver.call;
end

% *************************************************************************
% DID WE SELECT THE INTERNAL SDPMILP SOLVER
% This solver solves MISDP problems by solving MILP problems and adding SDP
% cuts based on the infasible MILP solution.
% *************************************************************************
if strcmpi(solver.tag,'cutsdp')

    % Relax problem for lower solver
    tempProblemClass = ProblemClass;
    tempProblemClass.constraint.inequalities.elementwise.linear =  tempProblemClass.constraint.inequalities.elementwise.linear |     tempProblemClass.constraint.inequalities.semidefinite.linear | tempProblemClass.constraint.inequalities.secondordercone;
    tempProblemClass.constraint.inequalities.semidefinite.linear = 0;
    tempProblemClass.constraint.inequalities.secondordercone = 0;

    temp_options = options;
    temp_options.solver = options.cutsdp.solver;

    [cutsolver,problem] = selectsolver(temp_options,tempProblemClass,solvers,socp_are_really_qc);

    if isempty(cutsolver) | strcmpi(cutsolver.tag,'cutsdp') |strcmpi(cutsolver.tag,'bmibnb') | strcmpi(cutsolver.tag,'bnb')
        diagnostic.solvertime = 0;
        diagnostic.info = yalmiperror(-2,'YALMIP');
        diagnostic.problem = -2;
        return
    end
    solver.cutsolver = cutsolver;
end

showprogress(['Solver chosen : ' solver.tag],options.showprogress);

% *************************************************************************
% CONVERT MAXDET TO SDP USING GEOMEAN?
% *************************************************************************
% MAXDET using geometric mean construction
if ~isempty(P) & solver.objective.maxdet==0
    t = sdpvar(1,1);
    if isequal(P,sdpvar(F(end)))
        F = F(1:end-1);
    end
    F = F + detset(t,P);
    if isempty(h)
        h = -t;
    else
        h = h-t;
        % Warn about logdet -> det^1/m
        if options.verbose>0 & options.warning>0
            disp(' ')
            disp('Objective c''x-logdet(P) has been changed to c''x-det(P)^(1/(2^ceil(log2(length(X))))).')
            disp('See the MAXDET section in the manual for details.')
            disp(' ')
        end
    end
end

% *************************************************************************
% Change binary variables to integer?
% *************************************************************************
old_binary_variables = binary_variables;
if ~isempty(binary_variables) & (solver.constraint.binary==0)
    x_bin = recover(binary_variables(ismember(binary_variables,unique([getvariables(h) getvariables(F)]))));
    F = F + set(x_bin<1)+set(x_bin>0);
    integer_variables = union(binary_variables,integer_variables);
    binary_variables = [];
end

% *************************************************************************
% Model quadratics using SOCP?
% Should not be done when using PENNLP or BMIBNB, or if we have relaxed the
% monmial terms or...Ouch, need to clean up all special cases, this sucks.
% *************************************************************************
convertQuadraticObjective = ~strcmpi(solver.tag,'pennlp-standard');
convertQuadraticObjective = convertQuadraticObjective & ~strcmpi(solver.tag,'bmibnb');
convertQuadraticObjective = convertQuadraticObjective & ( ~(options.relax==1 | options.relax==3)  &(~isempty(quad_info) & solver.objective.quadratic.convex==0 | (~isempty(quad_info) & strcmp(solver.tag,'bnb') & localsolver.objective.quadratic.convex==0)));
if strcmpi(solver.tag,'bnb') & ~isempty(quad_info)
    if solver.lower.objective.quadratic.convex==0
        convertQuadraticObjective = 1;
    end
end
if convertQuadraticObjective
    t = sdpvar(1,1);
    x = quad_info.x;
    R = quad_info.R;
    if ~isempty(R)
        c = quad_info.c;
        f = quad_info.f;
        F = F + lmi(cone([2*R*x;1-(t-c'*x-f)],1+t-c'*x-f));
        h = t;
    end
    quad_info = [];
end
if solver.constraint.inequalities.rotatedsecondordercone == 0
    [F,changed] = convertlorentz(F);
    if changed
        options.saveduals = 0; % We cannot calculate duals since we change the problem
    end
end
% Whoa, horrible tests to find out when to convert SOCP to SDP
% This should not be done if :
%   1. Problem is actually a QCQP and solver supports this
%   2. Problem is integer, local solver supports SOCC
%   3. Solver supports SOCC
if ~((solver.constraint.inequalities.elementwise.quadratic.convex == 1) & socp_are_really_qc)
    if ~(strcmp(solver.tag,'bnb') & socp_are_really_qc & localsolver.constraint.inequalities.elementwise.quadratic.convex==1 )
        if ((solver.constraint.inequalities.secondordercone == 0) | (strcmpi(solver.tag,'bnb') & localsolver.constraint.inequalities.secondordercone==0))
            [F,changed] = convertsocp(F);
            if changed
                options.saveduals = 0; % We cannot calculate duals since we change the problem
            end
        end
    end
end

% *************************************************************************
% Add logaritmic barrier cost/constraint for MAXDET and SDPT3-4. Note we
% have to add it her ein order for a complex valued matrix to be converted.
% *************************************************************************
if ~isempty(P) & solver.objective.maxdet==1
    F = F + set(P > 0);
end

if ((solver.complex==0) & ProblemClass.complex) | ((strcmp(solver.tag,'bnb') & localsolver.complex==0) & ProblemClass.complex)
    showprogress('Converting to real constraints',options.showprogress)
    F = imag2reallmi(F);
    options.saveduals = 0; % We cannot calculate duals since we change the problem
end

% *************************************************************************
% CREATE OBJECTIVE FUNCTION c'*x+x'*Q*x
% *************************************************************************
showprogress('Processing objective h(x)',options.showprogress);
try
    if strcmpi(solver.tag,'bmibnb') | strcmpi(solver.tag,'pennlp-standard')
        tempoptions = options;
        tempoptions.relax = 1;
        [c,Q,f]=createobjective(h,P,tempoptions,quad_info);
    else
        [c,Q,f]=createobjective(h,P,options,quad_info);
    end
catch
    error(lasterr)
end

% *************************************************************************
% Convert {F(x),G(x)} to a numerical SeDuMi-like format
% *************************************************************************
showprogress('Processing F(x)',options.showprogress);
[F_struc,K,KCut] = lmi2sedumistruct(F);
% We add a field to remember the dimmesion of the logarithmic cost.
% Actually, the actually value is not interesting, we know that the
% logarithmic cost corresponds to the last SDP constraint anyway
K.m = length(P);

% *************************************************************************
% SOME HORRIBLE CODE TO DETERMINE USED VARIABLES
% *************************************************************************
% Which sdpvar variables are actually in the problem
used_variables_LMI = find(any(F_struc(:,2:end),1));
used_variables_obj = find(any(c',1) | any(Q));
used_variables = uniquestripped([used_variables_LMI used_variables_obj]);
% The problem is that linear terms might be missing in problems with only
% nonlinear expressions
[monomtable,variabletype] = yalmip('monomtable');
if (options.relax==1)|(options.relax==3)
    monomtable = [];
    nonlinearvariables = [];
    linearvariables = used_variables;
else
    nonlinearvariables = find(variabletype);
    linearvariables = used_variables(find(variabletype(used_variables)==0));
end
needednonlinear = nonlinearvariables(ismembc(nonlinearvariables,used_variables));
linearinnonlinear = find(sum(abs(monomtable(needednonlinear,:)),1));
missinglinear = setdiff(linearinnonlinear(:),linearvariables);
used_variables = uniquestripped([used_variables(:);missinglinear(:)]);

% *************************************************************************
% So are we done now? No... What about variables hiding inside so called
% evaluation variables. We detect these, and at the same time set up the
% structures needed to support general functions such as exp, log, etc
% NOTE : This is experimental code
% FIX  : Clean up...
% *************************************************************************
[evalMap,evalVariables,used_variables,nonlinearvariables,linearvariables] = detectHiddenNonlinear(used_variables,options,nonlinearvariables,linearvariables);

% *************************************************************************
% REMOVE UNNECESSARY VARIABLES FROM PROBLEM
% *************************************************************************
if length(used_variables)<yalmip('nvars')
    c = c(used_variables);
    Q = Q(:,used_variables);Q = Q(used_variables,:);
    if ~isempty(F_struc)
        F_struc = sparse(F_struc(:,[1 1+used_variables]));
    end
end

% *************************************************************************
% Map variables and constraints in low-rank definition to local stuff
% *************************************************************************
if ~isempty(lowrankdetails)