octave_packages/control-2.3.52/@tf/__sys2ss__.m

   1 ## Copyright (C) 2009, 2011, 2012   Lukas F. Reichlin
   2 ##
   3 ## This file is part of LTI Syncope.
   4 ##
   5 ## LTI Syncope is free software: you can redistribute it and/or modify
   6 ## it under the terms of the GNU General Public License as published by
   7 ## the Free Software Foundation, either version 3 of the License, or
   8 ## (at your option) any later version.
   9 ##
  10 ## LTI Syncope is distributed in the hope that it will be useful,
  11 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13 ## GNU General Public License for more details.
  14 ##
  15 ## You should have received a copy of the GNU General Public License
  16 ## along with LTI Syncope.  If not, see <http://www.gnu.org/licenses/>.
  17
  18 ## -*- texinfo -*-
  19 ## TF to SS conversion.
  20 ## Reference:
  21 ## Varga, A.: Computation of irreducible generalized state-space realizations.
  22 ## Kybernetika, 26:89-106, 1990
  23
  24 ## Special thanks to Vasile Sima and Andras Varga for their advice.
  25 ## Author: Lukas Reichlin <lukas.reichlin@gmail.com>
  26 ## Created: October 2009
  27 ## Version: 0.4.1
  28
  29 function [retsys, retlti] = __sys2ss__ (sys)
  30
  31   ## TODO: determine appropriate tolerance from number of inputs
  32   ##       (since we multiply all denominators in a row), index, ...
  33   ##       default tolerance from TB01UD is TOLDEF = N*N*EPS
  34
  35   ## SECRET WISH: a routine which accepts individual denominators for
  36   ##              each channel and which supports descriptor systems
  37
  38   [p, m] = size (sys);
  39   [num, den] = tfdata (sys);
  40
  41   len_num = cellfun (@length, num);
  42   len_den = cellfun (@length, den);
  43
  44   ## check for properness
  45   ## tfpoly ensures that there are no leading zeros
  46   tmp = len_num > len_den;
  47   if (any (tmp(:)))       # non-proper transfer function
  48     ## separation into strictly proper and polynomial part
  49     [numq, numr] = cellfun (@deconv, num, den, "uniformoutput", false);
  50     numq = cellfun (@__remove_leading_zeros__, numq, "uniformoutput", false);
  51     numr = cellfun (@__remove_leading_zeros__, numr, "uniformoutput", false);
  52
  53     ## minimal state-space realization for the proper part
  54     [a1, b1, c1] = __proper_tf2ss__ (numr, den, p, m);
  55     e1 = eye (size (a1));
  56
  57     ## minimal realization for the polynomial part
  58     [e2, a2, b2, c2] = __polynomial_tf2ss__ (numq, p, m);
  59
  60     ## assemble irreducible descriptor realization
  61     e = blkdiag (e1, e2);
  62     a = blkdiag (a1, a2);
  63     b = vertcat (b1, b2);
  64     c = horzcat (c1, c2);
  65     retsys = dss (a, b, c, [], e);
  66   else                    # proper transfer function
  67     [a, b, c, d] = __proper_tf2ss__ (num, den, p, m);
  68     retsys = ss (a, b, c, d);
  69   endif
  70
  71   retlti = sys.lti;       # preserve lti properties such as tsam
  72
  73 endfunction
  74
  75
  76 ## transfer function to state-space conversion for proper models
  77 function [a, b, c, d] = __proper_tf2ss__ (num, den, p, m)
  78
  79   ## new cells for the TF of same row denominators
  80   numc = cell (p, m);
  81   denc = cell (p, 1);
  82
  83   ## multiply all denominators in a row and
  84   ## update each numerator accordingly
  85   ## except for single-input models and those
  86   ## with equal denominators in a row
  87   for i = 1 : p
  88     if (m == 1 || isequal (den{i,:}))
  89       denc(i) = den{i,1};
  90       numc(i,:) = num(i,:);
  91     else
  92       denc(i) = __conv__ (den{i,:});
  93       for j = 1 : m
  94         idx = setdiff (1:m, j);
  95         numc(i,j) = __conv__ (num{i,j}, den{i,idx});
  96       endfor
  97     endif
  98   endfor
  99
 100   len_numc = cellfun (@length, numc);
 101   len_denc = cellfun (@length, denc);
 102
 103   ## check for properness
 104   ## tfpoly ensures that there are no leading zeros
 105   ## tmp = len_numc > repmat (len_denc, 1, m);
 106   ## if (any (tmp(:)))
 107   ##   error ("tf: tf2ss: system must be proper");
 108   ## endif
 109
 110   ## create arrays and fill in the data
 111   ## in a way that Slicot TD04AD can use
 112   max_len_denc = max (len_denc(:));
 113   ucoeff = zeros (p, m, max_len_denc);
 114   dcoeff = zeros (p, max_len_denc);
 115   index = len_denc-1;
 116
 117   for i = 1 : p
 118     len = len_denc(i);
 119     dcoeff(i, 1:len) = denc{i};
 120     for j = 1 : m
 121       ucoeff(i, j, len-len_numc(i,j)+1 : len) = numc{i,j};
 122     endfor
 123   endfor
 124
 125   tol = min (sqrt (eps), eps*prod (index));
 126   [a, b, c, d] = sltd04ad (ucoeff, dcoeff, index, tol);
 127
 128 endfunction
 129
 130
 131 ## realization of the polynomial part according to Andras' paper
 132 function [e2, a2, b2, c2] = __polynomial_tf2ss__ (numq, p, m)
 133
 134   len_numq = cellfun (@length, numq);
 135   max_len_numq = max (len_numq(:));
 136   numq = cellfun (@(x) prepad (x, max_len_numq, 0, 2), numq, "uniformoutput", false);
 137   f = @(y) cellfun (@(x) x(y), numq);
 138   s = 1 : max_len_numq;
 139   D = arrayfun (f, s, "uniformoutput", false);
 140
 141   e2 = diag (ones (p*(max_len_numq-1), 1), -p);
 142   a2 = eye (p*max_len_numq);
 143   b2 = vertcat (D{:});
 144   c2 = horzcat (zeros (p, p*(max_len_numq-1)), -eye (p));
 145
 146   ## remove uncontrollable part
 147   [a2, e2, b2, c2] = sltg01jd (a2, e2, b2, c2, 0.0, true, 1, 2);
 148
 149 endfunction
 150
 151
 152 ## convolution for more than two arguments
 153 function vec = __conv__ (vec, varargin)
 154
 155   if (nargin == 1)
 156     return;
 157   else
 158     for k = 1 : nargin-1
 159       vec = conv (vec, varargin{k});
 160     endfor
 161   endif
 162
 163 endfunction
 164
 165
 166 ## remove leading zeros from polynomial vector
 167 function p = __remove_leading_zeros__ (p)
 168
 169   idx = find (p != 0);
 170
 171   if (isempty (idx))
 172     p = 0;
 173   else
 174     p = p(idx(1) : end);  # p(idx) would remove all zeros
 175   endif
 176
 177 endfunction