sdmsol.m

这是matlab解2阶锥工具包

function [quiz,sdmdata] = sdmsol(quiz,pars,Rayon)
% SDMPB/SDMSOL - solve a linear matrix problem with SeDuMi
%  
% quiz = sdmsol(quiz,pars,Radius);
%  
% solve an optimization problem defined by the SDMPB object <quiz>
% <pars> containts optional parameters on the solver
%      .alg    : algorithm (1: NL cones, 2: linear) 
%      .beta   : in (0 1), neighborhood parameter (default 0.5)
%      .eps    : precision (default 1e-9)
%      .fid    : 0 to suppress SeDuMi screen output (default 1)
%      .maxiter: max number of iterations (default 50)
%      .numtol : tolerance to numerical problems (default 5e-7)   
%      .theta  : in (0 1], wide region parameter (default 0.25) 
% <Radius> : feasibility radius on the decision variables
%        0 > ... > Inf : over-bounding on the norm of the decision variables
%        = 0 (default) : a fictive radius constraint is added for conditioning
%        = Inf         : no radius and no fictive conditioning constraint
% 
% Once SDMSOL has run on the optimization problem 
% <quiz> contains the solution 
%  
% [quiz,sdmdata]=sdmsol(quiz,pars,Radius)
%  
% <sdmdata> contains the data used as inputs for the SeDuMi solver
%  
% SEE ALSO sedumi, sdmpb/get, sdmvar, sdmlmi, sdmineq and sdmobj
  
