QPpassive.m

Go to the documentation of this file.

function [SERCsubnew,SERDnew]=...
          QPpassive(bw,SERA,SERC,SERD,SERE,s,weightflag,...
                    auxflag,weightfactor,SERAsub,SERCsub,s2,s3,TOL);
%
%      ======================================================
%      =   Routine: QPpassive.m                             = 
%      =   Version 1.2                                      =
%      =   Last revised: 11.04.2007                         = 
%      =   Bjorn Gustavsen                                  =
%      =   SINTEF Energy Research, N-7465 Trondheim, NORWAY =
%      =   bjorn.gustavsen@sintef.no                        =
%      ======================================================
%
%function [SERCsubnew,SERDnew]=...
%          QPpassive(bw,SERA,SERC,SERD,SERE,s,weightflag,...
%                    auxflag,weightfactor,SERAsub,SERCsub,s2,s3,TOL);
%  
%  PURPOSE : Modify elements in SERC and SERD of rational approximation 
%                  
%                   N   SERCijm
%           Yij(s)=SUM(---------  ) +SERDij +s*SEREij ("ij" denotes element i,j)
%                  m=1 (s-SERAm)        
%
%            to enforce passivity at given frequency samples. 
%
% =================================================
%     INPUT :
% =================================================
%
%  bw  =0 --> Will modify elements on the diagonal of Rm, D              |0  1  2  3| 
%      =1 --> will modify elements on diagonal plus 1 band on each side  |1  0  1  2| 
%      =2 --> will modify elements on diagonal plus 2 bands on each side |2  1  0  1|
%      etc..                                                             |3  2  1  0|
% -----------------------------------------------------------------------------------                       
%  LEAST SQUARES PART:
%
%  SERA(N,1)     : Poles    
%  SERC(Nc,Nc,N) : Residues  
%  SERD(Nc,Nc)   : Constant term  
%  SERE(Nc,Nc)   : Capacitive term     
%
%  s(1,Ns)       : Frequency samples (jw [rad/sec])
%  
%  weightflag=1 --> weight=1 for all elements in Least Squares problem, at all freq.
%            =2 --> weight=1/sqrt(abs(Yijw)); indvidual element weighting
%            =3 --> weight=1/abs(Yijw)      ; indvidual element weighting
%            =4 --> weight=1/sqrt(norm(Yw)) ; common weighting for all matrix elements
%            =5 --> weight=1/norm(Yw)       ; common weighting for all matrix elements
%
%  auxflag   =1 --> Routine inserts additional samples for LS problem outside the
%                   fitting range. One sample at each pole location (SERAsub).
%
%  weightfactor : Factor defining the total weighting of additional samples outside
%                 the fitting range, relative to the total weighting of samples in the
%                 fitting range. Is only used if auxflag=1. 
%                 Default value: 1e-3
% -----------------------------------------------------------------------------------  
%  CONSTRAINT PART:
%
%  SERAsub(Nsub,1)     : Poles to be included in QP
%  SERCsub(Nc,Nc,Nsub) : Associated residues
%
%  s2(1,Ns2) (jw [rad/sec])  : List of frequencies. Will include in constraint problem
%                              eigenvalues which violate the passivity criterion at 
%                              these frequencies
%
%  s3(1,Ns3) (jw [rad/sec])  : List of frequencies. Will include in constraint problem 
%                              all eigenvalues at these frequencies.
%                               
% =================================================
%     OUTPUT :
% =================================================
%  SERCnew(Nc,Nc,Nsub) : Updated residues  
%  SERDnew(Nc,Nc)      : Updated constant term  
%
% *********************************************************************************** 
%  NOTE: This program is in the public domain and may be used by anyone. If the     *
%        program code (or a modified version) is used in a scientific work or       *
%        in a commercial program, then reference should be made to the following:   *
%        B. Gustavsen and A. Semlyen, "Enforcing Passivity for Admittance Matrices  *
%        Approximated by Rational Functions, IEEE Trans. Power Systems,             *  
%        vol. 16, no. 1, pp. 97-104, Feb. 2001.                                     *  
% ***********************************************************************************
%
% This file is part of the QPpassive.zip toolbox
%
% ***********************************************************************************
%
% Bug fixes:
% 21.10.2004 - Several bugs related to usage of cindexsub corrected. 
% 11.03.07   - Input parameter weighfactor correcctly included in the code.
%            - Proper dimensioning of array bigA2 
% ***********************************************************************************
 
