home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21fb.zip / octave / SCRIPTS.ZIP / scripts.fat / control / syscnct.m < prev    next >
Text File  |  1999-12-24  |  8KB  |  275 lines
  1. ## Copyright (C) 1996,1998 Auburn University.  All Rights Reserved.
  2. ##
  3. ## This file is part of Octave. 
  4. ##
  5. ## Octave is free software; you can redistribute it and/or modify it 
  6. ## under the terms of the GNU General Public License as published by the 
  7. ## Free Software Foundation; either version 2, or (at your option) any 
  8. ## later version. 
  9. ## 
  10. ## Octave is distributed in the hope that it will be useful, but WITHOUT 
  11. ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 
  12. ## FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License 
  13. ## for more details.
  14. ## 
  15. ## You should have received a copy of the GNU General Public License 
  16. ## along with Octave; see the file COPYING.  If not, write to the Free 
  17. ## Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111 USA. 
  18.  
  19. ## -*- texinfo -*- 
  20. ## @deftypefn {Function File } {@var{retsys} =} syscnct (@var{sys}, @var{out_idx},@var{in_idx}@{,@var{order}, @var{tol}@})
  21. ## Close the loop from specified outputs to respective specified inputs
  22. ## 
  23. ## @strong{Inputs}
  24. ## @table @var
  25. ## @item   sys
  26. ## system data structure
  27. ## @item   out_idx, in_idx
  28. ## list of connections indices; @math{y(out_idx(ii))}
  29. ## is connected to @math{u(in_idx(ii))}.
  30. ## @item   order
  31. ## logical flag (default = 0)
  32. ## @table @code
  33. ## @item    0
  34. ## leave inputs and outputs in their original order
  35. ## @item    1
  36. ## permute inputs and outputs to the order shown in the diagram below
  37. ## @end table
  38. ## @item     tol
  39. ## tolerance for singularities in algebraic loops default: 200@var{eps}
  40. ## @end table
  41. ## 
  42. ## @strong{Outputs}
  43. ##  @var{sys}: resulting closed loop system.
  44. ## 
  45. ## @strong{Method}
  46. ## @code{syscnct} internally permutes selected inputs, outputs as shown
  47. ##  below, closes the loop, and then permutes inputs and outputs back to their
  48. ##  original order
  49. ## @example
  50. ## @group
  51. ##                  ____________________
  52. ##  u_1       ----->|                  |----> y_1
  53. ##                  |        sys       |
  54. ##          old u_2 |                  |
  55. ## u_2* ---->(+)--->|                  |----->y_2 
  56. ## (in_idx)   ^     -------------------|    | (out_idx)
  57. ##            |                             |
  58. ##            -------------------------------
  59. ## @end group
  60. ## @end example
  61. ## The input that has the summing junction added to it has an * added to the end 
  62. ## of the input name.
  63. ## 
  64. ## @end deftypefn
  65.  
  66. function sys = syscnct (sys, output_list, input_list, order, tol)
  67.  
  68.   ## A. S. Hodel August 1995
  69.   ## modified by John Ingram July 1996
  70.  
  71.   if( (nargin < 3) | (nargin > 5) )
  72.     usage("retsys = syscnct(sys,output_list,input_list[,order,tol])");
  73.   endif
  74.  
  75.   ## check order
  76.   if(nargin <= 3)
  77.     order = 0;
  78.   elseif( (order != 0) & (order != 1) )
  79.     error("syscnct: order must be either 0 or 1")
  80.   endif
  81.  
  82.   if (nargin <= 4)
  83.     tol = 200*eps;
  84.   elseif( !is_sampl(tol) )
  85.     error("syscnct: tol must be a positive scalar");
  86.   elseif(tol > 1e2*sqrt(eps))
  87.     warning(["syscnct: tol set to large value=",num2str(tol), ...
  88.     ", eps=",num2str(eps)])
  89.   endif
  90.  
  91.   ## verify sizes,format of input, output lists
  92.   if( min(size(output_list))*min(size(input_list)) != 1)
  93.     error("output_list and input_list must be vectors");
  94.   else
  95.     lo = length(output_list);
  96.     li = length(input_list);
  97.     if(lo != li)
  98.       error("output_list and input_list must be of the same length")
  99.     endif
  100.     
  101.     if(is_dupl(output_list) | is_dupl(input_list) )
  102.       error("duplicate entry in input_list and/or output_list");
  103.     endif
  104.   endif
  105.   
  106.   [nc,nz,mm,pp] = sysdimen(sys);
  107.   nn = nc+nz;
  108.  
  109.   if( !is_struct(sys))
  110.     error("sys must be in structured system form")
  111.   elseif(pp < li)
  112.     error(["length(output_list)=",num2str(li),", sys has only ", ...
  113.     num2str(pp),"system outputs"])
  114.   elseif(mm < li)
  115.     error(["length(input_list)=",num2str(li),", sys has only ", ...
  116.     num2str(mm),"system inputs"])
  117.   endif
  118.  
  119.   ## check that there are enough inputs/outputs in the system for the lists
  120.   if(max(input_list) > mm) 
  121.     error("max(input_list) exceeds the number of inputs");
  122.   elseif(max(output_list) > pp)
  123.     error("max(output_list) exceeds the number of outputs");
  124.   endif
  125.  
  126.   output_list = reshape(output_list,1,length(output_list));
  127.  
  128.   ## make sure we're in state space form
  129.   sys = sysupdat(sys,'ss');
  130.  
  131.   ## permute rows and columns of B,C,D matrices into pseudo-dgkf form...
  132.   all_inputs = sysreord(mm,input_list);
  133.   all_outputs = sysreord(pp,output_list);
  134.  
  135.   [aa,bb,cc,dd] = sys2ss(sys);
  136.   bb = bb(:,all_inputs);
  137.   cc = cc(all_outputs,:);
  138.   dd = dd(all_outputs,all_inputs);
  139.  
  140.   yd = sysgetsg(sys,"yd");
  141.   yd = yd(all_outputs);
  142.  
  143.   ## m1, p1 = number of inputs, outputs that are not being connected
  144.   m1 = mm-li;
  145.   p1 = pp-li;
  146.  
  147.   ## m2, p2: 1st column, row of B, C that is being connected
  148.   m2 = m1+1;
  149.   p2 = p1+1;
  150.  
  151.   ## partition system into a DGKF-like form; the loop is closed around
  152.   ## B2, C2
  153.   if(m1 > 0)
  154.     B1 = bb(:,1:m1);
  155.     D21= dd(p2:pp,1:m1);
  156.   endif
  157.   B2 = bb(:,m2:mm);
  158.   if(p1 > 0)
  159.     C1 = cc(1:p1,:);
  160.     D12= dd(1:p1,m2:mm);
  161.   endif
  162.   C2 = cc(p2:pp,:);
  163.   if(m1*p1 > 0)
  164.     D11= dd(1:p1,1:m1);
  165.   endif
  166.   D22= dd(p2:pp,m2:mm);
  167.  
  168.   if(norm(D22))
  169.     warning("syscnct: possible algebraic loop, D22 non-zero");
  170.     D22i = (eye(size(D22))-D22);
  171.     C2h = D22i\C2;
  172.     if(m1 > 0)
  173.       D21h = D22i\D21;
  174.     endif
  175.     D22h = D22i\D22;
  176.   else
  177.     C2h = C2;
  178.     if(m1 > 0)
  179.       D21h = D21;
  180.     endif
  181.     D22h = D22;
  182.  
  183.   endif
  184.  
  185.   ## check cont state -> disc output -> cont state
  186.   dyi = find(yd(p2:pp));
  187.  
  188.   ## disp("syscnct: dyi=")
  189.   ## dyi
  190.   ## nc
  191.   ## disp("/syscnct");
  192.  
  193.   if( (nc > 0) & find(dyi > 0) )
  194.     B2con = B2(1:nc,dyi);    # connection to cont states
  195.     C2hd = C2h(dyi,1:nc);    # cont states -> outputs
  196.   else
  197.     B2con = C2hd = [];
  198.   endif
  199.  
  200.   if(max(size(B2con)) & max(size(C2hd)) )
  201.     if(norm(B2con*C2hd))
  202.       warning("syscnct: cont-state -> disc output -> cont state derivative");
  203.       warning("    connection made; resulting system may not be meaningful");
  204.     endif
  205.   endif
  206.  
  207.   Ac = aa+B2*C2h;
  208.   if(m1 > 0)
  209.     B1c = B1 + B2*D21h;
  210.   endif
  211.   B2c = B2*(eye(size(D22h)) + D22h);
  212.   if(p1*m1 > 0)
  213.     D11c = D11 + D12*D21h;
  214.   endif
  215.   if(p1 > 0)
  216.     C1c  = C1+D12*C2h;
  217.     D12c = D12*(eye(size(D22h))+D22h);
  218.   endif
  219.  
  220.   ## construct system data structure
  221.   if(m1 > 0)
  222.    Bc = [B1c, B2c];
  223.   else
  224.    Bc = B2c;
  225.   endif
  226.  
  227.   if(p1 > 0)
  228.     Cc = [C1c;C2h];
  229.   else
  230.     Cc = C2h;
  231.   endif
  232.  
  233.   if(m1*p1 > 0)
  234.     Dc = [D11c,D12c; D21h,D22h];
  235.   elseif(m1 > 0)
  236.     Dc = [D21h, D22h];
  237.   elseif(p1 > 0)
  238.     Dc = [D12c; D22h];
  239.   else
  240.     Dc = D22h;
  241.   endif 
  242.  
  243.   ## permute rows and columns of Bc, Cc, Dc back into original order
  244.   Im = eye(mm,mm);
  245.   Pi = Im(:,all_inputs);
  246.   back_inputs = Pi*[1:mm]';
  247.  
  248.   Ip = eye(pp,pp);
  249.   Po = Ip(:,all_outputs);
  250.   back_outputs = Po*[1:pp]';
  251.  
  252.   Bc = Bc(:,back_inputs);
  253.   Cc = Cc(back_outputs,:);
  254.   Dc = Dc(back_outputs,back_inputs);
  255.   yd = yd(back_outputs);
  256.  
  257.   ## rebuild system
  258.   Ts = sysgetts(sys);
  259.   [stnam,innam,outnam] = sysgetsg(sys);
  260.   sys = ss2sys(Ac,Bc,Cc,Dc,Ts,nc,nz,stnam,innam,outnam,find(yd));
  261.  
  262.   ## update connected input names
  263.   for ii = 1:length(input_list)
  264.     idx = input_list(ii);
  265.     strval = sprintf("%s*",nth(sysgetsg(sys,"in",idx),1) );
  266.     sys = syssetsg(sys,"in",strval,idx);
  267.   endfor
  268.   
  269.   ## maintain original system type if it was SISO
  270.   if    (strcmp(sysgetty(sys),"tf") )       sysupdat(sys,'tf');
  271.   elseif(strcmp(sysgetty(sys),"zp") )       sysupdat(sys,'zp');
  272.   endif
  273.  
  274. endfunction
  275.