%   This file is part of SeDuMi Interface 1.04 (JUL2002)
%   Last Update 6th September 2002
%   Copyright (C) 2002 Dimitri Peaucelle & Krysten Taitz
%   LAAS-CNRS, Toulouse, France   
%  
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
% 
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
% 
%   You should have received a copy of the GNU General Public License
%   along with this program; if not, write to the Free Software
%   Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

  %%% check inputs
  if nargin<1
    error('not enough input arguments');
  elseif nargin>3
    error('too many input arguments');
  elseif ~isa(quiz,'sdmpb')
    error('1st input argument must be a SDMPB object');
  end;        
  %%% build inputs for SeDuMi
  %%% eliminate extra variables due to the structure of matrices 
  %%% (ex : symmetric terms)
  bt=quiz.obj.bt*quiz.var.struc;
  btconj=quiz.obj.btconj*quiz.var.strucconj;
  bt=bt+btconj;
  iAt=quiz.ineq.At*quiz.var.struc;
  iAtconj=quiz.ineq.Atconj*quiz.var.strucconj;
  iAt=iAt+iAtconj;
  eAt=quiz.eq.At*quiz.var.struc+quiz.eq.Atconj*quiz.var.strucconj; %Terms
                                                                   %depending on y
  eAtconj=quiz.eq.At*quiz.var.strucconj+quiz.eq.Atconj*quiz.var.struc;%Terms depending on conj(y)
  ic=quiz.ineq.c;
  ec=quiz.eq.c;
  %%% add as much lines as there are equalities involving both variables
  %%% and conjugate variables (in these cases it is necessary to write
  %%% the contraints on both the real part and the imaginary part)
  % first eliminate trivial equations
  ecAtAtconj=[ec,eAt,eAtconj];
  nonzerorows=find(max(abs(ecAtAtconj),[],2));
  ec=ec(nonzerorows,:);
  eAt=eAt(nonzerorows,:);
  eAtconj=eAtconj(nonzerorows,:);
  % initialize the fields 
  xcomplex=[];
  jj=0;
  ecnew=[];
  eAtnew=[];
  % then for each line test if it depends (or may be written as if) only
  % on the decision variables and not on their conjugate. In the last
  % case it is necessary to solve two real valued constraints on the real
  % and imaginary part
  for ii=1:length(ec)
    if ~any(eAtconj(ii,:))
      ecnew=[ecnew;ec(ii,:)];
      eAtnew=[eAtnew;eAt(ii,:)];
      if (~isreal(ec(ii,:)) | ~isreal(eAt(ii,:)))                      
        xcomplex=[xcomplex,length(ecnew)];
  %if ec and eAt are real,we look at the columns corresponding to the complex variables    
      elseif norm(eAt(ii,quiz.K.ycomplex),1)
	xcomplex=[xcomplex,length(ecnew)];
      end;
    elseif ~any(eAt(ii,:))
      ecnew=[ecnew;conj(ec(ii,:))];
      eAtnew=[eAtnew;conj(eAtconj(ii,:))];
      xcomplex=[xcomplex,length(ecnew)];
    else
      ecnew=[ecnew;ec(ii,:)];
      eAtnew=[eAtnew;eAt(ii,:)+conj(eAtconj(ii,:))];
      ecnew=[ecnew;-i*ec(ii,:)];
      eAtnew=[eAtnew;-i*eAt(ii,:)+i*conj(eAtconj(ii,:))];
    end;
  end;
  ec=ecnew;eAt=eAtnew;
  % eliminate trivial equations again (if appeared)
  ecAt=[ec,eAt];
  nonzerorows=find(max(abs(ecAt),[],2));
  ec=ec(nonzerorows,:);
  eAt=eAt(nonzerorows,:);
  %%% build the cone field
  K=[];
  K.scomplex=[];
  K.f=length(ec);
  if any(xcomplex)
    K.xcomplex=xcomplex;
  end;
  nonremoved = find(quiz.K.s);
  K.s=quiz.K.s(nonremoved);
  for ii=1:quiz.ineq.nb
    im=quiz.ineq.m(ii);
    iM=quiz.ineq.M(ii);
    if (~isreal(ic(im:iM,:)) | ~isreal(iAt(im:iM,:)))                      
      K.scomplex=[K.scomplex,ii];
  %if ic(LMI_m:LMI_M) and iAt(LMI_m:LMI_M) are real,we look at the columns corresponding to the complex variables    
    elseif norm(iAt(im:iM,quiz.K.ycomplex),1)
      K.scomplex=[K.scomplex,ii];
    end;
  end;
  K.q=[];
  K.ycomplex=quiz.K.ycomplex;
  %%% add a radius constraint on the variables
  if nargin<3
    Rayon=0;
  elseif ( ~isnumeric(Rayon) | any(size(Rayon)-[1,1]) | Rayon<0 )
    error(['3rd input argument (radius on decision variables) must be a positive scalar']);
  end;
  if Rayon>0
    if Rayon~=Inf
      nbdecvar=size(iAt,2);
      c=[ec;Rayon;sparse(nbdecvar,1);ic];
      At=[eAt;sparse(1,nbdecvar);speye(nbdecvar);iAt];
      K.q=nbdecvar+1;
    else  
      c=[ec;ic];
      At=[eAt;iAt];
    end;
  else
    nbdecvar=size(iAt,2);
    c=[ec;sparse(nbdecvar+1,1);ic];
    At=[eAt,sparse(size(eAt,1),1);
	sparse(1,nbdecvar),1;
	speye(nbdecvar),sparse(nbdecvar,1);
	iAt,sparse(size(iAt,1),1)];
    bt=sparse([bt,0]);  
    K.q=nbdecvar+1;
  end;
  
  %%% Default set default SeDuMi options
  if nargin<2
    pars.eps=1e-9;
    pars.w=[1 1e4];
  elseif isempty(pars)
    pars.eps=1e-9;
    pars.w=[1 1e4];    
  elseif ~isa(pars,'struct')
    error('2nd input argument must be a structure for SeDuMi options');
  end;
  if ~isfield(pars,'eps')
    pars.eps=1e-9;
  end;
  if ~isfield(pars,'w')
    pars.w=[1 1e4];
  end;  
  if ~isfield(pars,'fid')
    pars.fid=1;
  end;
  %%% solve the optimization problem with SeDuMi
  sdmdata.At=At;
  sdmdata.bt=bt;
  sdmdata.c=c;
  sdmdata.K=K;
  [x,y,info]=sedumi(At,bt,c,K,pars);
  %%%%% get solution (build matrices)
  if Rayon<=0
    ywR=y(1:(size(y,1)-1),1);
  else
    ywR=y;
  end;
  quiz.var.ysol=ywR;
  ytot=quiz.var.struc*quiz.var.ysol;
  ytotconj=quiz.var.strucconj*quiz.var.ysol;
  %%%%% get minimal eigenvalue of each LMI constraint
  if quiz.ineq.nb
    iz=quiz.ineq.c-quiz.ineq.At*ytot-quiz.ineq.Atconj*ytotconj;
    for ii=1:quiz.ineq.nb
      Z=full(mat(iz(quiz.ineq.m(ii):quiz.ineq.M(ii))));
      vapi=eig(Z+Z')/2;
      mvapi=min(real(vapi));
      if isempty(mvapi)
	mvapi=0;
      end;
      if (mvapi & abs(mvapi)<pars.eps)
	if mvapi>0
	  quiz.ineq.meig{ii}='eps';
	else
	  quiz.ineq.meig{ii}='-eps';
	end;
      else
	quiz.ineq.meig{ii}=mvapi;
      end;
    end;
  end;
  %%%%% get minimal norm of each LME constraint
  if quiz.eq.nb
    ez=quiz.eq.c-quiz.eq.At*(ytot+conj(ytotconj))-quiz.eq.Atconj*(conj(ytot)+ytotconj);
    for ii=1:quiz.eq.nb
      Z=full(sdmmat(ez(quiz.eq.m(ii):quiz.eq.M(ii)),quiz.eq.row(ii)));
      mnorm=norm(Z,1);
      if isempty(mnorm)
	mnorm=0;
      end;
      if (mnorm & mnorm<pars.eps)
	quiz.eq.mnorm{ii}='eps';
      else
	quiz.eq.mnorm{ii}=mnorm;
	quiz.feas;
      end;
    end;
  else 
    ez=0;
  end;
  %%%%% get the optimum
  quiz.obj.opt=quiz.obj.bt*ytot+quiz.obj.btconj*ytotconj;
  %%% result analysis 
  quiz.info=info;
  z=c-At*y;
  z=z((size(eAt,1)+1):size(z,1));
  allLMIfeas=min(real(eigK(z,K)));
  allLMEfeas=max(abs(ez));
  if ( allLMIfeas>=0 & allLMEfeas<=pars.eps )
    quiz.feas='feasible';
  elseif ( allLMIfeas>=-pars.eps & allLMEfeas<=pars.eps & abs(info.feasratio-1)<1e-1 )
    quiz.feas='feasible';
  elseif (abs(info.feasratio+1)<1e-1 | info.pinf | info.dinf)
    quiz.feas='infeasible';
  else
    quiz.feas='marginal feasible';
  end;
  if pars.fid
    fprintf('%s',quiz.feas);
    if quiz.feas(1)=='m'
      if allLMEfeas>pars.eps
        fprintf('  :  %g  should be <= %g',allLMEfeas,pars.eps);
      end;
      if allLMIfeas < -pars.eps
	fprintf('  :  %g  should be >= 0',allLMIfeas);
      end;
    end;
    fprintf('\n\n');
  end;