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