solvesdp.m

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

function diagnostic = solvesdp(varargin)
%SOLVESDP Computes solution to optimization problem
%
%   DIAGNOSTIC = SOLVESDP(F,h,options) is the common command to
%   solve optimization problems of the following kind
%
%    min        h
%    subject to
%            F >(=) 0
%
%   NOTES
%    Desptite the name, SOLVESDP is the interface for solving all
%    supported problem classes (LP, QP, SOCP, SDP, BMI, MILP,MIQP,...)
%
%    To obtain solution for a variable, use DOUBLE.
%
%    To obtain dual variable for a constraint, use DUAL.
%
%    See YALMIPERROR for error codes returned in output.
%
%   OUTPUT
%     diagnostic : Diagnostic information
%
%   INPUT
%     F          : SET object describing the constraints. Can be [].
%     h          : SDPVAR object describing the objective h(x). Can be [].
%     options    : Options structure. See SDPSETTINGS. Can be [].
%
%   EXAMPLE
%    x = sdpvar(5,1);A = randn(15,5);b = rand(15,1)*5;c = randn(5,1);
%    solvesdp(set(A*x<b),c'*x);double(x)
%
%   See also @SDPVAR/SET, DUAL, @SDPVAR/DOUBLE, SDPSETTINGS, YALMIPERROR

% Author Johan L鰂berg
% $Id: solvesdp.m,v 1.40 2005/06/08 13:52:18 joloef Exp $

yalmiptime = clock; % Let us see how much time we spend

% Arrrgh, new format with logdet much better, but we have to
% take care of old code, requires some testing...
varargin = combatible({varargin{:}});
nargin = length(varargin);
% *********************************
% CHECK INPUT
% *********************************
if nargin<1
    help solvesdp
    return
else
    F = varargin{1};
    if ~(isempty(F) | isa(F,'lmi'))
        error('First argument (F) should be an lmi object.');
    end
end

if nargin>=2
    h = varargin{2};
    if isa(h,'double')
        h = [];
    end
    if ~(isempty(h) | isa(h,'sdpvar') | isa(h,'logdet'))
        error('Second argument (the objective function h) should be an sdpvar or logdet object (or empty).');
    end
    if isa(h,'logdet')
        P  = getP(h);
        g  = getgain(h);
        if g>0
            warning('Perhaps you mean -logdet(P)...')
            diagnostic.yalmiptime = etime(clock,yalmiptime);
            diagnostic.solvertime = 0;
            diagnostic.info = yalmiperror(-2,'YALMIP');
            diagnostic.problem = -2;
            return
        end
        h = getcx(h);
        G = lmi(P>0);
        if isempty(F)
            F = G;
        else
            F = F + G;
        end
    else
        P = [];
        G = [];
    end
else
    P = [];
    h = [];
    G = [];
end

if nargin>=3
    options = varargin{3};
    if ~(isempty(options) | isa(options,'struct'))
        error('Third argument (options) should be an sdpsettings struct (or empty).');
    end
    if isempty(options)
        options = sdpsettings;
    end
else
    options = sdpsettings;
end
options.solver = lower(options.solver);

% Just for safety
if isempty(F) & isempty(G)
    F = lmi;
end

if any(is(F,'sos'))
    error('You have SOS constraints. Perhaps you meant to call SOLVESOS.');
end

% ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
% DID WE SELECT THE MOMENT SOLVER
% ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
if isequal(options.solver,'moment')
    options.solver = options.moment.solver;
    [diagnostic,x,momentdata] = solvemoment(F,h,options,options.moment.order);
    diagnostic.momentdata = momentdata;
    diagnostic.xoptimal = x;
    return
end

% ******************************************
% COMPILE IN GENERALIZED YALMIP FORMAT
% ******************************************
[interfacedata,recoverdata,solver,diagnostic,F] = compileinterfacedata(F,G,P,h,options,0);

% ******************************************
% FAILURE?
% ******************************************
if ~isempty(diagnostic)
    diagnostic.yalmiptime = etime(clock,yalmiptime);
    return
end

% ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
% DID WE SELECT THE LMILAB SOLVER WITH A KYP
% ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
if  strcmpi(solver.tag,'lmilab') & any(is(F,'kyp'))
    [diagnostic,failed] = calllmilabstructure(F,h,options);
    if ~failed % Did this problem pass (otherwise solve using unstructured call)
        diagnostic.yalmiptime = etime(clock,yalmiptime)-diagnostic.solvertime;
        return
    end
end

% ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
% DID WE SELECT THE KYPD SOLVER
% ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
if strcmpi(solver.tag,'kypd')
    diagnostic = callkypd(F,h,options);
    diagnostic.yalmiptime = etime(clock,yalmiptime)-diagnostic.solvertime;
    return
end

% ******************************************
% DID WE SELECT THE BMILIN SOLVER (obsolete)
% ******************************************
if strcmpi(solver.tag,'bmilin')
    diagnostic = callbmilin(F,h,options);
    return
end

% ******************************************
% DID WE SELECT THE BMIALT SOLVER (obsolete)
% ******************************************
if strcmp(solver.tag,'bmialt')
    diagnostic = callbmialt(F,h,options);
    return
end

%******************************************
% DID WE SELECT THE MPT solver (backwards comb)
%******************************************
if strcmpi(solver.tag,'mpt') | strcmp(solver.tag,'mpcvx')
    actually_save_output = interfacedata.options.savesolveroutput;
    interfacedata.options.savesolveroutput = 1;
    if isempty(interfacedata.parametric_variables)
        if (nargin < 4 | ~isa(varargin{4},'sdpvar'))
            error('You must specify parametric variables.')
        else
            interfacedata.parametric_variables = find(ismember(recoverdata.used_variables,getvariables(varargin{4})));
        end
    end
end

% ********************************
% TRY TO SOLVE PROBLEM
% ********************************
try
    eval(['output = ' solver.call '(interfacedata);']);
catch
    output.Primal = zeros(length(interfacedata.c),1)+NaN;
    output.Dual  = [];
    output.Slack = [];    
    output.solvertime   = nan;
    output.solverinput  = [];
    output.solveroutput = [];
    output.problem = 9;
    output.infostr = yalmiperror(output.problem,lasterr);
end

if options.dimacs
    try       
        b = -interfacedata.c;
        c = interfacedata.F_struc(:,1);
        A = -interfacedata.F_struc(:,2:end)';
        x = output.Dual;
        y = output.Primal;  
        % FIX this nonlinear crap (return variable type in
        % compileinterfacedata)
        if options.relax == 0 & any(full(sum(interfacedata.monomtable,2)~=0))
            if ~isempty(find(sum(interfacedata.monomtable | interfacedata.monomtable,2)>1))
                for i = 1:size(interfacedata.monomtable,1)
                    z(i,1) = prod(y(find(interfacedata.monomtable(i,:))));
                end
                y = z;
            end
        end
        
        if isfield(output,'Slack')
            s = output.Slack;
        else
            s = [];
        end
                    
        dimacs = computedimacs(b,c,A,x,y,s,interfacedata.K);
    catch
        dimacs = [nan nan nan nan nan nan];
    end
else
    dimacs = [nan nan nan nan nan nan];
end


% ********************************
% ORIGINAL COORDINATES
% ********************************
output.Primal = recoverdata.x_equ+recoverdata.H*output.Primal;

% ********************************
% OUTPUT
% ********************************
diagnostic.yalmiptime = etime(clock,yalmiptime)-output.solvertime;
diagnostic.solvertime = output.solvertime;
diagnostic.info = yalmiperror(output.problem,solver.tag);
diagnostic.problem = output.problem;
diagnostic.dimacs = dimacs;

% Some more info is saved internally
solution_internal = diagnostic;
solution_internal.variables = recoverdata.used_variables(:);
solution_internal.optvar = output.Primal;

if ~isempty(interfacedata.parametric_variables)
    diagnostic.mpsol = output.solveroutput;
    options.savesolveroutput = actually_save_output;
end;

if interfacedata.options.savesolveroutput
    diagnostic.solveroutput = output.solveroutput;
end
if interfacedata.options.savesolverinput
    diagnostic.solverinput = output.solverinput;
end

if options.warning & warningon & isempty(findstr(output.infostr,'No problems detected'))
    disp(['Warning: ' output.infostr]);
end

if ismember(output.problem,options.beeponproblem)
    try
        beep; % does not exist on all ML versions
    catch
    end
end

% And we are done! Save the result
if ~isempty(output.Primal)
    yalmip('setsolution',solution_internal);
end
if interfacedata.options.saveduals & solver.dual
    setduals(F,output.Dual,interfacedata.K);
end



function newinputformat = combatible(varargin)

varargin = varargin{1};

classification = 0;
% 0 : Ambigious
% 1 : Old
% 2 : New

% Try some fast methods to determine...
m = length(varargin);
if m==1
    classification = 2;
elseif m>=3 & isstruct(varargin{3})
    classification = 2;
elseif m>=4 & isstruct(varargin{4})
    classification = 1;
elseif m>=2 & isa(varargin{2},'lmi')
    classification = 1;
elseif m>=3 & isa(varargin{3},'sdpvar')
    classification = 1;
elseif m>=2 & isa(varargin{2},'sdpvar') & min(size(varargin{2}))==1
    classification = 2;
elseif m>=2 & isa(varargin{2},'sdpvar') & prod(size(varargin{2}))>=1
    classification = 1;
elseif m>=2 & isa(varargin{2},'logdet')
    classification = 2;
elseif m==2 & isempty(varargin{2})
    classification = 2;
end

if classification==0
    warning('I might have interpreted this problem wrong due to the new input format in version 3. To get rid of this warning, use an options structure');
    classification = 2;
end

if classification==2
    newinputformat = varargin;
else
    newinputformat = varargin;
    P = varargin{2};
    % 99.9% of the cases....
    if isempty(P)
        newinputformat = {newinputformat{[1 3:end]}};
    else
        if isa(P,'lmi')
            P = sdpvar(P);
        end
        if m>=3
            cxP = newinputformat{3}-logdet(P);
            newinputformat{3}=cxP;
        else
            cxP = -logdet(P);
            newinputformat{3}=cxP;
        end
        newinputformat = {newinputformat{[1 3:end]}};
    end
end

function yesno = warningon

s = warning;
yesno = isequal(s,'on');