Mercurial > forge

changeset 10201:6f1e79f2d62e octave-forge

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
control-devel: fix arx argument checking for degrees of numerator and denominator
author paramaniac
date Wed, 09 May 2012 10:36:17 +0000
parents 6e3b26da118c
children 46f08d64ad0e
files extra/control-devel/devel/CDplayerARX.m extra/control-devel/devel/GlassFurnaceARX.m extra/control-devel/inst/arx.m
diffstat 3 files changed, 112 insertions(+), 37 deletions(-) [+]
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/extra/control-devel/devel/CDplayerARX.m	Wed May 09 10:36:17 2012 +0000
@@ -0,0 +1,78 @@
+%{
+Contributed by:
+	Favoreel
+	KULeuven
+	Departement Electrotechniek ESAT/SISTA
+Kardinaal Mercierlaan 94
+B-3001 Leuven
+Belgium
+	wouter.favoreel@esat.kuleuven.ac.be
+Description:
+	Data from the mechanical construction of a CD player arm.  
+	The inputs are the forces of the mechanical actuators
+  	while the outputs are related to the tracking accuracy of the arm.
+  	The data was measured in closed loop, and then through a two-step
+  	procedure converted to open loop equivalent data
+    	The inputs are highly colored.
+Sampling:
+Number:
+	2048
+Inputs:
+	u: forces of the mechanical actuators
+Outputs:
+	y: tracking accuracy of the arm
+References:
+  	We are grateful to R. de Callafon of the
+    	Mechanical Engineering Systems and Control group of Delft, who
+    	provided us with these data.
+	
+	- Van Den Hof P., Schrama R.J.P., An Indirect Method for Transfer 
+	  Function Estimation From Closed Loop Data. Automatica, Vol. 29, 
+	  no. 6, pp. 1523-1527, 1993.
+
+Properties: 
+Columns:
+	Column 1: input u1
+	Column 2: input u2
+	Column 1: output y1
+	Column 2: output y2
+Category:
+	mechanical systems
+
+%}
+
+clear all, close all, clc
+
+load CD_player_arm-1.dat
+U=CD_player_arm_1(:,1:2);
+Y=CD_player_arm_1(:,3:4);
+
+
+dat = iddata (Y, U)
+
+% [sys, x0] = ident (dat, 15, 8)     % s=15, n=8
+sys = arx (dat, 4, [4 4])
+
+%[y, t] = lsim (sys, U, [], x0);
+%[y, t] = lsim (sys(:,1:2), U);
+
+[A, B] = filtdata (sys);
+%[A, B] = tfdata (sys);
+
+
+y1 = filter (B{1,1}, A{1,1}, U(:,1)) + filter (B{1,2}, A{1,2}, U(:,2));
+y2 = filter (B{2,1}, A{2,1}, U(:,1)) + filter (B{2,2}, A{2,2}, U(:,2));
+y = [y1, y2];
+
+t = 0:length(U)-1;
+
+err = norm (Y - y, 1) / norm (Y, 1)
+
+figure (1)
+p = columns (Y);
+for k = 1 : p
+  subplot (2, 1, k)
+  plot (t, Y(:,k), 'b', t, y(:,k), 'r')
+endfor
+
+
--- a/extra/control-devel/devel/GlassFurnaceARX.m	Wed May 09 03:49:09 2012 +0000
+++ b/extra/control-devel/devel/GlassFurnaceARX.m	Wed May 09 10:36:17 2012 +0000
@@ -49,7 +49,7 @@
 dat = iddata (Y, U)
 
 %[sys, x0] = ident (dat, 10, 5)     % s=10, n=5
-sys = arx (dat, 5, [5 5 5])
+sys = arx (dat, 5, 5)
 
 %[y, t] = lsim (sys, U, [], x0);
 [y, t] = lsim (sys(:, 1:3), U);
--- a/extra/control-devel/inst/arx.m	Wed May 09 03:49:09 2012 +0000
+++ b/extra/control-devel/inst/arx.m	Wed May 09 10:36:17 2012 +0000
@@ -9,25 +9,36 @@
   endif
   
   if (! isa (dat, "iddata"))
