newconvexitypropagation.m

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

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

% Author Johan L鰂berg
% $Id: newconvexitypropagation.m,v 1.1 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));

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

% All analyzed variables
propagated = [];

% ***************************************************
% At the top level, we have 
% min sum c'x+sum(di hi(x))
%  s.t    Ax+sum(ei8i(x)) > 0
%
% All we have to check is that di>0,hi convex
%                           or di<0,hi concave
%                              ei<0,fi convex
%                              ei>0,fi concave
% ***************************************************
is_elementwise = is(F,'element-wise');

[F_extend,convexconcave,defining_sets,arguments] = model(recover(allvars(extended)));

ok=1;
% Original constraints
for j = 1:length(F)      
    Fj = sdpvar(F(j));
    if any(ismembc(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) or don''t both (0)');
                        end
                    else
                        ok = 0;
                    end
                end
            end
        end
    end
end

% Objective function
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

% OK, the problem is correct on top level
% We now need to expand the arguments in each nonlinear operator
% to see that composition ruls are allowed
% Example, f(g(x)) < 1
% We need to ensure that g(x) is convex and positive
%                        g(x) is convex and non-decreasing


if ok
    for i = 1:length(extended)
        [F_extend,ok] = checkCompositionRules(arguments{i},convexconcave{i},F_extend);
    end
end

F=F+F_extend;


    failure = ~ok;

  
function [PropagationConstraints,ok] = checkCompositionRules(arguments,fcnClass,PropagationConstraints);

ok = 1;
i  = 1;
while ok & (i<=length(arguments))
    argi = arguments{i};
    allvars = getvariables(argi);           
    allextvars = yalmip('extvariables');
    % Index to variables modeling operators
    extended = find(ismembc(allvars,allextvars));
    if ~isempty(extended)
        j = 1;
        while ok & (j <=length(extended))
            [F_extend,convexconcave,newdefining_sets,newarguments] = model(recover(allvars(extended(j))));
            PropagationConstraints=PropagationConstraints+F_extend;
            argiBase = getbase(arguments{i});argiBase=argiBase(:,2:end); 
            ok =  0;
            % f = h(g(x))
            % h convex nondecreasing, g multiplied by positive, g convex
            ok = ok | ((fcnClass(1)==1) & (argiBase(:,extended)>=0) & (convexconcave{1}(1)==1) & (fcnClass(2)==1));
            % h convex nonincreasing, g multiplied by negative, g concave
            % ok = ok | ((fcnClass(1)==1) & (argiBase(:,extended)<=0) & (convexconcave{1}(1)==-1) & (fcnClass(3)==-1));
            % h concave, g multiplied by negative, g concave
            ok = ok | ((fcnClass(1)==1) & (argiBase(:,extended)<=0) & (convexconcave{1}(1)==-1) & (fcnClass(3)==1));
            ok=all(ok);
            ok=all(ok);
            
            
%             ok =  0;
%             % f convex, multiplied by positive, g convex non-decreasing
%             ok = ok | ((fcnClass(1)==1) & (argiBase(:,extended)>=0) & (convexconcave{1}(1)==1) & (convexconcave{1}(2)==1));
%             % f convex, multiplied by positive, g convex positive
%             ok = ok | ((fcnClass(1)==1) & (argiBase(:,extended)>=0) & (convexconcave{1}(1)==1) & (convexconcave{1}(3)==1));
%             % f convex, multiplied by positive, g convex non-decreasing
%             ok = ok | ((fcnClass(1)==1) & (argiBase(:,extended)<=0) & (convexconcave{1}(1)==-1) & (convexconcave{1}(2)==1));            
%             % f convex, multiplied by negative, g convex positive
%             ok = ok | ((fcnClass(1)==1) & (argiBase(:,extended)<=0) & (convexconcave{1}(1)==-1) & (convexconcave{1}(3)==1));
%             ok=all(ok);
            if ok
                k=1;
                while ok & (k <= length(newarguments))
                    [PropagationConstraints,ok] = checkCompositionRules(newarguments{k},convexconcave{1},PropagationConstraints);
                    k=k+1;
                end
            end
            j = j + 1;
        end
    else
        % Nothing to check, its an affine expression
        ok = 1;
    end     
    i=i+1;
end