solvesos.m

求解线性矩阵不等式简单方便--与LMI工具箱相比

for constraint = 1:length(p)
    % *********************************************
    %% FIND THE VARIABLES IN p, SORT, UNIQUE ETC
    % *********************************************
    if options.verbose>1;disp(['Creating SOS-description ' num2str(constraint) '/' num2str(length(p)) ]);end

    pVariables = depends(p{constraint});
    AllVariables = uniquestripped([pVariables ParametricVariables]);
    MonomVariables = setdiff1D(pVariables,ParametricVariables);
    x = recover(MonomVariables);
    z = recover(AllVariables);
    MonomIndicies = find(ismember(AllVariables,MonomVariables));
    ParametricIndicies = find(ismember(AllVariables,ParametricVariables));

    if isempty(MonomIndicies)
        % This is the case set(sos(t)) where t is a parametric (matrix) variable
        % This used to create an error message befgore to avoid some silly
        % bug in the model generation. Creating this error message is
        % stupid, but at the same time I can not remember where the bug was
        % and I have no regression test  for this case. To avoid
        % introducing same bug again by mistake, I create all data
        % specifically for this case        
        previous_exponent_p_monoms = [];%exponent_p_monoms;
        n = length(p{constraint});
        A_basis = getbase(sdpvar(n,n,'full'));d = find(triu(ones(n)));A_basis = A_basis(d,2:end);
        BlockedA{constraint} = {A_basis};
        Blockedb{constraint} = p{constraint}(d);        
        BlockedN{constraint} = {zeros(1,0)};
        Blockedx{constraint} = x;
        Blockedvarchange{constraint}=zeros(1,0);    
        continue
        %  error('You have constraints of the type set(sos(f(parametric_variables))). Please use set(f(parametric_variables) > 0) instead')
    end

    % *********************************************
    %% Express p in monimials and coefficients  
    % *********************************************
    [exponent_p,p_base] = getexponentbase(p{constraint},z);
    
    % *********************************************
    %% Powers for user defined candidate monomials
    % (still experimental)
    % *********************************************
    if ~all(cellfun('isempty',candidateMonomials))
        exponent_c = [];
        if isa(candidateMonomials{constraint},'cell')            
            for i = 1:length(candidateMonomials{constraint})
                 exponent_c{i} = getexponentbase(candidateMonomials{constraint}{i},z);
                 exponent_c{i} = exponent_c{i}(:,MonomIndicies);
            end
        else
            exponent_c{1} = getexponentbase(candidateMonomials{constraint},z);
            exponent_c{1} = exponent_c{1}(:,MonomIndicies);
        end
    else
        exponent_c = [];
    end

    % *********************************************
    %% STUPID PROBLEM WITH ODD HIGHEST POWER?...
    % *********************************************
    if isempty(ParametricIndicies)
        max_degrees = max(exponent_p(:,MonomIndicies),[],1);
        bad_max = any(max_degrees-fix((max_degrees/2))*2);
        if bad_max
            for i = 1:length(p)
                Q{i}=[];
                m{i}=[];
            end
            residuals=[];
            everything = [];
            sol.yalmiptime = etime(clock,yalmip_time);
            sol.solvertime = 0;
            sol.info = yalmiperror(1,'YALMIP');
            sol.problem = 2;
            return
        end
    end

    % *********************************************
    %% Can we make a smart variable change (no code)
    % *********************************************
    exponent_p_monoms = exponent_p(:,MonomIndicies);
    varchange = ones(1,size(MonomIndicies,2));

    % *********************************************
    %% Unique monoms (copies due to parametric terms)
    % *********************************************
    exponent_p_monoms = uniquesafe(exponent_p_monoms,'rows');

    if options.sos.reuse & constraint > 1 & isequal(previous_exponent_p_monoms,exponent_p_monoms)
        % We don't have to do anything, candidate monomials can be-used
        if options.verbose>1;disp(['Re-using all candidate monomials (same problem structure)']);end
    else


        % *********************************************
        %% PRUNE W.R.T USER DEFINED
        %*********************************************
        %         if ~isempty(exponent_c)
        %             hash = randn(size(exponent_m{1},2),1);
        %             for i = 1:length(exponent_m)
        %                 hash = randn(size(exponent_m{i},2),1);
        %                 keep = find(ismember(exponent_m{i}*hash,exponent_c{1}*hash));
        %                 if length(keep)>0
        %                     exponent_m{i} = exponent_m{i}(keep,:);
        %                 else
        %                     exponent_m{i}=[];
        %                 end
        %             end
        %         end

        % *********************************************
        % User has supplied the whole candidate structure
        % Don't process this
        % *********************************************
        if ~isempty(exponent_c)
            exponent_m{1} = [];
            N = {};
            for i = 1:length(exponent_c)
                exponent_m{i} = [exponent_m{1};exponent_c{i}];
                N{i,1} = exponent_c{i};
            end
        else
            % *********************************************
            %% CORRELATIVE SPARSITY PATTERN
            % *********************************************
            [C,csclasses] = corrsparsity(exponent_p_monoms,options);

            % *********************************************
            %% GENERATE MONOMIALS
            % *********************************************
            exponent_m = monomialgeneration(exponent_p_monoms,csclasses);

            % *********************************************
            %% REDUCE #of MONOMIALS
            % *********************************************
            % Fix for matrix case, perform newton w.r.t
            % diagonal polynomials only. This can be
            % improved, but for now, keep it simple...
            n = length(p{constraint});diag_elements = 1:(n+1):n^2;used_diagonal = find(any(p_base(diag_elements,:),1));
            exponent_p_monoms_diag = exponent_p(used_diagonal,MonomIndicies);
            exponent_m = monomialreduction(exponent_m,exponent_p_monoms_diag,options,csclasses,LPmodel);

            % *********************************************
            %% BLOCK PARTITION THE MONOMIALS BY CONGRUENCE
            % *********************************************
            N = congruenceblocks(exponent_m,exponent_p_monoms,options,csclasses);
            % *********************************************
            %% REDUCE FURTHER BY EXPLOITING BLOCK-STRUCTURE
            % *********************************************
            N = blockmonomialreduction(exponent_p_monoms_diag,N,options);
            
        end


        % *********************************************
        %% PREPARE FOR SDP FORMULATION BY CALCULATING ALL
        % POSSIBLE MONOMIAL PRODUCS IN EACH BLOCK
        % *********************************************
        [exponent_m2,N_unique] = monomialproducts(N);

        % *********************************************
        %% CHECK FOR BUG/IDIOT PROBLEMS IN FIXED PROBLEM
        % *********************************************
        if isempty(ParametricIndicies)
            if ~isempty(setdiff(exponent_p_monoms,N_unique(:,3:end),'rows'))
                for i = 1:length(p)
                    Q{i} = [];
                    m{i} = [];
                end
                residuals = [];everything = [];
                sol.problem = 2;
                sol.info = yalmiperror(1,'YALMIP');
                warning('Problem is trivially infeasible (odd highest power?)');
                return
            end
        end
    end

    previous_exponent_p_monoms = exponent_p_monoms;

    % *********************************************
    %% GENERATE DATA FOR SDP FORMULATIONS
    % *********************************************
    p_base_parametric = [];
    n = length(p{constraint});
    for i=1:length(p{constraint})
        for j = 1:length(p{constraint})
            p_base_parametric = [p_base_parametric parameterizedbase(p{constraint}(i,j),z,params,ParametricIndicies,exponent_p,p_base((i-1)*n+j,:))];
        end
    end
    [BlockedA{constraint},Blockedb{constraint}] = generate_kernel_representation_data(N,N_unique,exponent_m2,exponent_p,p{constraint},options,p_base_parametric,ParametricIndicies,MonomIndicies,constraint==1);

    % SAVE FOR LATER
    BlockedN{constraint} = N;
    Blockedx{constraint} = x;
    Blockedvarchange{constraint}=varchange;
