Menu
▾
▴
octave-cvsupdate
|
From: <par...@us...> - 2012-04-30 15:57:28
|
Revision: 10335
http://octave.svn.sourceforge.net/octave/?rev=10335&view=rev
Author: paramaniac
Date: 2012-04-30 15:57:22 +0000 (Mon, 30 Apr 2012)
Log Message:
-----------
control-devel: add arx function
Added Paths:
-----------
trunk/octave-forge/extra/control-devel/inst/arx.m
Added: trunk/octave-forge/extra/control-devel/inst/arx.m
===================================================================
--- trunk/octave-forge/extra/control-devel/inst/arx.m (rev 0)
+++ trunk/octave-forge/extra/control-devel/inst/arx.m 2012-04-30 15:57:22 UTC (rev 10335)
@@ -0,0 +1,36 @@
+## Author: Lukas Reichlin <luk...@gm...>
+## Created: April 2012
+## Version: 0.1
+
+function sys = arx (dat, na, nb)
+
+ if (nargin != 3)
+ print_usage ();
+ endif
+
+ if (! isa (dat, "iddata"))
+ error ("arx: ");
+ endif
+
+ if (! is_real_scalar (na, nb))
+ error ("arx: ");
+ ## Test for integers
+ endif
+
+ ## TODO: handle MIMO and multi-experiment data
+
+ Y = dat.y{1};
+ U = dat.u{1};
+
+ PhiY = toeplitz ([0; Y(1:end-1, :)], zeros (1, na));
+ PhiU = toeplitz ([0; U(1:end-1, :)], zeros (1, nb));
+ Phi = [PhiY, PhiU];
+
+ Theta = Phi \ Y;
+
+ A = Theta(1:n);
+ B = Theta(n+1:end);
+
+ sys = tf ({B, 1}, {A, A}, dat.tsam);
+
+endfunction
\ No newline at end of file
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <par...@us...> - 2012-04-30 17:13:25
|
Revision: 10337
http://octave.svn.sourceforge.net/octave/?rev=10337&view=rev
Author: paramaniac
Date: 2012-04-30 17:13:19 +0000 (Mon, 30 Apr 2012)
Log Message:
-----------
control-devel: fix arx function partially, error term needs still some attention
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-04-30 16:51:51 UTC (rev 10336)
+++ trunk/octave-forge/extra/control-devel/inst/arx.m 2012-04-30 17:13:19 UTC (rev 10337)
@@ -23,14 +23,15 @@
U = dat.u{1};
Ts = dat.tsam{1};
- PhiY = toeplitz (Y(1:end-1, :), zeros (1, na));
- PhiU = toeplitz (U(1:end-1, :), zeros (1, nb));
- Phi = [-PhiY, PhiU]
+ ## avoid warning: toeplitz: column wins anti-diagonal conflict
+ PhiY = toeplitz (Y(1:end-1, :), [Y(1, :); zeros(na-1, 1)]);
+ PhiU = toeplitz (U(1:end-1, :), [U(1, :); zeros(nb-1, 1)]);
+ Phi = [-PhiY, PhiU];
Theta = Phi \ Y(2:end, :);
- A = [1; Theta(1:na)]; % ???
- B = Theta(na+1:end);
+ A = [1; Theta(1:na)]; # a0 = 1, a1 = Theta(1), an = Theta(n)
+ B = Theta(na+1:end); # b0 = 0 (leading zeros are removed by tf, no need to add one)
sys = tf ({B, 1}, {A, A}, Ts);
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <par...@us...> - 2012-05-06 15:53:41
|
Revision: 10367
http://octave.svn.sourceforge.net/octave/?rev=10367&view=rev
Author: paramaniac
Date: 2012-05-06 15:53:35 +0000 (Sun, 06 May 2012)
Log Message:
-----------
control-devel: use filt instead of tf
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-06 14:32:55 UTC (rev 10366)
+++ trunk/octave-forge/extra/control-devel/inst/arx.m 2012-05-06 15:53:35 UTC (rev 10367)
@@ -28,13 +28,13 @@
PhiU = toeplitz (U(1:end-1, :), [U(1, :); zeros(nb-1, 1)]);
Phi = [-PhiY, PhiU];
- Theta = Phi \ Y(2:end, :) # naive formula
+ ## Theta = Phi \ Y(2:end, :); # naive formula
## solve linear least squares problem by pseudoinverse
## the pseudoinverse is computed by singular value decomposition
## M = U S V* ---> M+ = V S+ U*
## Th = Ph \ Y = Ph+ Y
- ## Th = V S+ U' Y, S+ = 1 ./ diag (S)
+ ## Th = V S+ U* Y, S+ = 1 ./ diag (S)
[U, S, V] = svd (Phi, 0); # 0 for "economy size" decomposition, U overwrites input U
S = diag (S); # extract main diagonal
@@ -42,11 +42,11 @@
V = V(:, 1:r);
S = S(1:r);
U = U(:, 1:r);
- Theta = V * (S .\ (U' * Y(2:end, :))) # U' is the conjugate transpose
+ Theta = V * (S .\ (U' * Y(2:end, :))); # U' is the conjugate transpose
- A = [1; Theta(1:na)]; # a0 = 1, a1 = Theta(1), an = Theta(n)
- B = Theta(na+1:end); # b0 = 0 (leading zeros are removed by tf, no need to add one)
+ A = [1; Theta(1:na)]; # a0 = 1, a1 = Theta(1), an = Theta(n)
+ B = [0; Theta(na+1:end)]; # b0 = 0 (leading zero required by filt)
- sys = tf ({B, 1}, {A, A}, Ts);
+ sys = filt ({B, 1}, {A, A}, Ts);
endfunction
\ No newline at end of file
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <par...@us...> - 2012-05-07 12:40:00
|
Revision: 10373
http://octave.svn.sourceforge.net/octave/?rev=10373&view=rev
Author: paramaniac
Date: 2012-05-07 12:39:54 +0000 (Mon, 07 May 2012)
Log Message:
-----------
control-devel: fix arx initial conditions
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-07 04:39:42 UTC (rev 10372)
+++ trunk/octave-forge/extra/control-devel/inst/arx.m 2012-05-07 12:39:54 UTC (rev 10373)
@@ -22,13 +22,16 @@
Y = dat.y{1};
U = dat.u{1};
Ts = dat.tsam{1};
+
+ n = max (na, nb);
## avoid warning: toeplitz: column wins anti-diagonal conflict
PhiY = toeplitz (Y(1:end-1, :), [Y(1, :); zeros(na-1, 1)]);
PhiU = toeplitz (U(1:end-1, :), [U(1, :); zeros(nb-1, 1)]);
Phi = [-PhiY, PhiU];
+ Phi = Phi(n:end, :)
- ## Theta = Phi \ Y(2:end, :); # naive formula
+ ## Theta = Phi \ Y(n+1:end, :); # naive formula
## solve linear least squares problem by pseudoinverse
## the pseudoinverse is computed by singular value decomposition
@@ -42,7 +45,7 @@
V = V(:, 1:r);
S = S(1:r);
U = U(:, 1:r);
- Theta = V * (S .\ (U' * Y(2:end, :))); # U' is the conjugate transpose
+ Theta = V * (S .\ (U' * Y(n+1:end, :))); # U' is the conjugate transpose
A = [1; Theta(1:na)]; # a0 = 1, a1 = Theta(1), an = Theta(n)
B = [0; Theta(na+1:end)]; # b0 = 0 (leading zero required by filt)
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <par...@us...> - 2012-05-07 13:22:07
|
Revision: 10374
http://octave.svn.sourceforge.net/octave/?rev=10374&view=rev
Author: paramaniac
Date: 2012-05-07 13:22:01 +0000 (Mon, 07 May 2012)
Log Message:
-----------
control-devel: support miso arx identification
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-07 12:39:54 UTC (rev 10373)
+++ trunk/octave-forge/extra/control-devel/inst/arx.m 2012-05-07 13:22:01 UTC (rev 10374)
@@ -15,20 +15,24 @@
if (! is_real_scalar (na, nb))
error ("arx: ");
## Test for integers
+ ## numel (nb) == size (dat, 3)
endif
## TODO: handle MIMO and multi-experiment data
+ [~, p, m, ex] = size (dat);
Y = dat.y{1};
U = dat.u{1};
Ts = dat.tsam{1};
- n = max (na, nb);
+ max_nb = max (nb);
+ n = max (na, max_nb);
## avoid warning: toeplitz: column wins anti-diagonal conflict
PhiY = toeplitz (Y(1:end-1, :), [Y(1, :); zeros(na-1, 1)]);
- PhiU = toeplitz (U(1:end-1, :), [U(1, :); zeros(nb-1, 1)]);
- Phi = [-PhiY, PhiU];
+ ## 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, :)
## Theta = Phi \ Y(n+1:end, :); # naive formula
@@ -48,8 +52,15 @@
Theta = V * (S .\ (U' * Y(n+1:end, :))); # U' is the conjugate transpose
A = [1; Theta(1:na)]; # a0 = 1, a1 = Theta(1), an = Theta(n)
- B = [0; Theta(na+1:end)]; # b0 = 0 (leading zero required by filt)
+ ## B = [0; Theta(na+1:end)]; # b0 = 0 (leading zero required by filt)
+ ThetaB = Theta(na+1:end);
+ B = mat2cell (ThetaB, nb);
+ B = reshape (B, 1, []);
+ B = cellfun (@(x) [0; x], B, "uniformoutput", false);
- sys = filt ({B, 1}, {A, A}, Ts);
+ ## sys = filt ({B, 1}, {A, A}, Ts);
+ num = [B, {1}];
+ den = repmat ({A}, 1, m+1);
+ sys = filt (num, den, Ts);
endfunction
\ No newline at end of file
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <par...@us...> - 2012-05-08 13:13:49
|
Revision: 10381
http://octave.svn.sourceforge.net/octave/?rev=10381&view=rev
Author: paramaniac
Date: 2012-05-08 13:13:38 +0000 (Tue, 08 May 2012)
Log Message:
-----------
control-devel: mimo support for arx
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-08 12:06:50 UTC (rev 10380)
+++ trunk/octave-forge/extra/control-devel/inst/arx.m 2012-05-08 13:13:38 UTC (rev 10381)
@@ -22,46 +22,96 @@
## TODO: handle MIMO and multi-experiment data
[~, p, m, ex] = size (dat);
- Y = dat.y{1};
- U = dat.u{1};
+ Y = dat.y;
+ U = dat.u;
Ts = dat.tsam{1};
max_nb = max (nb);
n = max (na, max_nb);
+
+ num = cell (p, m+1);
+ den = cell (p, m+1);
+
+ for i = 1 : p # for every output
+ Phi = cell (1, ex);
+ 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, :);
+ 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);
+ B = reshape (B, 1, []);
+ B = cellfun (@(x) [0; x], B, "uniformoutput", false);
+
+ num(i, :) = [B, {1}];
+ den(i, :) = repmat ({A}, 1, m+1);
+ endfor
+
+ sys = filt (num, den, Ts);
+
+endfunction
+
+
+function theta = __theta__ (phi, y, i, n)
+
+ ## Theta = Phi \ Y(n+1:end, :); # naive formula
+
+ ex = numel (phi); # number of experiments
+
+ if (ex == 1)
+ ## single-experiment dataset
+
+ ## solve linear least squares problem by pseudoinverse
+ ## the pseudoinverse is computed by singular value decomposition
+ ## M = U S V* ---> M+ = V S+ U*
+ ## Th = Ph \ Y = Ph+ Y
+ ## Th = V S+ U* Y, S+ = 1 ./ diag (S)
+
+ [U, S, V] = svd (phi{1}, 0); # 0 for "economy size" decomposition, U overwrites input U
+ S = diag (S); # extract main diagonal
+ r = sum (S > eps*S(1));
+ 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
+ else
+ ## multi-experiment dataset
+ ## TODO: find more sophisticated formula than
+ ## Theta = (Phi1' Phi + Phi2' Phi2 + ...) \ (Phi1' Y1 + Phi2' Y2 + ...)
+
+ ## covariance matrix C = (Phi1' Phi + Phi2' Phi2 + ...)
+ tmp = cellfun (@(Phi) Phi.' * Phi, phi, "uniformoutput", false);
+ C = plus (tmp{:});
+
+ ## PhiTY = (Phi1' Y1 + Phi2' Y2 + ...)
+ tmp = cellfun (@(Phi, Y) Phi.' * Y(n+1:end, i), phi, y, "uniformoutput", false);
+ PhiTY = plus (tmp{:});
+
+ ## pseudoinverse Theta = C \ Phi'Y
+ theta = C \ PhiTY;
+ 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, :);
-
- ## Theta = Phi \ Y(n+1:end, :); # naive formula
-
- ## solve linear least squares problem by pseudoinverse
- ## the pseudoinverse is computed by singular value decomposition
- ## M = U S V* ---> M+ = V S+ U*
- ## Th = Ph \ Y = Ph+ Y
- ## Th = V S+ U* Y, S+ = 1 ./ diag (S)
-
- [U, S, V] = svd (Phi, 0); # 0 for "economy size" decomposition, U overwrites input U
- S = diag (S); # extract main diagonal
- r = sum (S > eps*S(1));
- V = V(:, 1:r);
- S = S(1:r);
- U = U(:, 1:r);
- Theta = V * (S .\ (U' * Y(n+1:end, :))); # U' is the conjugate transpose
-
- A = [1; Theta(1:na)]; # a0 = 1, a1 = Theta(1), an = Theta(n)
- ## B = [0; Theta(na+1:end)]; # b0 = 0 (leading zero required by filt)
- ThetaB = Theta(na+1:end);
- B = mat2cell (ThetaB, nb);
- B = reshape (B, 1, []);
- B = cellfun (@(x) [0; x], B, "uniformoutput", false);
-
- ## sys = filt ({B, 1}, {A, A}, Ts);
- num = [B, {1}];
- den = repmat ({A}, 1, m+1);
- sys = filt (num, den, Ts);
-endfunction
\ No newline at end of file
+endfunction
+%}
\ No newline at end of file
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <par...@us...> - 2012-05-14 19:29:38
|
Revision: 10439
http://octave.svn.sourceforge.net/octave/?rev=10439&view=rev
Author: paramaniac
Date: 2012-05-14 19:29:32 +0000 (Mon, 14 May 2012)
Log Message:
-----------
control-devel: warning for rank-deficient coefficient matrix
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-14 17:23:29 UTC (rev 10438)
+++ trunk/octave-forge/extra/control-devel/inst/arx.m 2012-05-14 19:29:32 UTC (rev 10439)
@@ -112,6 +112,11 @@
[U, S, V] = svd (phi{1}, 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");
+ warning ("sampling time too small");
+ warning ("persistence of excitation");
+ endif
V = V(:, 1:r);
S = S(1:r);
U = U(:, 1:r);
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <par...@us...> - 2012-05-15 06:43:02
|
Revision: 10440
http://octave.svn.sourceforge.net/octave/?rev=10440&view=rev
Author: paramaniac
Date: 2012-05-15 06:42:56 +0000 (Tue, 15 May 2012)
Log Message:
-----------
control-devel: add todo
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-14 19:29:32 UTC (rev 10439)
+++ trunk/octave-forge/extra/control-devel/inst/arx.m 2012-05-15 06:42:56 UTC (rev 10440)
@@ -128,6 +128,7 @@
## covariance matrix C = (Phi1' Phi + Phi2' Phi2 + ...)
tmp = cellfun (@(Phi) Phi.' * Phi, phi, "uniformoutput", false);
+ rc = cellfun (@rcond, tmp); # C auch noch testen? QR oder SVD?
C = plus (tmp{:});
## PhiTY = (Phi1' Y1 + Phi2' Y2 + ...)
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <par...@us...> - 2012-05-15 14:40:21
|
Revision: 10444
http://octave.svn.sourceforge.net/octave/?rev=10444&view=rev
Author: paramaniac
Date: 2012-05-15 14:40:12 +0000 (Tue, 15 May 2012)
Log Message:
-----------
control-devel: style fix
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-15 14:15:20 UTC (rev 10443)
+++ trunk/octave-forge/extra/control-devel/inst/arx.m 2012-05-15 14:40:12 UTC (rev 10444)
@@ -97,10 +97,8 @@
function theta = __theta__ (phi, y, i, n)
## Theta = Phi \ Y(n+1:end, :); # naive formula
-
- ex = numel (phi); # number of experiments
-
- if (ex == 1)
+
+ if (numel (phi) == 1)
## single-experiment dataset
theta = __ls_svd__ (phi{1}, y{1}(n(i)+1:end, i));
else
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <par...@us...> - 2012-05-15 15:22:10
|
Revision: 10445
http://octave.svn.sourceforge.net/octave/?rev=10445&view=rev
Author: paramaniac
Date: 2012-05-15 15:22:00 +0000 (Tue, 15 May 2012)
Log Message:
-----------
control-devel: arx: use Ljung's QR method
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-15 14:40:12 UTC (rev 10444)
+++ trunk/octave-forge/extra/control-devel/inst/arx.m 2012-05-15 15:22:00 UTC (rev 10445)
@@ -95,14 +95,17 @@
function theta = __theta__ (phi, y, i, n)
-
- ## Theta = Phi \ Y(n+1:end, :); # naive formula
- if (numel (phi) == 1)
- ## single-experiment dataset
- theta = __ls_svd__ (phi{1}, y{1}(n(i)+1:end, i));
- else
- ## multi-experiment dataset
+ if (numel (phi) == 1) # single-experiment dataset
+ A = horzcat (phi{1}, y{1}(n(i)+1:end, i)); # [Phi, Y]
+ R0 = triu (qr (A, 0)); # 0 for economy-size R (without zero rows)
+ R1 = R0(1:end-1, 1:end-1); # R1 is triangular - can we exploit this in R1\R2
+ R2 = R0(1:end-1, end);
+ theta = __ls_svd__ (R1, R2); # R1 \ R2
+
+ ## Theta = Phi \ Y(n+1:end, :); # naive formula
+ ## theta = __ls_svd__ (phi{1}, y{1}(n(i)+1:end, i));
+ else # multi-experiment dataset
## TODO: find more sophisticated formula than
## Theta = (Phi1' Phi + Phi2' Phi2 + ...) \ (Phi1' Y1 + Phi2' Y2 + ...)
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <par...@us...> - 2012-05-15 20:01:27
|
Revision: 10446
http://octave.svn.sourceforge.net/octave/?rev=10446&view=rev
Author: paramaniac
Date: 2012-05-15 20:01:20 +0000 (Tue, 15 May 2012)
Log Message:
-----------
control-devel: minor changes
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-15 15:22:00 UTC (rev 10445)
+++ trunk/octave-forge/extra/control-devel/inst/arx.m 2012-05-15 20:01:20 UTC (rev 10446)
@@ -26,6 +26,8 @@
function sys = arx (dat, na, nb)
+ ## TODO: delays
+
if (nargin != 3)
print_usage ();
endif
@@ -97,9 +99,10 @@
function theta = __theta__ (phi, y, i, n)
if (numel (phi) == 1) # single-experiment dataset
+ ## use "square-root algorithm"
A = horzcat (phi{1}, y{1}(n(i)+1:end, i)); # [Phi, Y]
R0 = triu (qr (A, 0)); # 0 for economy-size R (without zero rows)
- R1 = R0(1:end-1, 1:end-1); # R1 is triangular - can we exploit this in R1\R2
+ R1 = R0(1:end-1, 1:end-1); # R1 is triangular - can we exploit this in R1\R2?
R2 = R0(1:end-1, end);
theta = __ls_svd__ (R1, R2); # R1 \ R2
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <par...@us...> - 2012-05-15 20:40:36
|
Revision: 10447
http://octave.svn.sourceforge.net/octave/?rev=10447&view=rev
Author: paramaniac
Date: 2012-05-15 20:40:30 +0000 (Tue, 15 May 2012)
Log Message:
-----------
control-devel: minor changes
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-15 20:01:20 UTC (rev 10446)
+++ trunk/octave-forge/extra/control-devel/inst/arx.m 2012-05-15 20:40:30 UTC (rev 10447)
@@ -130,6 +130,10 @@
function x = __ls_svd__ (A, b)
+ ## solve the problem Ax=b
+ ## x = A\b would also work,
+ ## but this way we have better control and warnings
+
## solve linear least squares problem by pseudoinverse
## the pseudoinverse is computed by singular value decomposition
## M = U S V* ---> M+ = V S+ U*
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <par...@us...> - 2012-05-16 09:02:51
|
Revision: 10449
http://octave.svn.sourceforge.net/octave/?rev=10449&view=rev
Author: paramaniac
Date: 2012-05-16 09:02:41 +0000 (Wed, 16 May 2012)
Log Message:
-----------
control-devel: add fixme comments to arx
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-16 08:00:43 UTC (rev 10448)
+++ trunk/octave-forge/extra/control-devel/inst/arx.m 2012-05-16 09:02:41 UTC (rev 10449)
@@ -27,6 +27,13 @@
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 ();
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <par...@us...> - 2012-05-20 09:18:59
|
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.
|
|
From: <par...@us...> - 2012-05-22 15:43:57
|
Revision: 10496
http://octave.svn.sourceforge.net/octave/?rev=10496&view=rev
Author: paramaniac
Date: 2012-05-22 15:43:48 +0000 (Tue, 22 May 2012)
Log Message:
-----------
control-devel: improve code documentation
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-22 15:24:41 UTC (rev 10495)
+++ trunk/octave-forge/extra/control-devel/inst/arx.m 2012-05-22 15:43:48 UTC (rev 10496)
@@ -64,12 +64,14 @@
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)
-
+ ## create empty cells for numerator and denominator polynomials
num = cell (p, m+p);
den = cell (p, m+p);
-
+
+ ## MIMO (p-by-m) models are identified as p MISO (1-by-m) models
+ ## For multi-experiment data, minimize the trace of the error
for i = 1 : p # for every output
- Phi = cell (ex, 1);
+ Phi = cell (ex, 1); # one regression matrix per 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)
@@ -78,19 +80,24 @@
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
-
+
+ ## compute parameter vector Theta
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)
- B = mat2cell (ThetaB, nb(i,:));
- B = reshape (B, 1, []);
- B = cellfun (@(x) [0; x], B, "uniformoutput", false);
- Be = repmat ({0}, 1, p);
- Be(i) = 1;
- num(i, :) = [B, Be];
- den(i, :) = repmat ({A}, 1, m+p);
+ ## extract polynomial matrices A and B from Theta
+ ## A is a scalar polynomial for output i, i=1:p
+ ## B is polynomial row vector (1-by-m) for output i
+ A = [1; Theta(1:na(i))]; # a0 = 1, a1 = Theta(1), an = Theta(n)
+ ThetaB = Theta(na(i)+1:end); # all polynomials from B are in one column vector
+ B = mat2cell (ThetaB, nb(i,:)); # now separate the polynomials, one for each input
+ B = reshape (B, 1, []); # make B a row cell (1-by-m)
+ B = cellfun (@(x) [0; x], B, "uniformoutput", false); # b0 = 0 (leading zero required by filt)
+
+ ## add error inputs
+ Be = repmat ({0}, 1, p); # there are as many error inputs as system outputs (p)
+ Be(i) = 1; # inputs m+1:m+p are zero, except m+i which is one
+ num(i, :) = [B, Be]; # numerator polynomials for output i, individual for each input
+ den(i, :) = repmat ({A}, 1, m+p); # in a row (output i), all inputs have the same denominator polynomial
endfor
## A(q) y(t) = B(q) u(t) + e(t)
@@ -100,8 +107,10 @@
## 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);
-
+ sys = filt (num, den, tsam); # filt creates a transfer function in z^-1
+
+ ## compute initial state vector x0 if requested
+ ## this makes only sense for state-space models, therefore convert TF to SS
if (nargout > 1)
sys = prescale (ss (sys(:,1:m)));
x0 = slib01cd (Y, U, sys.a, sys.b, sys.c, sys.d, 0.0)
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
×
Want the latest updates on software, tech news, and AI?
Get latest updates about software, tech news, and AI from SourceForge directly in your inbox once a month.
Thanks for helping keep SourceForge clean.
X