[Octave-cvsupdate] SF.net SVN: octave:[10466] trunk/octave-forge/extra/control-devel/inst/arx. m

[<par...@us...>]

Revision: 10466
          http://octave.svn.sourceforge.net/octave/?rev=10466&view=rev
Author:   paramaniac
Date:     2012-05-20 09:18:53 +0000 (Sun, 20 May 2012)
Log Message:
-----------
control-devel: ARX models have individual error inputs for each output

Modified Paths:
--------------
    trunk/octave-forge/extra/control-devel/inst/arx.m

Modified: trunk/octave-forge/extra/control-devel/inst/arx.m
===================================================================
--- trunk/octave-forge/extra/control-devel/inst/arx.m	2012-05-19 22:15:57 UTC (rev 10465)
+++ trunk/octave-forge/extra/control-devel/inst/arx.m	2012-05-20 09:18:53 UTC (rev 10466)
@@ -27,13 +27,6 @@
 function sys = arx (dat, na, nb)
 
   ## TODO: delays
-  
-  ## FIXME: MIMO models have p error inputs, not just 1
-  
-  ## FIXME: Unlike SISO or MISO models, the transfer function
-  ##        B(z) over A(z) is wrong for MIMO models.
-  ##        We need a Matrix Fraction Description (MFD)
-  ##        y = A^-1(q) B(q) u(t) + A^-1(q) e(t)
 
   if (nargin != 3)
     print_usage ();
@@ -60,24 +53,24 @@
 
   
   if (is_real_scalar (na, nb))
-    na = repmat (na, p, 1);     # na(p-by-1)
-    nb = repmat (nb, p, m);     # nb(p-by-m)
+    na = repmat (na, p, 1);                         # na(p-by-1)
+    nb = repmat (nb, p, m);                         # nb(p-by-m)
   elseif (! (is_real_vector (na) && is_real_matrix (nb) \
           && rows (na) == p && rows (nb) == p && columns (nb) == m))
     error ("arx: require na(%dx1) instead of (%dx%d) and nb(%dx%d) instead of (%dx%d)", \
             p, rows (na), columns (na), p, m, rows (nb), columns (nb));
   endif
 
-  max_nb = max (nb, [], 2);     # one maximum for each row/output, max_nb(p-by-1)
-  n = max (na, max_nb);         # n(p-by-1)
+  max_nb = max (nb, [], 2);                         # one maximum for each row/output, max_nb(p-by-1)
+  n = max (na, max_nb);                             # n(p-by-1)
 
 
-  num = cell (p, m+1);
-  den = cell (p, m+1);
+  num = cell (p, m+p);
+  den = cell (p, m+p);
   
-  for i = 1 : p     # for every output
+  for i = 1 : p                                     # for every output
     Phi = cell (ex, 1);
-    for e = 1 : ex      # for every experiment  
+    for e = 1 : ex                                  # for every experiment  
       ## avoid warning: toeplitz: column wins anti-diagonal conflict
       ## therefore set first row element equal to y(1)
       PhiY = toeplitz (Y{e}(1:end-1, i), [Y{e}(1, i); zeros(na(i)-1, 1)]);
@@ -88,16 +81,25 @@
   
     Theta = __theta__ (Phi, Y, i, n);
       
-    A = [1; Theta(1:na(i))];                     # a0 = 1, a1 = Theta(1), an = Theta(n)
-    ThetaB = Theta(na(i)+1:end);                 # b0 = 0 (leading zero required by filt)
+    A = [1; Theta(1:na(i))];                        # a0 = 1, a1 = Theta(1), an = Theta(n)
+    ThetaB = Theta(na(i)+1:end);                    # b0 = 0 (leading zero required by filt)
     B = mat2cell (ThetaB, nb(i,:));
     B = reshape (B, 1, []);
     B = cellfun (@(x) [0; x], B, "uniformoutput", false);
 
-    num(i, :) = [B, {1}];
-    den(i, :) = repmat ({A}, 1, m+1);
+    Be = repmat ({0}, 1, p);
+    Be(i) = 1;
+    num(i, :) = [B, Be];
+    den(i, :) = repmat ({A}, 1, m+p);
   endfor
 
+  ## A(q) y(t) = B(q) u(t) + e(t)
+  ## there is only one A per row
+  ## B(z) and A(z) are a Matrix Fraction Description (MFD)
+  ## y = A^-1(q) B(q) u(t) + A^-1(q) e(t)
+  ## since A(q) is a diagonal polynomial matrix, its inverse is trivial:
+  ## the corresponding transfer function has common row denominators.
+
   sys = filt (num, den, tsam);
 
 endfunction
@@ -147,8 +149,8 @@
   ## Th = Ph \ Y = Ph+ Y
   ## Th = V S+ U* Y,   S+ = 1 ./ diag (S)
 
-  [U, S, V] = svd (A, 0);                         # 0 for "economy size" decomposition
-  S = diag (S);                                   # extract main diagonal
+  [U, S, V] = svd (A, 0);                           # 0 for "economy size" decomposition
+  S = diag (S);                                     # extract main diagonal
   r = sum (S > eps*S(1));
   if (r < length (S))
     warning ("arx: rank-deficient coefficient matrix");
@@ -158,6 +160,6 @@
   V = V(:, 1:r);
   S = S(1:r);
   U = U(:, 1:r);
-  x = V * (S .\ (U' * b));                    # U' is the conjugate transpose
+  x = V * (S .\ (U' * b));                          # U' is the conjugate transpose
 
 endfunction

This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.