sdmineq.m

这是matlab解2阶锥工具包

function quiz = sdmineq(quiz,lmidata,Xindex,inL,inR,inK,kside);  
% SDMPB/SDMINEQ - add term in some inequality constraint
%  
% quiz = sdmineq(quiz,[lmiindex blrow blcol],Xindex,L,R);
%
%   if <R> is non empty adds the term (L*X*R)  to the (blrow,blcol) block
%    and the transpose conjugate term (L*X*R)' to the (blcol,blrow) block
%   
%   else if <R=[]> (or ommited) 
%    adds the term (L*X*L') to the (blrow,blcol) diagonnal block (blrow=blcol)
%    
%
%   INPUTS:  
%     lmiindex = +n : left-hand side of the n-th LMI (left < right)
%              = -n : right-hand side of the n-th LMI (left < right)
%     blrow         : row index of the block 
%     blcol         : column index of the block   
%     Xindex   =  0 : constant term <L*R>
%              = +m : term on the m-th variable <X>
%              = -m : term on the transpose of the m-th variable <X.'>
%     L : matrix of appropriate dimensions, scalars accepted (default L=1 ).
%     R : matrix of appropriate dimensions, scalars accepted
%
%  
% quiz = sdmineq(quiz,lmiindex,Xindex,L,R);
%
%   Same use, but the terms are added regardless to the block decomposition.
%   This can be used for LMIs composed of a single block, but also for any other LMI. 
%   The terms (L*X*R)+(L*X*R)' or (L*X*L') are added to the whole LMI.
%   <L> and <R> must have respectively their number of rows and columns
%   equal to the LMI dimension (sum of the blocks dimensions)  
%
%
% quiz = sdmineq(quiz,[lmiindex blrow blcol],Xindex,L,R,K);
%
%   same use, but adds terms such that :
%    (L*(K kron X)*R) + (L*(K kron X)*R)'  or  (L*(K kron X)*L')  
%   where (K kron X) is the Kronecker product of K and X  
%
%
% quiz = sdmineq(quiz,[lmiindex blrow blcol],Xindex,L,R,K,-1);
%
%   same use, but adds terms such that :
%    (L*(X kron K)*R) + (L*(X kron K)*R)'  or  (L*(X kron K)*L')  
%
%  
% Note that complex valued L, R and K matrices are accepted. 
%  
% SEE ALSO sdmvar, sdmlmi, sdmeq and sdmsol.

