convexitypropagation.m

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

function [F,failure] = convexitypropagation(F,h)

% (not used) Checkfor = 1  : Propagate convexity
%                       -1 : Propagate concavity

% Author Johan L鰂berg
% $Id: convexitypropagation.m,v 1.16 2005/02/20 17:04:28 johanl Exp $

% FIX : Current code very experimental, un-necessarily complex and stupidly
% conservative

% All variables used in problem
allvars = unique(union(depends(h),depends(F)));
allextvars = yalmip('extvariables');

% Index to variables modeling operators
extended = find(ismembc(allvars,allextvars));

% Simple check for nonlinearity
% FIX : Robustify for non-integer and negative
mt = yalmip('monomtable');

% No extended variables at all
if isempty(extended)
    failure = 0;
    return
end

% % Bilinear or similiar
% if any (sum(mt(allvars,allextvars)~=0,2)>1)
%     failure = 1;
%     F = set([]);
%     return;
% end

% All analyzed variables
propagated = [];

% Run recursively by extending variables defining nonlinear operators
failure = 0;loop = 1;start = 1;
is_elementwise = is(F,'element-wise');
while ~isempty(extended)
    
    % ******************************************************
    % Create all candidate constraints
    % Example : "t models max(x,y)" generates the constraint
    %            set(x<t)+set(y<t)
    % ******************************************************
    [F_extend,convexconcave,defining_sets] = model(recover(allvars(extended)));
           
    ok  = 1;
    
%      % First, check for simple polynomial stuff, like norm^2    
%      nonlinearuse = find(sum(mt(allvars,allvars(extended)),2)>1);
%      if ~isempty(nonlinearuse)
%          for i = 1:length(extended)
%              powers = mt(allvars,allvars(extended(i)));
%              if isempty(setdiff(powers,[0 1]))
%                  
%              else
%                  if (convexconcave{i}(1)==1) & (convexconcave{i}(3)==1) & all(powers(powers~=0)>=1)                 
%                  else
%                      ok = 0;
%                  end
%              end
%          end                  
%      end

  
    % ******************************************************
    % Major computation goes here for convexity checks
    %
    % Go through all constraint to see that convexity w.r.t
    % these implied variables is ok
    %
    % We need to check that convex variables occur on the left
    % side of <, and that it is not used in any cones or equalities.
    %
    % Of-course, this is very conservative, but it is sufficient
    % for a lot of practical cases.
    % The reason for the simple approach below is due to the lack
    % of monotonocity info in the operators. Probable improvement later...
    % With this info, we would just check top layer constraints, and
    % then propagate the involved variables, using composition rules.
    %
    % The operators may not be used in anything but element-wise
    % constraints.
    %
    % When we add new constraints to model a variable, these
    % constraints may contain other operator variables (from
    % max(abs(x),max(z,t)) and similiar constructions.
    %
    % However, in these constraints, we should not check the rules
    % with respect to the variable that we are checking (for the model
    % "t models max(z,y)", we get z<t,y<t. t is on right-hand side here)
    % ******************************************************
     
    % Original constraints
    for j = start:length(F)      
        Fj = sdpvar(F(j));
        if any(ismembc(getvariables(Fj),allvars(extended)))%%~isempty(intersect(getvariables(Fj),allvars(extended)))
            %if ~isempty(intersect(getvariables(Fj),allvars(extended)))            
            for i = 1:length(extended)
                if (convexconcave{i}(1)~=0)
                    extvar = allvars(extended(i));
                    this_base = getbasematrix(Fj,extvar);
                    if nnz(this_base)>0
                        if is_elementwise(j)
                            switch convexconcave{i}(1)
                            case 0
                                % don't bother
                            case -1 % Concave
                                ok = ok & (all(this_base>=0)) & is_elementwise(j);
                            case 1  % Convex
                                ok = ok & (all(this_base<=0)) & is_elementwise(j);
                            otherwise
                                error('Nonlinear operator defined as neither convex (1) or concave (-1)');
                            end
                        else
                            ok = 0;
                        end
                    end
                end
            end
        end
    end
    
    % Constraints generated when expanding the model
    for j = 1:length(F_extend);
        Fj = sdpvar(F_extend(j));
        for i = 1:length(extended)
            if (convexconcave{i}(1)~=0)
                extvar = allvars(extended(i));
                if ~ismembc(j,defining_sets{i}) % don't check the sets that we generated to model this variable!
                    this_base = getbasematrix(Fj,extvar);
                    if nnz(this_base)>0
                        switch convexconcave{i}(1)
                        case 0
                            % don't bother
                        case -1 % Concave
                            ok = ok & (all(this_base>=0));
                        case 1  % Convex
                            ok = ok & (all(this_base<=0));
                        otherwise
                            error('Nonlinear operator defined as neither convex (1) or concave (-1)');
                        end
                    end
                end
            end
        end
    end
    
    % Check objective function (only needed in first pass)
    if loop == 1
        for i = 1:length(extended)
            if (convexconcave{i}(1)~=0)                
                extvar = allvars(extended(i));
                this_base = getbasematrix(h,extvar);
                switch convexconcave{i}(1)
                case 0
                    % don't bother
                case 1 % Convex
                    ok = ok & (all(this_base>=0));
                case -1 % Concave
                    ok = ok & (all(this_base<=0));
                otherwise
                    error('Nonlinear operator defined as neither convex (1) or concave (-1)');
                end
            end
        end
    end
    
    if ok
        start = start + length(F); % We only need to check the new constraints in the next pass
        loop = loop + 1;
        F = F + F_extend;        
        % We have now extended the model, might have genererated
        % new variables. Iterate!
        propagated = unique([propagated allvars(extended)]);
        allvars = depends(F_extend); % Only new variables here
        %newextended = find(ismember(allvars,yalmip('extvariables')));
        newextended = find(ismembc(allvars,yalmip('extvariables')));
        %extended = newextended(find(~ismember(allvars(newextended),propagated)));
        extended = newextended(find(~ismembc(allvars(newextended),propagated)));
        is_elementwise = is(F,'element-wise');
    else
        failure = 1;
        return
    end   
end