end

% *********************************************
%% And now get the SDP formulations
%
% The code above has generated matrices A and b
% in AQ == b(parametric)
%
% We use these to generate kernel or image models
% *********************************************
sol.problem = 0;
switch options.sos.model
    case 1
        % Kernel model
        [F,obj,BlockedQ,Primal_matrices,Free_variables] = create_kernelmodel(BlockedA,Blockedb,F_parametric,parobj,options,[]);
    case 2
        % Image model
        [F,obj,BlockedQ,sol] = create_imagemodel(BlockedA,Blockedb,F_parametric,parobj,options);

    case 3
        % Un-official model to solve bilinearly parameterized SOS using SDPLR
        [F,obj,options] = create_lrmodel(BlockedA,Blockedb,F_parametric,parobj,options,ParametricVariables);

    otherwise
end

% Unofficial fifth output with pseudo-numerical data
everything.BlockedA = BlockedA;
everything.Blockedb = Blockedb;
everything.F = F;
everything.obj = obj;

% % **********************************************
% %% SOLVE SDP
% % **********************************************
if sol.problem == 0
    if options.verbose > 0
        disp(' ');
    end
    if all(isinf(ranks))
        % Standard case
        %sol = solvesdp(F+set(-10<recover(depends(F)) < 10),obj,sdpsettings('bmibnb.maxit',200,'solver','bmibnb','bmibnb.lpred',1))
        sol =  solvesdp(F,obj,options);
    else
        % We have to alter the problem slightly if there are rank
        % constraints on the decompositions
        sol =  solveranksos(F,obj,options,ranks,BlockedQ);      
    end