%   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<3
    error('not enough input arguments');
  elseif nargin>7
    error('too many input arguments');
  elseif ~isa(quiz,'sdmpb')
    error('1st input argument must be a SDMPB object');
  end;      
  %%% check the input <lmidata>
  if isempty(lmidata)
    error('2nd input argument (LMI and block indexes) must contain at least the index of a LMI constraint')
  elseif min(size(lmidata))~=1
    error('2nd input argument (LMI and blocks indexes) must be a scalar or a vector');
  elseif any(round(lmidata)-lmidata)
    error('2nd input argument (LMI and blocks indexes) must be a vector of integers');
  end;
  lmiindex=lmidata(1);

  %%% check the input <lmiindex>
  if lmiindex<0
    side=-1;
  else
    side=1;
  end;  
  lmiindex=abs(lmiindex);
  if lmiindex>quiz.ineq.nb
    error('2nd input argument exceeds the number of inequality constraints');
  end;
  lmim=quiz.ineq.m(lmiindex);
  lmiM=quiz.ineq.M(lmiindex);
  nbrow=quiz.ineq.row(lmiindex);

  %%% check the inputs <blrow> and <blcol> if any
  if length(lmidata)<2
    % first case : no block partitionning
    lmidata(2)=0.5;
    lmidata(3)=0.5;
    % these two values are only added to remind that if there is no block 
    % then it means that all terms write "on the diagonal"
    mrow = 1;
    Mrow = nbrow;
    mcol = 1;
    Mcol = nbrow;
    % these four values are set to indicate that the block has the full
    % size of the LMI constraint
  else
    % second case : block partitionning
    if length(lmidata)>3
      warning('2nd input argument (LMI and blocks indexes) is too long! Only the first 3 will be used');
    elseif length(lmidata)==2
      error('2nd input argument (LMI and blocks indexes) must be either a scalar or a vector of 3 elements');
    end;
    if any(lmidata(2:3) <= 0)
      error('2nd and 3rd parameters of the 2nd input argument (line and column indexes of the block) must be positive');
    end;
    % Check if blocks indexes are defined
    if lmidata(2)>size(quiz.ineq.blockM{lmiindex})
      error('2nd parameter of the 2nd input argument (blocks'' row index) refers to a non-existing block.');
    elseif lmidata(3)>size(quiz.ineq.blockM{lmiindex})
      error('3rd parameter of the 2nd input argument (bolcks'' column index) refers to a non-existing block.');
    end;
    % Definition of the position of the block in the LMI rows and columns
    mrow = quiz.ineq.blockm{lmiindex};    mrow = mrow(lmidata(2));
    Mrow = quiz.ineq.blockM{lmiindex};    Mrow = Mrow(lmidata(2));
    mcol = quiz.ineq.blockm{lmiindex};    mcol = mcol(lmidata(3));
    Mcol = quiz.ineq.blockM{lmiindex};    Mcol = Mcol(lmidata(3));
  end;

  %%% check the input <Xindex>
  if ( ~isnumeric(Xindex) | (mod(abs(Xindex),2)~=1 & mod(abs(Xindex),2)~=0) | any(size(Xindex)-[1,1]))
    error('3rd input argument (index of the variable) does not fit the convensions');
  end;
  % check if the variable is transposed
  if Xindex<0, trans=1; else, trans=0; end;  
  Xindex=abs(Xindex);
  % check if the variable exists
  if Xindex>quiz.var.nb
    error('3rd input argument exceeds the number of variables');
  elseif Xindex>0
    varm=quiz.var.m(Xindex);
    varM=quiz.var.M(Xindex);
    varrow=quiz.var.row(Xindex);
    varcol=quiz.var.col(Xindex);    
  end;
  
  %%% check and append the input <L>
  if nargin<4, inL=1;  end;
  if isempty(inL), inL=1;  end;
  L=inL; 
  if ~isnumeric(L)
    error('4th input argument (left multiplier) must be a matrix or a scalar');
  % if scalar convert into a matrix with row block dimension
  elseif max(size(L))==1
    L=L*speye(Mrow-mrow+1);
  % check the dimensions fit with the specified block
  elseif ~isequal(size(L,1),(Mrow-mrow+1))
    error('4th input argument (left multiplier) must have the same number of rows as the specified block');
  end;
  % convert to a left multiplier regarless of blocks
  [ii,jj,s]=find(L);
  ii=ii+mrow-1;
  L=sparse(ii,jj,s,nbrow,size(L,2));
  
  %%% check and append the input <R>
  if nargin<5, inR=[]; end;
  % Default, if R is not defined
  if isempty(inR)  
    % first case : the block is on the diagonal
    if (lmidata(2) == lmidata(3)) 
      %%% check that the variable is symmetric
      if Xindex
        decvar=get(quiz,'vardec',Xindex);
        if size(decvar,1)~=size(decvar,2)
	  error(['the variable referenced by the 3rd input argument is not square : inconsistent with this syntax']);	
        end;
        if norm(decvar-decvar',1)
	  warning(['the variable referenced to by the 3rd argument is not Hermitian : inconsistent with this syntax']);	
        end;
      end;        
      inR=0.5*inL';  
    % second case : the block is off-diagonal
    else
      error('the 5th input argument (right multiplier) must be specified for off-diagonal blocks')
    end;
  end;
  R=inR;
  if ~isnumeric(R)
    error('5th input argument (right multiplier) must be a matrix or a scalar');
  % convert into a matrix with column block dimension
  elseif max(size(R))==1
    R=R*speye(Mcol-mcol+1); 
  % check the dimensions fit with the specified block
  elseif ~isequal(size(R,2),(Mcol-mcol+1))
    error('5th input argument (right multiplier) must have the same number of columns as the dimension of the specified block');
  end;
  % convert to a right multiplier regardless of blocks
  [ii,jj,s]=find(R);
  jj=jj+mcol-1;
  R=sparse(ii,jj,s,size(R,1),nbrow);
  % SeDuMi transforms the terms LXR into 0.5(LXR+R^*X^*L^*) 
  % to remove the 0.5 factor the R matrix is doubled
  R=2*R;
  
  %%% check the inputs <inK> and <kside> (only for variable terms)
  if Xindex
    if nargin<6, inK=[]; end;
    if isempty(inK), inK=1; end;
    K=inK;
    if ~isnumeric(K)
      error('6th input argument (matrix of the Kronecker product) must be a matrix');
    end;
    if nargin<7, kside=1; end;
    if abs(kside)~=1
      error('7th input argument (to invert Kronecker product) must be 1 or -1');
    end;
    % in case of a scalar variable convert the terms into
    % x*k*I=((k*I) kron x) to make it convenient with respect to the matrix dimensions
    if varrow==1 & varcol==1 & max(size(K))==1
      K=K*speye(size(L,2));
    end;
    [krow,kcol]=size(K);      
    %%% if L(X kron K)R , permute in L(K kron X)R
    if kside<0      
      L=L*sdmupq(varrow,krow);
      R=sdmupq(kcol,varcol)*R;
    end;
  end;
  
  %%% check that all multiplications and Kronecker products fit in size
  % for variable terms
  if Xindex
    % first in the case the variable is transposed
    if trans
      if size(L,2)~=varcol*krow
	if max(size(inK))==1
	  error('Term cannot be added: inconsistent number of columns of the variable');
	else
	  error('Term cannot be added: inconsistent number of rows in the resulting Kronecker product');
	end;
      elseif size(R,1)~=varrow*kcol
	if max(size(inK))==1
	  error('Term cannot be added: inconsistent number of rows of the variable');
	else
	  error('Term cannot be added: inconsistent number of columns in the resulting Kronecker product');
	end;
      end;
    % second in the case the variable is not transposed
    else
      if size(L,2)~=varrow*krow
	if max(size(inK))==1
	  error('Term cannot be added: inconsistent number of rows of the variable');
	else
	  error('Term cannot be added: inconsistent number of rows in the resulting Kronecker product');
	end;
      elseif size(R,1)~=varcol*kcol
	if max(size(inK))==1
	  error('Term cannot be added: inconsistent number of columns of the variable');
	else
	  error('Term cannot be added: inconsistent number of columns in the resulting Kronecker product');
	end;
      end;
    end;
  % for constant terms  
  else
    if size(L,2)~=size(R,1)
      error(['Number of columns of 4th input argument (left multiplier)' ...
	     ' does not fit the number of rows of 5th input argument (right multiplier)']);
      end;
  end;
    
  %%% All the tests are accomplished 
  %%% the remaining is generic
        
  %%% clear the previous solution if necessary
  quiz = sdmclear(quiz);
  
  %%% what to add if constant term
  if Xindex==0
    sic=length(quiz.ineq.c);
    [ii,jj,s]=find(side*sdmvec(L*R));
    ii=ii+lmim-1;
    add=sparse(ii,jj,s,sic,1);
    quiz.ineq.c = quiz.ineq.c-add;
  %%% if variable term
  else
    %%% add the term that depends on the variable terms (not conjugated)
    struclocal=quiz.var.struc;
    if norm(struclocal(varm:varM,:),1)
      %%% define the term for each value of K(ii,jj)
      add=sparse(kcol*varcol*krow*varrow,varcol*varrow);
      for ii=1:krow
	ER=kron(K(ii,:),speye(varcol));
	EL=kron(sparse(ii,1,1,krow,1),speye(varrow));
	add=add+kron(ER.',EL);      
      end;
      add=side*kron(R.',L)*add;
      %%% if X.' , permute decision variables of X
      if trans
	add=add*sdmupq(varrow,varcol);
      end;
      %%% add the term to the lmi
      siAt=size(quiz.ineq.At);
      [ii,jj,s]=find(add);
      ii=ii+lmim-1;
      jj=jj+varm-1;
      add=sparse(ii,jj,s,siAt(1),siAt(2));
      quiz.ineq.At=quiz.ineq.At+add;
    end;
    %%% now we do the same for the conjugate transpose elements of the variable
    struclocal=quiz.var.strucconj;
    if norm(struclocal(varm:varM,:),1)
      %%% conjugate what needs to be
      K=K';
      [krow,kcol]=size(K);
      varcol=quiz.var.row(Xindex);
      %varrow=(-varm+1+varM)/varcol;    
      varrow=quiz.var.col(Xindex);
      %%% define the term for each value of K(ii,jj)
      add=sparse(kcol*varcol*krow*varrow,varcol*varrow);
      for ii=1:krow
	ER=kron(K(ii,:),speye(varcol));
	EL=kron(sparse(ii,1,1,krow,1),speye(varrow));
	add=add+kron(ER.',EL);      
      end;
      add=side*kron((L').',R')*add;
      %%% if X , permute decision variables of X'
      if ~trans
	add=add*sdmupq(varrow,varcol);
      end;
      %%% add the term to the lmi
      siAtc=size(quiz.ineq.Atconj);
      [ii,jj,s]=find(add);
      ii=ii+lmim-1;
      jj=jj+varm-1;
      add=sparse(ii,jj,s,siAtc(1),siAtc(2));
      quiz.ineq.Atconj=quiz.ineq.Atconj+add;
    end;
  end;