decouple.m

内模控制器（IMC）工具箱。包括参数整定、PID控制器参数转换等

function [qnummt,qdenmt,qdeadmt,Gnum,Gdead,input_delay,n_dec,TR,TL]=decouple(mnummt,mdenmt,mdeadmt,y)
%function [q,qnum,qden,qdead,Gnum,Gden,Gdead,input_delay]=decouple(mnummt,mdenmt,mdeadmt,y)
%Fine a posible decouple control system. It return an invertible part of the model (qxxx),
%non-invertible part (Gxxx), and the input delay.
global Priority

[m,n]=size(mnummt);
if iscell(mdeadmt) % Conver a cell array of dead time to a matrix
temp=zeros(m,n);
   for i=1:m
      for j=1:n
         temp(i,j)=mt2dt(mdeadmt{i,j},y);
      end
   end
mdeadmt=temp;
end
[input_delay,dec_M,qdeadmt]=dtarray(mdeadmt,y);
common_RHPZ=find_common(mnummt);
if all(iscellempty(common_RHPZ))
   [n_rhpz]=detect_RHPZ(mnummt,mdenmt,qdeadmt,0);
else
   for i=1:m
      for j=1:n
         if ~isempty(common_RHPZ{i})
            mnummt{i,j}=deconv(mnummt{i,j},common_RHPZ{i});
         end
      end
   end
   n_rhpz=detect_RHPZ(mnummt,mdenmt,qdeadmt,0);