%  11.04.2007: Added option for using CPLEX over quadprog, see below.
%
 
solver='quadprog';
%solver='cplex   '; %Requires installment of CPLEX (http://tomopt.com/)
 
 
SERCsubnew=SERCsub; SERDnew=SERD; SEREnew=SERE;
 
 
N=length(SERA);
Nsub=length(SERAsub); 
 
bigA=sparse([]);
bigb=[];
bigB=[];
bigc=[];
Ns=length(s);
Ns2=length(s2);
Nc=length(SERD);
I=eye(Nc);
M2mat=[];
 
 
%=======================================================
% Finding out which poles are complex :
%=======================================================
% LS problem:
cindex=zeros(1,N);
for m=1:N 
  if imag(SERA(m))~=0  
    if m==1 
      cindex(m)=1;
    else
      if cindex(m-1)==0 | cindex(m-1)==2
        cindex(m)=1; cindex(m+1)=2; 
      else
        cindex(m)=2;
      end
    end 
  end
end
 
% Constraint problem:
cindexsub=zeros(1,Nsub);
for m=1:Nsub 
  if imag(SERAsub(m))~=0 
    if m==1 
      cindexsub(m)=1; 
    else
      if cindexsub(m-1)==0 | cindexsub(m-1)==2
        cindexsub(m)=1; cindexsub(m+1)=2; 
      else
        cindexsub(m)=2;
      end
    end 
  end
end
 
%========================================
% LOOP FOR LEAST SQUARES PROBLEM:
%========================================
nnn=Nc;
for row=1:min(bw,Nc)
  nnn=nnn+(Nc-row);  %nnn is the number of elements in Y to be put in solution vector
end
 
if weightflag==4 | weightflag==5
  weightarr=ones(1,Ns);
  for k=1:Ns
    sk=s(k);
    for row=1:Nc
      for col=row:min(row+bw,Nc)
        Y(row,col)=SERD(row,col)+sk*SERE(row,col);
        Y(row,col)=Y(row,col)+sum(squeeze(SERC(row,col,1:N))./(s(k)-SERA(1:N)));
      end
    end
    if weightflag==4 
      weightarr(k)=1/sqrt(norm(Y));
    elseif weightflag==5 
      weightarr(k)=1/norm(Y);
    end
  end %for k=1:Ns
end %if weightflag==4 | weightflag==5
 
bigA=zeros(Ns*nnn,nnn*(Nsub+1));
for k=1:Ns
  sk=s(k);
  % Calculating matrix Mmat (coefficients for LS-problem)
  tell=0;
  offs=0;
 
  for row=1:Nc
    ind2=0;
    for col=row:min(row+bw,Nc)
      y=SERD(row,col)+sk*SERE(row,col);
      y=y+sum(squeeze(SERC(row,col,1:N))./(s(k)-SERA(1:N)));
      if weightflag==1
        weight=1;
      elseif weightflag==2
        weight=1/sqrt(abs(y));
      elseif weightflag==3
        weight=1/abs(y); 
      elseif weightflag==4 | weightflag==5
        weight=weightarr(k);
      else
        'ERROR!!!!!!!!!'
        weight=1;
      end
 
      if row==col
        faktor1=weight*1;
      else
        faktor1=weight*sqrt(2);
      end
 
      tell=tell+1;
      for m=1:Nsub
       if cindexsub(m)==0      %real pole
         Mmat(tell,offs+m)=faktor1*1/(sk-SERAsub(m));
       elseif cindexsub(m)==1  %complex pole, 1st part
         Mmat(tell,offs+m)=  faktor1* (1/(sk-SERAsub(m)) + 1/(sk-SERAsub(m)') );
         Mmat(tell,offs+m+1)=faktor1* (i/(sk-SERAsub(m)) - i/(sk-SERAsub(m)') ); 
        end
      end %m
      Mmat(tell,offs+Nsub+1)=faktor1*1; 
      offs=offs + (Nsub+1);
    end %col
  end %row
       
  bigA((k-1)*nnn+1:k*nnn,:)=Mmat;
  
end %for k=1:Ns
bigb=zeros(length(bigA(:,1)),1);
 
 
%========================================
% INTRODUCING SAMPLES OUTSIDE LS REGION: ONE SAMPLE PER POLE (s4)
%========================================
if auxflag==1
s4=[];    
tell=0;
for m=1:length(SERAsub)
  if cindexsub(m)==0 %real pole
    if ( ( abs(SERAsub(m))>s(Ns)/j ) | ( abs(SERAsub(m))<s(1)/j ) )
      tell=tell+1;
      s4(tell)=j*abs(SERAsub(m));
    end
  elseif cindexsub(m)==1 %complex pole, first part
    if ( ( abs(imag(SERAsub(m)))>s(Ns)/j ) | ( abs(imag(SERAsub(m)))<s(1)/j ) )
      tell=tell+1;
      s4(tell)=j*abs(imag(SERAsub(m)));
    end
  end
end
Ns4=length(s4);
 
 
%Weighting:
if weightflag==4 | weightflag==5
  weightarr4=ones(1,Ns4);
  for k=1:Ns4
    sk=s4(k);
    for row=1:Nc
      for col=row:min(row+bw,Nc)
        Y(row,col)=SERD(row,col)+sk*SERE(row,col);
        Y(row,col)=Y(row,col)+sum(squeeze(SERC(row,col,1:N))./(s(k)-SERA(1:N)));
      end
    end
    if weightflag==4 
      weightarr4(k)=1/sqrt(norm(Y));
    elseif weightflag==5 
      weightarr4(k)=1/norm(Y);
    end
  end %for k=1:Ns4
end %if weightflag==4 | weightflag==5
 
 
 
nnn=Nc;
for row=1:min(bw,Nc)
  nnn=nnn+(Nc-row);
end
 
bigA2=zeros(Ns4*nnn,nnn*(Nsub+1));
for k=1:Ns4
  sk=s4(k);
 
  % Calculating matrix Mmat (coefficients for LS-problem)
  tell=0;
  offs=0;
 
  for row=1:Nc
    ind2=0;
    for col=row:min(row+bw,Nc)
 
      if weightflag==1
        weight=1;
      elseif weightflag==2
        weight=1/sqrt(abs(y));
      elseif weightflag==3
        weight=1/abs(y); 
      elseif weightflag==4 | weightflag==5
        weight=weightarr4(k);
      else
        'ERROR!!!!!!!!!'
        weight=1;
      end
 
      if row==col
        faktor1=weight;
      else
        faktor1=weight*sqrt(2);
      end
      faktor1=weightfactor*faktor1;
      tell=tell+1;
      for m=1:Nsub
       if cindexsub(m)==0      %real pole
         Mmat(tell,offs+m)=faktor1*1/(sk-SERAsub(m));
       elseif cindexsub(m)==1  %complex pole, 1st part
         Mmat(tell,offs+m)=  faktor1* (1/(sk-SERAsub(m)) + 1/(sk-SERAsub(m)') );
         Mmat(tell,offs+m+1)=faktor1* (i/(sk-SERAsub(m)) - i/(sk-SERAsub(m)') ); %!!!!!
        end
      end %m
      Mmat(tell,offs+Nsub+1)=faktor1*1; 
      offs=offs + (Nsub+1);
 
    end %col
  end %row
       
  bigA2((k-1)*nnn+1:k*nnn,:)=Mmat;
  
end %for k=1:Ns
bigA=[bigA;bigA2];
bigb=zeros(length(bigA(:,1)),1);
 
end %if auxflag==1
 
%========================================
% LOOP FOR CONSTRAINT PROBLEM, TYPE #1 (violating eigenvalues in s2):
%========================================
for k=1:Ns2
  sk=s2(k);
  for row=1:Nc
    for col=1:Nc
      Y(row,col)=SERD(row,col)+sk*SERE(row,col);
      Y(row,col)=Y(row,col)+sum(squeeze(SERC(row,col,1:N))./(sk-SERA(1:N)));
    end
  end 
 
  %Calculating eigenvalues and eigenvectors: 
  [V,Z] = eig(real(Y)); Z=diag(Z);
  EE(:,k)=real(Z);
  if min(real(Z))<0  %any violations
 
    % Calculating matrix M2mat; matrix of partial derivatives :
    tell=0;
    offs=0;
 
    for row=1:Nc
      ind2=0;
      for col=row:min(row+bw,Nc)
        if row==col
          faktor2=1;
        else
          faktor2=2;
        end
        tell=tell+1;
        for m=1:Nsub
          if cindexsub(m)==0      %real pole
            Mmat2(tell,offs+m)=faktor2*1/(sk-SERAsub(m));
          elseif cindexsub(m)==1  %complex pole, 1st part
            Mmat2(tell,offs+m)=  faktor2* (1/(sk-SERAsub(m)) + 1/(sk-SERAsub(m)') );
            Mmat2(tell,offs+m+1)=faktor2* (i/(sk-SERAsub(m)) - i/(sk-SERAsub(m)') );
          end                  
        end %m
        Ndum=Nsub;
        Mmat2(tell,offs+Ndum+1)=faktor2*1; 
        offs=offs + (Nsub+1);
      end %col
    end %row
 
    for n=1:Nc
      tell=0;
      V1=V(:,n); 
      for row=1:Nc
        for col=row:min(row+bw,Nc) 
          if row==col
            qij=V1(row)^2;
          else
            qij=V1(row)*V1(col); 
          end
          tell=tell+1;
          Q(n,tell)=qij;
        end
      end
    end
 
 
    B=Q*Mmat2; 
    delz=real(Z);
    for n=1:Nc   %ustabilitet?
      if delz(n)<0
        bigB=[bigB;B(n,:)];
        bigc=[bigc;-TOL+delz(n)];
      end
    end
          
  end %if max(real(Z))>0  %any violations
  
end %k
 
 
 
%========================================
% LOOP FOR CONSTRAINT PROBLEM, TYPE #2: (all eigenvalues in s3):
%========================================
Ns3=length(s3);
for k=1:1:Ns3
  sk=s3(k);
  for row=1:Nc
    for col=1:Nc
      Y(row,col)=SERD(row,col)+sk*SERE(row,col);
      Y(row,col)=Y(row,col)+sum(squeeze(SERC(row,col,1:N))./(sk-SERA(1:N)));
    end
  end
  %Calculating eigenvalues and eigenvectors: 
  [V,Z] = eig(real(Y)); Z=diag(Z);
  EE(:,k)=real(Z);
 
    % Calculating matrix M2mat; matrix of partial derivatives :
    tell=0;
    offs=0;
 
    for row=1:Nc
      ind2=0;
      for col=row:min(row+bw,Nc)
        if row==col
          faktor2=1;
        else
          faktor2=2;
        end
        tell=tell+1;
        for m=1:Nsub
          if cindexsub(m)==0      %real pole
            Mmat2(tell,offs+m)=faktor2*1/(sk-SERA(m));
          elseif cindexsub(m)==1  %complex pole, 1st part
            Mmat2(tell,offs+m)=  faktor2* (1/(sk-SERAsub(m)) + 1/(sk-SERAsub(m)') );
            Mmat2(tell,offs+m+1)=faktor2* (i/(sk-SERAsub(m)) - i/(sk-SERAsub(m)') );
          end                  
        end %m
        Ndum=Nsub;
        Mmat2(tell,offs+Ndum+1)=faktor2*1; 
        offs=offs + (Nsub+1);
      end %col
    end %row
 
    for n=1:Nc
      tell=0;
      V1=V(:,n); 
      for row=1:Nc
        for col=row:min(row+bw,Nc) 
          if row==col
            qij=V1(row)^2;
          else
            qij=V1(row)*V1(col); 
          end
          tell=tell+1;
          Q(n,tell)=qij;
        end
      end
    end
 
 
    B=Q*Mmat2; 
    delz=real(Z);
    for n=1:Nc   
      bigB=[bigB;B(n,:)];
      bigc=[bigc;-TOL+delz(n)];
    end
  
end %for k=1:1:Ns3
 
if length(bigB)==0 
  return   %No passivity violations
end
 
 
c=bigc;
 
bigA=[real(bigA); imag(bigA)];
bigb=[real(bigb); imag(bigb)];
bigB=[real(bigB)];
bigA=sparse(bigA);
 
 
Acol=length(bigA(1,:));
for col=1:Acol
  Escale(col)=norm(bigA(:,col),2);
  bigA(:,col)=bigA(:,col)./Escale(col);
  if length(bigB)>0
    bigB(:,col)=bigB(:,col)./Escale(col);
  end
end
 
H=bigA.'*bigA;
ff=bigA.'*bigb;
H=full(H);
 
 
clear bigA;
 
warning off
tic
 
if solver=='quadprog'
     disp('Starting quadprog...')
    [dx,lambda]=quadprog(H,ff,-bigB,bigc,[],[],[],[],[],optimoptions('quadprog','MaxIterations',2e4));
elseif solver=='cplex   '
  c=0*ff;
  H0=c.';
 
  %F=sparse(H);
  A=-bigB; 
  H=sparse(H); %A=sparse(A);
  x_0=zeros(length(H(:,1)));
  b_U=bigc;
  b_L = -inf*ones(length(b_U),1);
  x_L = -inf*ones(length(x_0),1);
  x_U = inf*ones(length(x_0),1);
 
  disp('Starting CPLEX...') 
  Prob = qpAssign(H,c,A,b_L,b_U,x_L,x_U,x_0,'dust');
  Prob.MIP.cpxControl.QPMETHOD = 4;
  Prob.MIP.cpxControl.BARALG = 2; %2
  PriLev = 0;
  Prob.PriLevOpt = 0;
  [x, slack, v, rc, f_k] = cplex(c, A, x_L, x_U, b_L, b_U, ...
           Prob.MIP.cpxControl, [], PriLev, Prob, [], [], [], [], ...
           [], [], H); 
  dx=x;
 
else
  'ERROR: invalid value for solver'   
end
 
disp(['Elapsed time: ' num2str(toc) ' sec'])
warning on
 
dx=dx./Escale.';
 
 
 
%Converting dx so that the residues become complex:
tell=0;
for row=1:Nc
  for col=row:min(row+bw,Nc)
    for m=1:Nsub
      tell=tell+1;
      if cindexsub(m)==1
        r1=dx(tell); r2=dx(tell+1);
        dx(tell)=r1+i*r2; dx(tell+1)=r1-i*r2; 
      end
    end %m
    tell=tell+1; %Skip SERD
  end %col
end %row
 
 
%Updating SERC, SERD:
tell=0;
for row=1:Nc
  for col=row:min(row+bw,Nc)
    Ndum=Nsub;
    for m=1:Nsub 
      tell=tell+1;
      SERCsubnew(row,col,m)=SERCsubnew(row,col,m)+dx(tell);
      if row~=col
        SERCsubnew(col,row,m)=SERCsubnew(col,row,m)+dx(tell); 
      end
    end
    SERDnew(row,col)=SERDnew(row,col)+dx(tell+1);
    if row~=col
      SERDnew(col,row)=SERDnew(col,row)+dx(tell+1);
    end
    tell=tell+1;
  end
end