-    error ("arx: ");
-  endif
-  
-##  if (! is_real_scalar (na, nb))
-  if (! is_real_vector (na, nb))
-    error ("arx: ");
-    ## Test for integers
-    ## numel (nb) == size (dat, 3)
+    error ("arx: first argument must be an iddata dataset");
   endif
 
-  ## TODO: handle MIMO and multi-experiment data
+  ## p: outputs,  m: inputs,  ex: experiments
   [~, p, m, ex] = size (dat);
 
+  ## extract data  
   Y = dat.y;
   U = dat.u;
-  Ts = dat.tsam{1};
+  tsam = dat.tsam;
+
+  ## multi-experiment data requires equal sampling times  
+  if (ex > 1 && ! isequal (tsam{:}))
+    error ("arx: require equally sampled experiments");
+  else
+    tsam = tsam{1};
+  endif
+
   
-  max_nb = max (nb);
-  n = max (na, max_nb);
+  if (is_real_scalar (na, nb))
+    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)
 
 
   num = cell (p, m+1);
@@ -38,16 +49,17 @@
     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-1, 1)]); % TODO: multiple na
-      PhiU = arrayfun (@(x) toeplitz (U{e}(1:end-1, x), [U{e}(1, x); zeros(nb(x)-1, 1)]), 1:m, "uniformoutput", false);
-      Phi{e} = (horzcat (-PhiY, PhiU{:}))(n:end, :);
+      PhiY = toeplitz (Y{e}(1:end-1, i), [Y{e}(1, i); zeros(na(i)-1, 1)]);
+      ## create MISO Phi for every experiment
+      PhiU = arrayfun (@(x) toeplitz (U{e}(1:end-1, x), [U{e}(1, x); zeros(nb(i,x)-1, 1)]), 1:m, "uniformoutput", false);
+      Phi{e} = (horzcat (-PhiY, PhiU{:}))(n(i):end, :);
     endfor
   
     Theta = __theta__ (Phi, Y, i, n);
       
-    A = [1; Theta(1:na)];                     # a0 = 1, a1 = Theta(1), an = Theta(n)
-    ThetaB = Theta(na+1:end);                 # b0 = 0 (leading zero required by filt)
-    B = mat2cell (ThetaB, nb);
+    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);
 
@@ -55,7 +67,7 @@
     den(i, :) = repmat ({A}, 1, m+1);
   endfor
 
-  sys = filt (num, den, Ts);
+  sys = filt (num, den, tsam);
 
 endfunction
 
@@ -81,7 +93,7 @@
     V = V(:, 1:r);
     S = S(1:r);
     U = U(:, 1:r);
-    theta = V * (S .\ (U' * y{1}(n+1:end, i)));     # U' is the conjugate transpose
+    theta = V * (S .\ (U' * y{1}(n(i)+1:end, i)));     # U' is the conjugate transpose
   else
     ## multi-experiment dataset
     ## TODO: find more sophisticated formula than
@@ -92,7 +104,7 @@
     C = plus (tmp{:});
     
     ## PhiTY = (Phi1' Y1 + Phi2' Y2 + ...)
-    tmp = cellfun (@(Phi, Y) Phi.' * Y(n+1:end, i), phi, y, "uniformoutput", false);
+    tmp = cellfun (@(Phi, Y) Phi.' * Y(n(i)+1:end, i), phi, y, "uniformoutput", false);
     PhiTY = plus (tmp{:});
     
     ## pseudoinverse  Theta = C \ Phi'Y
@@ -100,18 +112,3 @@
   endif
   
 endfunction
-
-%{
-function Phi = __phi__ (dat, na, nb, ex)
-
-  ## avoid warning: toeplitz: column wins anti-diagonal conflict
-  ## therefore set first row element equal to y(1)
-  PhiY = toeplitz (Y(1:end-1, :), [Y(1, :); zeros(na-1, 1)]);
-  
-  ## PhiU = toeplitz (U(1:end-1, :), [U(1, :); zeros(nb-1, 1)]);
-  PhiU = arrayfun (@(x) toeplitz (U(1:end-1, x), [U(1, x); zeros(nb(x)-1, 1)]), 1:m, "uniformoutput", false);
-  Phi = horzcat (-PhiY, PhiU{:});
-  Phi = Phi(n:end, :);
-
-endfunction
-%}
\ No newline at end of file