cutsdp.m

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

function output = cutsdp(p)
% CUTSDP
%
% See also SOLVESDP, BNB, BINVAR, INTVAR, BINARY, INTEGER, LMI

% Author Johan L鰂berg 
% $Id: cutsdp.m,v 1.6 2005/05/10 15:04:52 joloef Exp $

% ********************************
%% INITIALIZE DIAGNOSTICS IN YALMIP
% ********************************
bnbsolvertime = clock; 
showprogress('Cutting plane solver started',p.options.showprogress);

% ********************************
%% If we want duals, we may not 
%  extract bounds.
% However, bounds must be extracted
% in discrete problems.
% ********************************
if p.options.cutsdp.recoverdual
    warning('Dual recovery not implemented yet in CUTSDP')
end
if 1%~isempty(p.binary_variables) | ~isempty(p.integer_variables) | (p.options.cutsdp.recoverdual==0)

    p.options.cutsdp.recoverdual=0;
    % ********************************
    %% Define infinite bounds
    % ********************************
    if isempty(p.ub)
        p.ub = repmat(inf,length(p.c),1);
    end
    if isempty(p.lb)
        p.lb = repmat(-inf,length(p.c),1);
    end
    
    % ********************************
    %% Extract better bounds from model
    % ********************************
    if ~isempty(p.F_struc)
        [lb,ub,used_rows] = findulb(p.F_struc,p.K);
        if ~isempty(used_rows)
            lower_defined = find(~isinf(lb));
            if ~isempty(lower_defined)
                p.lb(lower_defined) = max(p.lb(lower_defined),lb(lower_defined));
            end
            upper_defined = find(~isinf(ub));
            if ~isempty(upper_defined)
                p.ub(upper_defined) = min(p.ub(upper_defined),ub(upper_defined));
            end
            p.F_struc(p.K.f+used_rows,:)=[];
            p.K.l = p.K.l - length(used_rows);
        end
    end

    % ********************************
    %% ADD CONSTRAINTS 0<x<1 FOR BINARY
    % ********************************
    if ~isempty(p.binary_variables)
        p.ub(p.binary_variables) =  min(p.ub(p.binary_variables),1);
        p.lb(p.binary_variables) =  max(p.lb(p.binary_variables),0);
    end
    
    p.ub = min(p.ub,p.options.cutsdp.variablebound');
    p.lb = max(p.lb,-p.options.cutsdp.variablebound');
    
    % *******************************
    %% PRE-SOLVE (nothing fancy coded)
    % *******************************
    if isempty(find(isinf([p.ub;p.lb]))) & p.K.l>0 
        [p.lb,p.ub] = tightenbounds(-p.F_struc(1+p.K.f:p.K.f+p.K.l,2:end),p.F_struc(1+p.K.f:p.K.f+p.K.l,1),p.lb,p.ub,p.integer_variables);  
    end
end

% *******************************
%% PERTURBATION OF LINEAR COST
% *******************************
p.corig = p.c;
if nnz(p.Q)~=0
    g = randn('seed');
    randn('state',1253); %For my testing, I keep this the same...   
    % This perturbation has to be better. Crucial for many real LP problems
    p.c = (p.c).*(1+randn(length(p.c),1)*1e-4); 
    randn('seed',g);
end

% *******************************
%% We don't need this
% *******************************
p.options.savesolverinput  = 0;
p.options.savesolveroutput = 0;

% *******************************
%% Display logics
% 0 : Silent
% 1 : Display cut progress
% 2 : Display node solver prints
% *******************************
switch max(min(p.options.verbose,3),0)
case 0
    p.options.cutsdp.verbose = 0;
case 1
    p.options.cutsdp.verbose = 1;
    p.options.verbose = 0;
case 2
    p.options.cutsdp.verbose = 2;
    p.options.verbose = 0;
case 3
    p.options.cutsdp.verbose = 2;
    p.options.verbose = 1;
otherwise 
    p.options.cutsdp.verbose = 0;
    p.options.verbose = 0;    
end

% *******************************
%% START CUTTING
% *******************************
[x_min,solved_nodes,lower,feasible,D_struc] = cutting(p);
%% --

% **********************************
%% CREATE SOLUTION
% **********************************
output.problem = 0;
if ~feasible
    output.problem = 1;
end
if solved_nodes == p.options.cutsdp.maxiter
    output.problem = 3;
end
output.solved_nodes = solved_nodes;
output.Primal       = x_min;
output.Dual = D_struc;
output.Slack = [];
output.solverinput  = 0;
output.solveroutput =[];
output.solvertime   = etime(clock,bnbsolvertime);
%% --

function [x,solved_nodes,lower,feasible,D_struc] = cutting(p)

% *******************************
%% Sanity check
% *******************************
if any(p.lb>p.ub)
    x = zeros(length(p.c),1);
    solved_nodes = 0;
    lower = inf;
    feasible = 0;
    D_struc = [];
    return
end

% *******************************
%% Create function handle to solver
% *******************************
cutsolver = p.solver.cutsolver.call;    

% *******************************
%% Create copy of model without
%  the SDP part
% *******************************
p_lp = p;
p_lp.F_struc = p_lp.F_struc(1:p.K.l+p.K.f,:);
%p_lp.F_struc(end-sum((p.K.s).^2)+1:end,:) = [];
p_lp.K.s = 0;    
     
% *******************************
%% DISPLAY HEADER
% *******************************
if p.options.cutsdp.verbose    
    disp('* Starting YALMIP cutting plane for MISDP based on MILP'); 
    disp(['* Lower solver   : ' p.solver.cutsolver.tag]);    
    disp(['* Max iterations : ' num2str(p.options.cutsdp.maxiter)]);        
end

if p.options.bnb.verbose;            disp(' Node       Infeasibility.     Lower    LP cuts');end;

%% Initialize diagnostic
infeasibility = -inf;
solved_nodes = 0;
feasible = 1;
lower = -inf;
saveduals = 1;

%% Add diagonal cuts to begin with
savedCuts = [];
savedIndicies = [];
if p.K.s(1)>0         
    top = p.K.f+p.K.l+1;    
    for i = 1:length(p.K.s)        
        n = p.K.s(i);
        newF=[];
        for m = 1:p.K.s(i)
            d = eyev(p.K.s(i),m);
          %  dd = d*d';dd=dd(:);dd = find(dd);
            index = (1+(m-1)*(p.K.s(i)+1));
            newF = [newF;p.F_struc(top+index-1,:)];
          %  for j = 1:size(p.F_struc,2)
          %      newF(m,j)= d'*reshape(p.F_struc(top:top+n^2-1,j),n,n)*d;
          %  end            
        end
        % Clean
        newF(abs(newF)<1e-12) = 0;
        keep=find(any(newF(:,2:end),2));
        newF = newF(keep,:);        
        % Save info for recovering duals in SDP problems
        if p.options.cutsdp.recoverdual
            savedCuts{i}=eye(n);savedCuts{i}=savedCuts{i}(:,keep);;
            savedIndicies{i}=(size(p_lp.F_struc,1)+1:1:size(p_lp.F_struc,1)+size(newF,1))';
        end
        p_lp.F_struc = [p_lp.F_struc;newF];
        p_lp.K.l = p_lp.K.l + size(newF,1);               
        top = top+n^2;
    end
end
  
goon = 1;
while goon
    % Solve SDP relaxation  
    output = feval(cutsolver,p_lp);
    
    if p.options.cutsdp.recoverdual
        saver.s=savedCuts;
        saver.i=savedIndicies;
    end
    
    if output.problem == 1
        % LP relaxation was infeasible, hence problem is infeasible
        feasible = 0;
        lower = inf;
        goon = 0;
        x = zeros(length(p.c),1);
        lower = inf;        
    else
        % Relaxed solution
        x = output.Primal;
        lower = p.f+p.c'*x+x'*p.Q*x;
        
        infeasibility = 0;
        if p.K.s(1)>0
            % Add cuts
            top = p.K.f+p.K.l+1;           
            for i = 1:1:length(p.K.s)
                n = p.K.s(i);
                X = p.F_struc(top:top+n^2-1,:)*[1;x];
                X = reshape(X,n,n);
                newcuts = 1;
%               newF = zeros(n,size(p.F_struc,2));
                [d,v] = eig(X);
                infeasibility = min(infeasibility,min(diag(v)));
                dummy=[];
                newF = [];
                if infeasibility<0
                    [ii,jj] = sort(diag(v));
                    for m = jj(1:min(length(jj),p.options.cutsdp.cutlimit))'%find(diag(v<0))%1:1%length(v)
                        if v(m,m)<0
                            dummy = [dummy d(:,m)];
 %                           d(:,m)'*kron(d(:,m),speye(n))'*p.F_struc(top:top+n^2-1,:);
                            newF = [newF;d(:,m)'*(kron(d(:,m),speye(n))'*p.F_struc(top:top+n^2-1,:))];
%                            for j = 1:length(x)+1;
%                                S = reshape(p.F_struc(top:top+n^2-1,j),n,n);
%                                newF(newcuts,j)= d(:,m)'*S*d(:,m);                                
%                            end
                            newcuts = newcuts + 1;
                        end
                    end
                end
                newF(abs(newF)<1e-12) = 0;
                keep=find(any(newF(:,2:end),2));
                newF = newF(keep,:);
                if size(newF,1)>0
                    if p.options.cutsdp.recoverdual
                        savedCuts{i} = [savedCuts{i} dummy(:,keep)];
                        savedIndicies{i}=[savedIndicies{i};((size(p_lp.F_struc,1)+1):1:(size(p_lp.F_struc,1)+size(newF,1)))'];
                    end
                    p_lp.F_struc = [p_lp.F_struc;newF];
                    p_lp.K.l = p_lp.K.l + size(newF,1);
                   [i,j] = sort(p_lp.F_struc*[1;x]);
                   %p_lp.F_struc = p_lp.F_struc(j(1:20),:);
                   %p_lp.K.l = 20;
                end
                
                top = top+n^2;
            end           
        end
        goon = infeasibility < p.options.cutsdp.feastol;
        goon = goon & feasible;
        goon = goon & (solved_nodes < p.options.cutsdp.maxiter-1);
    end
    solved_nodes = solved_nodes + 1;
    if p.options.cutsdp.verbose;fprintf(' %4.0f :      %12.3E      %12.3E      %2.0f\n',solved_nodes,infeasibility,lower,p_lp.K.l-p.K.l);end    
end

D_struc = [];
% if feasible & p.options.cutsdp.recoverdual
%     if ~isempty(output.Dual)
%         D_struc = output.Dual(1:1:p.K.f+p.K.l);
%         if p.K.s(1)>0
%             lambda=output.Dual;
%             savedIndicies=saver.i;
%             savedCuts=saver.s;
%             for i = 1:1:length(p.K.s)
%                 X=0;
%                 for j=1:length(savedIndicies{i})
%                     X=X+lambda(savedIndicies{i}(j))*savedCuts{i}(:,j)*savedCuts{i}(:,j)';
%                 end
%                 D_struc = [D_struc;X(:)];
%             end 
%         end
%     end
% end