end

% **********************************************
%% Recover solution (and rescale model+solution)
% **********************************************
for constraint = 1:length(p)
    for i = 1:length(BlockedQ{constraint})
        doubleQ{constraint}{i} = normp(constraint)*double(BlockedQ{constraint}{i});
    end
    doubleb{constraint}=normp(constraint)*double(Blockedb{constraint});
end

% **********************************************
%% Minor post-process
% **********************************************
if all(sizep<=1)
    [doubleQ,residuals] = postprocesssos(BlockedA,doubleb,doubleQ,[],options);
else
    residuals = 0;
end

% **********************************************
%% Safety check for bad solvers...
% **********************************************
if any(residuals > 1e-3) & (sol.problem == 0) & options.verbose>0
    disp(' ');
    disp('-> Although the solver indicates no problems,')
    disp('-> the residuals in the problem are really bad.')
    disp('-> My guess: the problem is probably infeasible.')
    disp('-> Make sure to check how well your decomposition')
    disp('-> matches your polynomial (see manual)')
    disp('-> You can also try to change the option sos.model')
    disp('-> or use another SDP solver.')
    disp(' ');
end

% **********************************************
%% Confused users. Primal dual kernel image?...
% **********************************************
if options.verbose > 0
    if sol.problem == 1
        if options.sos.model == 1
            disp(' ')
            disp('-> Solver reported infeasible dual problem.')
            disp('-> Your SOS problem is probably unbounded.')
        elseif options.sos.model == 2
            disp(' ')
            disp('-> Solver reported infeasible primal problem.')
            disp('-> Your SOS problem is probably infeasible.')
        end
    elseif sol.problem == 2
        if options.sos.model == 1
            disp(' ')
            disp('-> Solver reported unboundness of the dual problem.')
            disp('-> Your SOS problem is probably infeasible.')
        elseif options.sos.model == 2
            disp(' ')
            disp('-> Solver reported unboundness of the primal problem.')
            disp('-> Your SOS problem is probably unbounded.')            
        end
    elseif sol.problem == 12
            disp(' ')
            disp('-> Solver reported unboundness or infeasibility of the primal problem.')
            disp('-> Your SOS problem is probably unbounded.')            
    end
end

% **********************************************
%% De-block
% **********************************************
for constraint = 1:length(p)
    Qtemp = [];
    for i = 1:length(BlockedQ{constraint})
        Qtemp = blkdiag(Qtemp,doubleQ{constraint}{i});
    end
    Q{constraint} = full(Qtemp);
end

% **********************************************
%% Experimental code for declaring sparsity
% **********************************************
if options.sos.impsparse == 1
    somesmall = 0;
    for i = 1:length(BlockedQ)
        for j =  1:length(BlockedQ{i})
            small = find(abs(double(BlockedQ{i}{j}))<options.sos.sparsetol);
            nullThese{i}{j} = small;
            if ~isempty(small)
                somesmall = 1;
            end
        end
    end
    if somesmall
        [F,obj,BlockedQ,Primal_matrices,Free_variables] = create_kernelmodel(BlockedA,Blockedb,F_parametric,parobj,options,nullThese);
        sol =  solvesdp(F,obj,options);
        for constraint = 1:length(p)
            for i = 1:length(BlockedQ{constraint})
                doubleQ{constraint}{i} = normp(constraint)*double(BlockedQ{constraint}{i});
            end
            doubleb{constraint}=normp(constraint)*double(Blockedb{constraint});
        end
        % for i = 1:length(doubleQ)
        %     for j =  1:length(doubleQ{i})
        %         doubleQ{i}{j}(nullThese{i}{j}) = 0;
        %     end
        % end
        % **********************************************
        %% Post-process
        % **********************************************
        [doubleQ,residuals] = postprocesssos(BlockedA,doubleb,doubleQ,nullThese,options);
        for constraint = 1:length(p)
            Qtemp = [];
            for i = 1:length(BlockedQ{constraint})
                Qtemp = blkdiag(Qtemp,doubleQ{constraint}{i});
            end
            Q{constraint} = Qtemp;
        end
    end
end

% *********************************************
%% EXTRACT DECOMPOSITION
% *********************************************
switch sol.problem
    case {0,1,2,3,4,5} % Well, it didn't f**k up completely at-least

        %         % If SDPLR was sucessful, the last dual should have rank 1
        %         Z = dual(F(end));