compileinterfacedata.m

国外专家做的求解LMI鲁棒控制的工具箱,可以相对高效的解决LMI问题

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

persistent CACHED_SOLVERS
persistent EXISTTIME
persistent NCHECKS

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

% ******************************************
% 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; % We cannot calculate duals since we changed the problem
end

% ******************************************
% Take care of the compound objects
% by convexifying stuff like t = max(x,y)
% ******************************************
if ~isempty(yalmip('extvariables'))
    if 1
    [F,failure] = convexitypropagation(F,h);
    else
    [F,failure] = newconvexitypropagation(F,h);
    end
    if failure
        interfacedata = [];
        recoverdata = [];
        solver = '';
        diagnostic.solvertime = 0;
        diagnostic.problem = 14;
        diagnostic.info = yalmiperror(14,'YALMIP');
        return
    end
end

% *********************************************
% CONVERT CONVEX QUADRATIC CONSTRAINTS
% *********************************************
if (options.convertconvexquad) & ~(strcmpi(options.solver,'pennlp')) & ~(strcmpi(options.solver,'fmincon-geometric')) & ~(strcmpi(options.solver,'bmibnb')) &~(options.relax==1) & ~(options.relax==3)
    [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

% ******************************************
% LOOK FOR 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)  & ((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
    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

% 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
    
% ******************************************
% TRY TO SELECT SUITABLE SOLVER
% ******************************************
[ProblemClass,integer_variables,binary_variables,parametric_variables,quad_info] = newcatsdp(F,G,h,options.relax);

[solver,problem] = newselectsolver(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 = newselectsolver(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.subcall = localsolver.call;
    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 = newselectsolver(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.subcall = localsolver.call;
    solver.initialsolver = initialsolver;
end


% *****************************************
% CONVERT POLYNOMIAL PROGRAM TO BILINEAR
% *****************************************
replacements = [];
violates = ProblemClass.objective.polynomial & (~(solver.objective.polynomial | solver.objective.sigmonial));
violates = violates | (ProblemClass.constraint.inequalities.elementwise.polynomial & (~(solver.constraint.inequalities.elementwise.polynomial | solver.constraint.inequalities.elementwise.sigmonial)));
violates = violates | (ProblemClass.constraint.inequalities.semidefinite.polynomial & (~solver.constraint.inequalities.semidefinite.polynomial));
violates = violates | (ProblemClass.constraint.equalities.polynomial & (~solver.constraint.equalities.polynomial));
if violates
    % Bilinearize?
    mt = yalmip('monomtable');
    vars = union(getvariables(h),getvariables(F));
    if any((sum(mt(vars,:),2))>2)
        % Bummer, high order stuff
        F_new = set([]);
        violates = find((sum(mt(vars,:),2))>2);
        replacements = [];
        for i = 1:length(violates)
            these = find(mt(vars(violates(i)),:));
            pow = [];
            for j = 1:length(these)
                pow = [pow  repmat(these(j),1,mt(vars(violates(i)),these(j)))];
            end
            [w,F_new] = recdef(pow,F_new,[],options.usex0);
            replacements = [replacements;vars(violates(i)) getvariables(w)];
        end
        F = F + F_new;
    end
    ProblemClass.objective.polynomial = 0;
    ProblemClass.objective.inequalities.elementswise.polynomial = 0;
    ProblemClass.objective.equality.quadratic = 1;
end

% *****************************************
% DID WE SELECT THE INTERNAL BMIBNB SOLVER
% IN THAT CASE, SELECT UPPER/LOWER SOLVERs
% (UNLESS ALREADY SPECIFIED IN OPTIONS)
% *****************************************
if strcmpi(solver.tag,'bmibnb')
        
    % Relax problem for lower solver
    tempProblemClass = ProblemClass;    

    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;
    
    tempProblemClass.objective.quadratic.convex = 0;
    tempProblemClass.objective.quadratic.nonconvex = 0;
    
    tempProblemClass.constraint.inequalities.rank  = 0;
        
    temp_options = options;
    temp_options.solver = options.bmibnb.lowersolver;    

    [lowersolver,problem] = newselectsolver(temp_options,tempProblemClass,solvers,socp_are_really_qc);

    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] = newselectsolver(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] = newselectsolver(temp_options,temp_ProblemClass,solvers,socp_are_really_qc);       
    if isempty(uppersolver) | strcmp(lower(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;

    [lpsolver,problem] = newselectsolver(temp_options,tempProblemClass,solvers,socp_are_really_qc);
    
    if isempty(lowersolver)| strcmp(lower(lowersolver.tag),'bmibnb')
        tempbinary = tempProblemClass.constraint.binary;
        tempProblemClass.constraint.binary = 0;
        [lpsolver,problem] = newselectsolver(temp_options,tempProblemClass,solvers,socp_are_really_qc);
        tempProblemClass.constraint.binary = tempbinary;
    end
    
    if isempty(lpsolver) | strcmp(lower(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
% *****************************************
if strcmp(lower(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] = newselectsolver(temp_options,tempProblemClass,solvers,socp_are_really_qc);

    if isempty(cutsolver) | strcmp(lower(cutsolver.tag),'cutsdp') |strcmp(lower(cutsolver.tag),'bmibnb') | strcmp(lower(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 TO STANDARD PROBlEMS