end
if n_rhpz
   % seek all (n-1)x(n-1) submatrix
   ButtonName=questdlg(['The whole (m x n) system could not be fully decoupled.?',...
         ' The software can seek to decouple a square sub-system. Do you want to proceed?'],...
      'A right half plane tranmission zero is detected in the given process model.');
   switch ButtonName
   case 'Yes'
      OK=0;
      count=1;
      prompt={'Least important output --> Most important output (separate by space or comma)',...
            'Least preferred control effort --> Most preferred control effort'};
      while ~OK
         if count>1
            dlgTitle='The priorities of ALL outputs and efforts are required!';
         else
            dlgTitle='Please sort the priorities of all outputs and all control efforts';
         end
         Priority=inputdlg(prompt,dlgTitle,1,Priority);
         count=count+1;
         output_Priority=str2num(Priority{1});
         effort_Priority=str2num(Priority{2});
         if length(output_Priority)==m & length(effort_Priority)==n
            OK=1;
         end
      end
      for level=1:min([m n])-1
         %seeking for a submatrix that have no RHPZ.
         [n_rhpz,index]=seek(mnummt,mdenmt,mdeadmt,level,output_Priority,effort_Priority);
         if n_rhpz==0
            break
         end
      end
      if n_rhpz==0
         % calualte the triangular decoupling
         % Re-arrange the process model to an appropriate form
         TR=eye(m);TL=TR;
         [h,temp]=size(index);
         for i=1:h
         mnummt(m+1,:)=mnummt(index(i,1),:);mnummt(index(i,1),:)=[];
         mdenmt(m+1,:)=mdenmt(index(i,1),:);mdenmt(index(i,1),:)=[];
         mdeadmt(m+1,:)=mdeadmt(index(i,1),:);mdeadmt(index(i,1),:)=[];
         TL(m+1,:)=TL(index(i,1),:);TL(index(i,1),:)=[];
         mnummt(:,n+1)=mnummt(:,index(i,2));mnummt(:,index(i,2))=[];
         mdenmt(:,n+1)=mdenmt(:,index(i,2));mdenmt(:,index(i,2))=[];
         mdeadmt(:,n+1)=mdeadmt(:,index(i,2));mdeadmt(:,index(i,2))=[];
         TR(:,n+1)=TR(:,index(i,2));TR(:,index(i,2))=[];
         end
         [input_delay,dec_M,qdeadmt]=dtarray(mdeadmt(1:m-h,1:n-h),y);
         T=transform(qdeadmt);
         dec_M(m-h+1:m)=0;
         for j=n-h+1:n
            input_delay(j)=max(dec_M'-mdeadmt(:,j));
         end
         for i=1:m
            mdeadmt(i,:)= (mdeadmt(i,:)+input_delay)-dec_M(i);
         end
         input_delay=input_delay*inv(TR); % input delay applied to the original system
         mnummt(:,1:n-h)=sw_col(mnummt(:,1:n-h),T);
         mdenmt(:,1:n-h)=sw_col(mdenmt(:,1:n-h),T);
         mdeadmt(:,1:n-h)=sw_col(mdeadmt(:,1:n-h),T);
         TR(:,1:n-h)=sw_col(TR(:,1:n-h),T);
         % At this point we have a re-arranged process model with
         % a sub-system having no RHPZ on the upper left conner.
         % (mnummt, mdenmt and mdeadmt).
         % A decoupled matrix G+(s) having only the main diagonal
         % element containing common RHPZs and dead time in each row
         % of process model is factored out.
         % (common_RHPZ and dec_M)
         % n_rhpz=detect_RHPZ(mnummt,mdenmt,mdeadmt,1);
         % common_RHPZ{:}
         % dec_M
         %Gnum=eye(n,n);Gden=Gnum;Gdead=diag(dec_M);
         Gdead=dec_M; n_dec=m-h;
         for i=1:n
            for j=1:n
               if i==j
                  Gnum{i,j}=common_RHPZ{i};
               else
                  Gnum{i,j}=0;
              end
              if i==n | j==n
                 mnummt{i,j}=[];
                 mdenmt{i,j}=1;
              end
              if j==n & i==n
                 mnummt{i,j}=1;
              end
            end
         end
         qnummt=mnummt;qdenmt=mdenmt;qdeadmt=mdeadmt;
      else
         errordlg('The system is unlikely to be controlable.');
         qnummt={};qdenmt={};qdeadmt={};Gnum={};
         Gdead={};input_delay=[];TR=[];TL=[];n_dec=m;
      end
   otherwise
      qnummt={};qdenmt={};qdeadmt={};Gnum={};
      Gdead={};input_delay=[];TR=[];TL=[];
   end
else
   % fully decouple
   qnummt=mnummt;
   qdenmt=mdenmt;n_dec=m;
   Gnum={};Gdead=dec_M;
   TL=eye(n);TR=transform(qdeadmt);
end
if ~iscell(qdeadmt)
temp=cell(m,n);
   for i=1:m
      for j=1:n
         temp{i,j}=qdeadmt(i,j);
      end
   end
	qdeadmt=temp;
end  

function common_RHPZ=find_common(mnummt)
% Find common RHPZ in the same row of mnummt
[m,n]=size(mnummt);
for i=1:m
   temp=gcdpoly(mnummt{i,1},mnummt{i,2});
   if length(temp)>1
      for j=3:n
         temp=gcdpoly(temp,mnummt{i,j});
         if length(temp)==1
            break
         end
      end
   end
   if length(temp)>1
      temp=roots(temp);
      temp=temp(find(temp>0));
      temp=poly(temp);
      if length(temp)>1
         common_RHPZ{i}=temp/temp(end);
      else
         common_RHPZ{i}=[];
      end
   else
      common_RHPZ{i}=[];
   end
end

function [num,den,dead]=approximate(nummt,denmt,deadmt)
% Functuon to approximate the transfer function between
% un and yn, form the MIMO model arranged to have all
% RHPZ in the last output
[m,n]=size(nummt);
for i=1:m-1
   for j=1:n-1
      tmp(i,j)=length(denmt{i,j})-length(nummt{i,j});
      temp=['(' mtx2str(denmt{i,j},2) ')'];
      if isempty(mtx2str(deadmt{i,j},4))
         m{i,j}=[mtx2str(nummt{i,j},2) '/' temp];
      else
         m{i,j}=[mtx2str(nummt{i,j},2) '/' temp '*' mtx2str(deadmt{i,j},4)];
      end
   end
end
order=max(tmp);
model=cell2str(m);
s=0; K=inv(eval(model));
s=1e17; tau=diag(s.^order)*eval(model); tau=inv(tau);
tau=tau.*K;

function [n_rhpz,index]=seek(mnummt,mdenmt,mdeadmt,level,output_Priority,effort_Priority)
[m,n]=size(mnummt);
for i=1:m
   for j=1:n
      index=[output_Priority(i) effort_Priority(j)];
      tmpnum=mnummt;tmpden=mdenmt;tmpdead=mdeadmt;
      tmpnum(output_Priority(i),:)=[];
      tmpnum(:,effort_Priority(j))=[];
      tmpden(output_Priority(i),:)=[];tmpden(:,effort_Priority(j))=[];
      tmpdead(output_Priority(i),:)=[];tmpdead(:,effort_Priority(j))=[];
      if level > 1
         %seek a deeper level
         opt_Priority=output_Priority;
         temp=find(opt_Priority>opt_Priority(i));
         opt_Priority(temp)=opt_Priority(temp)-1;
         eff_Priority=effort_Priority;
         temp=find(eff_Priority>eff_Priority(j));
         eff_Priority(temp)=eff_Priority(temp)-1;
         opt_Priority(i)=[];eff_Priority(j)=[];
         [n_rhpz,temp]=seek(tmpnum,tmpden,tmpdead,level-1,opt_Priority,eff_Priority);
         index=[index;temp];
      else
         [input_delay,dec_M,qdeadmt]=dtarray(tmpdead,[]);
         n_rhpz=detect_RHPZ(tmpnum,tmpden,qdeadmt,0);
      end
      if n_rhpz==0
         break; % break out all for loops because the required invertible
      end		 % part of the model is found
   end
   if n_rhpz==0
      break; 
   end
end