compileinterfacedata.m

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

% ******************************************

% 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

% Relax binary valiables
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
convertQuadraticObjective = ~strcmp(lower(solver.tag),'pennlp-standard');
convertQuadraticObjective = convertQuadraticObjective & ~strcmp(lower(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 convertQuadraticObjective & 
if strcmp(lower(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
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
% **********************************************
if isequal(lower(solver.tag),'pennlp-standard')
    allx = recover(uniquestripped([depends(F) depends(h)]))
    dh  = jacobian(h,allx);
    ddh = jacobian(dh(:),allx);
    allconstraints = [];
    for i = 1:length(F)
        Fi = sdpvar(F(i));
        allconstraints = [allconstraints;Fi(:)];
    end
    dg = jacobian(allconstraints(:),allx);
    ddg = jacobian(dg(:),allx);
    
    gradientsHessians = [dh(:);ddh(:);dg(:);ddg(:)]; 
else
    gradientsHessians = [];
end

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

% **********************************************
% Convert {F(x),G(x)} to internal formats 
% **********************************************
showprogress('Processing F(x)',options.showprogress);
[F_struc,K,KCut] = lmi2sedumistruct(F,gradientsHessians);
K.m = length(sdpvar(G));
if ~isempty(replacements)
    n_orig_eq = K.f-length(F_new);  
    newcon = setdiff(1:K.f,1:n_orig_eq);
    replaceinthese = setdiff(1:size(F_struc,1),newcon);
    for i = 1:size(replacements,1)
        c(replacements(i,2)) = c(replacements(i,1));
        c(replacements(i,1)) = 0;
        F_struc(replaceinthese,1+replacements(i,2)) = F_struc(replaceinthese,1+replacements(i,1));
        F_struc(replaceinthese,1+replacements(i,1)) = 0;
    end
end

% **********************************************
% 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
monomtable = yalmip('monomtable');
if (options.relax==1)|(options.relax==3)
    monomtable = [];
    nonlinearvariables = [];
else
    nonlinearvariables = find(~(sum(monomtable~=0,2)==1 & sum(monomtable,2)==1));
end
linearvariables = setdiff(used_variables,nonlinearvariables);
needednonlinear = nonlinearvariables(ismember(nonlinearvariables,used_variables));
%linearinnonlinear = find(any(full(monomtable(needednonlinear,:)),1));
linearinnonlinear = find(sum(abs(monomtable(needednonlinear,:)),1));
missinglinear = setdiff(linearinnonlinear(:),linearvariables);
used_variables = unique(sort([used_variables(:);missinglinear(:)])');

% Check for unbounded variables
if 0%isempty(sqrList) WHY did I remove these
    unbounded_variables = setdiff(used_variables_obj,used_variables_LMI);
    if ~isempty(unbounded_variables)
        % Remove unbounded variable from problem
        used_variables = setdiff(used_variables,unbounded_variables);
    end
end

% ****************************************
% 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)
    % Identifiers of the SDP constraints
    lmiid = getlmiid(F);
    %lmiid = find(ismember(lmiid(is(F,'sdp'),);
    for i = 1:length(lowrankdetails)
        lowrankdetails{i}.id = find(ismember(lmiid,lowrankdetails{i}.id));
        if ~isempty(lowrankdetails{i}.variables)
            lowrankdetails{i}.variables = find(ismember(used_variables,lowrankdetails{i}.variables));
        end
    end
end

% ****************************************
% SPECIAL VARIABLES
% ****************************************
if (options.relax==1) | (options.relax==3)
    nonlins = [];
end

if (options.relax) | (options.relax==2)
    integer_variables = [];
    binary_variables  = [];
else
    integer_variables = find(ismember(used_variables,integer_variables));
    binary_variables  = find(ismember(used_variables,binary_variables));
end
parametric_variables  = find(ismember(used_variables,parametric_variables));
extended_variables =  find(ismember(used_variables,yalmip('extvariables')));

if (K.f>0) 
    % If user explicitely says remove, or user says nothing but solverdefinitions does
    % FIX : CHECK FOR NONLINEARS
    if ((options.removeequalities==1 | options.removeequalities==2) & isempty(intersect(used_variables,nonlinearvariables))) | ((options.removeequalities==0) & (solver.constraint.equalities.linear==-1))
        showprogress('Solving equalities',options.showprogress);   
        [x_equ,H,A_equ,b_equ] = solveequalities(F_struc,K,options.removeequalities==1);
        % Exit if no consistent solution exist
        if (norm(A_equ*x_equ-b_equ,'inf')>1e-5)%sqrt(eps)*size(A_equ,2))            
            diagnostic.solvertime = 0;
            diagnostic.info = yalmiperror(1,'YALMIP');
            diagnostic.problem = 1;
            solution = diagnostic;
            solution.variables = used_variables(:);
            solution.optvar = x_equ;
            % And we are done! Save the result
            sdpvar('setSolution',solution);
            return
        end
        % We dont need the rows for equalities anymore
        F_struc = F_struc(K.f+1:end,:);
        K.f = 0;
        Fold = F;
        [nlmi neq]=size(F);
        F = lmi;
        for i = 1:nlmi+neq
            if ~is(Fold(i),'equality')
                F = F + Fold(i);
            end
        end
        
        % No variables left!
        if size(H,2)==0            
            diagnostic.solvertime = 0;
            diagnostic.info = yalmiperror(0,'YALMIP');
            diagnostic.problem = 0;
            solution = diagnostic;
            solution.variables = used_variables(:);
            solution.optvar = x_equ;
            % And we are done! Save the result
            % Note, no dual is saved        
            sdpvar('setSolution',solution);
            p = checkset(F);
            if any(p<1e-5)
                diagnostic.info = yalmiperror(1,'YALMIP');
                diagnostic.problem = 1;
            end
            return
        end
        showprogress('Converting problem to new basis',options.showprogress)
        % objective in new basis
        c = H'*c;
        Q = H'*Q*H;
        % LMI in new basis
        F_struc = [F_struc*[1;x_equ] F_struc(:,2:end)*H];    
    else
        % Solver does not support equality constraints
        % and they hould be removed using double-sided
        % or user just want to remove them
         if (solver.constraint.equalities.linear==0 | options.removeequalities==-1)
            % Add equalities
            F_struc = [-F_struc(1:1:K.f,:);F_struc];
            K.l = K.l+K.f*2;
            % Keep this in mind...
            K.fold = K.f;   
            K.f = 0;
        end
        % For simpliciy we introduce a dummy coordinate change
        x_equ = 0;
        H     = 1;
    end
else
    x_equ = 0;
    H     = 1;
end    


% **************************************************
% INITIAL SOLUTION
% **************************************************
x0 = [];
if options.usex0
    if options.relax
        x0_used = relaxdouble(recover(used_variables));
    else
        x0_used = double(recover(used_variables));
    end
    x0 = zeros(sdpvar('nvars'),1);
    x0(used_variables)  = x0_used(:);   
    x0(isnan(x0))=0;
end

if ~isempty(x0)
    % Get a coordinate in the reduced space
    x0 = H\(x0(used_variables)-x_equ);
end

solveroutput = [];

% Monomial table for nonlinear variables
% FIX : Why here!!! mt handled above also
mt = yalmip('monomtable');
if size(mt,1)>size(mt,2)
    mt(size(mt,1),size(mt,1)) = 0;
end
% In local variables
mt = mt(used_variables,used_variables);
if (options.relax)|(options.relax==3)
    mt = eye(length(used_variables));
end

% Lower and upper bounds
lub = yalmip('getbounds',used_variables);
lb = lub(:,1);
ub = lub(:,2);
if 1%all(isinf([lb;ub]))
    lb = [];
    ub = [];
end
   

% ********************************
% GENERAL DATA EXCHANGE WITH SOLVER
% ********************************
interfacedata.F_struc = F_struc;
interfacedata.c = c;
interfacedata.Q = Q;
interfacedata.f = f;
interfacedata.K = K;
interfacedata.KCut = KCut;
interfacedata.integer_variables = integer_variables;
interfacedata.binary_variables  = binary_variables;
interfacedata.parametric_variables  = parametric_variables;
interfacedata.extended_variables    = extended_variables;
interfacedata.options = options;
interfacedata.solver = solver;
interfacedata.monomtable = mt;
interfacedata.x0 = x0;
interfacedata.solver = solver;
interfacedata.lb = lb;
interfacedata.ub = lb;
interfacedata.lowrankdetails = lowrankdetails;
interfacedata.solver = solver;
interfacedata.problemclass = ProblemClass;
interfacedata.getsolvertime = 1;

% ********************************
% GENERAL DATA EXCANGE SOLUTION
% ********************************
recoverdata.H = H;
recoverdata.x_equ = x_equ;
recoverdata.used_variables = used_variables;

function yesno = warningon

s = warning;
if isa(s,'char')
    yesno = isequal(s,'on');
else
    yesno = isequal(s(1).state,'on');
end