Menu
▾
▴
octave-cvsupdate
|
From: <par...@us...> - 2012-05-01 20:31:13
|
Revision: 10351
http://octave.svn.sourceforge.net/octave/?rev=10351&view=rev
Author: paramaniac
Date: 2012-05-01 20:31:07 +0000 (Tue, 01 May 2012)
Log Message:
-----------
control-devel: compute pseudoinverse by SVD instead of \ operator
Modified Paths:
--------------
trunk/octave-forge/extra/control-devel/inst/arx.m
Added Paths:
-----------
trunk/octave-forge/extra/control-devel/devel/test_arx.m
Removed Paths:
-------------
trunk/octave-forge/extra/control-devel/inst/test_arx.m
Copied: trunk/octave-forge/extra/control-devel/devel/test_arx.m (from rev 10337, trunk/octave-forge/extra/control-devel/inst/test_arx.m)
===================================================================
--- trunk/octave-forge/extra/control-devel/devel/test_arx.m (rev 0)
+++ trunk/octave-forge/extra/control-devel/devel/test_arx.m 2012-05-01 20:31:07 UTC (rev 10351)
@@ -0,0 +1,10 @@
+u = [ 0; 0.5; 1; 1; 1; 1; 1 ];
+y = [ 0; 0; 0.25; 0.62; 0.81; 0.90; 0.95 ];
+
+dat = iddata (y, u)
+
+sys = arx (dat, 1, 1)
+
+
+ysim = lsim (sys(1,1), u);
+
Modified: trunk/octave-forge/extra/control-devel/inst/arx.m
===================================================================
--- trunk/octave-forge/extra/control-devel/inst/arx.m 2012-05-01 20:26:49 UTC (rev 10350)
+++ trunk/octave-forge/extra/control-devel/inst/arx.m 2012-05-01 20:31:07 UTC (rev 10351)
@@ -28,8 +28,22 @@
PhiU = toeplitz (U(1:end-1, :), [U(1, :); zeros(nb-1, 1)]);
Phi = [-PhiY, PhiU];
- Theta = Phi \ Y(2:end, :);
+ 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)
+
+ [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(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)
Deleted: trunk/octave-forge/extra/control-devel/inst/test_arx.m
===================================================================
--- trunk/octave-forge/extra/control-devel/inst/test_arx.m 2012-05-01 20:26:49 UTC (rev 10350)
+++ trunk/octave-forge/extra/control-devel/inst/test_arx.m 2012-05-01 20:31:07 UTC (rev 10351)
@@ -1,10 +0,0 @@
-u = [ 0; 0.5; 1; 1; 1; 1; 1 ];
-y = [ 0; 0; 0.25; 0.62; 0.81; 0.90; 0.95 ];
-
-dat = iddata (y, u)
-
-sys = arx (dat, 1, 1)
-
-
-ysim = lsim (sys(1,1), u);
-
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <par...@us...> - 2012-05-19 08:39:11
|
Revision: 10464
http://octave.svn.sourceforge.net/octave/?rev=10464&view=rev
Author: paramaniac
Date: 2012-05-19 08:39:04 +0000 (Sat, 19 May 2012)
Log Message:
-----------
control-devel: fix dimension conflict for multi-experiment ARX identification
Modified Paths:
--------------
trunk/octave-forge/extra/control-devel/inst/arx.m
Added Paths:
-----------
trunk/octave-forge/extra/control-devel/devel/DestillationMEarx.m
trunk/octave-forge/extra/control-devel/devel/LakeErieARX.m
Added: trunk/octave-forge/extra/control-devel/devel/DestillationMEarx.m
===================================================================
--- trunk/octave-forge/extra/control-devel/devel/DestillationMEarx.m (rev 0)
+++ trunk/octave-forge/extra/control-devel/devel/DestillationMEarx.m 2012-05-19 08:39:04 UTC (rev 10464)
@@ -0,0 +1,92 @@
+%{
+This file describes the data in the destill.dat file.
+1. Contributed by:
+ Peter Van Overschee
+ K.U.Leuven - ESAT - SISTA
+ K. Mercierlaan 94
+ 3001 Heverlee
+ Pet...@es...
+2. Process/Description:
+ Data of a simulation (not real !) related to the identification
+ of an ethane-ethylene destillationcolumn. The series consists of 4
+ series:
+ U_dest, Y_dest: without noise (original series)
+ U_dest_n10, Y_dest_n10: 10 percent additive white noise
+ U_dest_n20, Y_dest_n20: 20 percent additive white noise
+ U_dest_n30, Y_dest_n30: 30 percent additive white noise
+3. Sampling time
+ 15 min.
+4. Number of samples:
+ 90 samples
+5. Inputs:
+ a. ratio between the reboiler duty and the feed flow
+ b. ratio between the reflux rate and the feed flow
+ c. ratio between the distillate and the feed flow
+ d. input ethane composition
+ e. top pressure
+6. Outputs:
+ a. top ethane composition
+ b. bottom ethylene composition
+ c. top-bottom differential pressure.
+7. References:
+ R.P. Guidorzi, M.P. Losito, T. Muratori, The range error test in the
+ structural identification of linear multivariable systems,
+ IEEE transactions on automatic control, Vol AC-27, pp 1044-1054, oct.
+ 1982.
+8. Known properties/peculiarities
+
+9. Some MATLAB-code to retrieve the data
+ !gunzip destill.dat.Z
+ load destill.dat
+ U=destill(:,1:20);
+ Y=destill(:,21:32);
+ U_dest=U(:,1:5);
+ U_dest_n10=U(:,6:10);
+ U_dest_n20=U(:,11:15);
+ U_dest_n30=U(:,16:20);
+ Y_dest=Y(:,1:3);
+ Y_dest_n10=Y(:,4:6);
+ Y_dest_n20=Y(:,7:9);
+ Y_dest_n30=Y(:,10:12);
+%}
+
+clear all, close all, clc
+
+% DaISy code is wrong,
+% first column is sample number
+load destill.dat
+U=destill(:,2:21);
+Y=destill(:,22:33);
+U_dest=U(:,1:5);
+U_dest_n10=U(:,6:10);
+U_dest_n20=U(:,11:15);
+U_dest_n30=U(:,16:20);
+Y_dest=Y(:,1:3);
+Y_dest_n10=Y(:,4:6);
+Y_dest_n20=Y(:,7:9);
+Y_dest_n30=Y(:,10:12);
+
+Y = {Y_dest; Y_dest_n10; Y_dest_n20; Y_dest_n30};
+U = {U_dest; U_dest_n10; U_dest_n20; U_dest_n30};
+
+dat = iddata (Y, U)
+
+[sys, x0] = ident (dat, 5, 4) % s=5, n=4
+sys2 = arx (dat, 4, 4);
+
+x0=x0{1};
+
+[y, t] = lsim (sys, U_dest, [], x0);
+[y2, t2] = lsim (sys(:, 1:5), U_dest);
+
+% ARX has no initial conditions, therefore the bad results
+
+err = norm (Y_dest - y, 1) / norm (Y_dest, 1)
+err2 = norm (Y_dest - y2, 1) / norm (Y_dest, 1)
+
+figure (1)
+%plot (t, Y_dest, 'b')
+plot (t, Y_dest, t, y, t, y2)
+legend ('y measured', 'y MOEN4', 'y ARX', 'location', 'southeast')
+
+
Added: trunk/octave-forge/extra/control-devel/devel/LakeErieARX.m
===================================================================
--- trunk/octave-forge/extra/control-devel/devel/LakeErieARX.m (rev 0)
+++ trunk/octave-forge/extra/control-devel/devel/LakeErieARX.m 2012-05-19 08:39:04 UTC (rev 10464)
@@ -0,0 +1,100 @@
+%{
+This file describes the data in the erie.dat file.
+1. Contributed by:
+ Peter Van Overschee
+ K.U.Leuven - ESAT - SISTA
+ K. Mercierlaan 94
+ 3001 Heverlee
+ Pet...@es...
+2. Process/Description:
+ Data of a simulation (not real !) related to the related to the
+ identification of the western basin of Lake Erie. The series consists
+ of 4 series:
+ U_erie, Y_erie: without noise (original series)
+ U_erie_n10, Y_erie_n10: 10 percent additive white noise
+ U_erie_n20, Y_erie_n20: 20 percent additive white noise
+ U_erie_n30, Y_erie_n30: 30 percent additive white noise
+3. Sampling time
+ 1 month
+4. Number of samples:
+ 57 samples
+5. Inputs:
+ a. water temperature
+ b. water conductivity
+ c. water alkalinity
+ d. NO3
+ e. total hardness
+6. Outputs:
+ a. dissolved oxigen
+ b. algae
+7. References:
+ R.P. Guidorzi, M.P. Losito, T. Muratori, On the last eigenvalue
+ test in the structural identification of linear multivariable
+ systems, Proceedings of the V European meeting on cybernetics and
+ systems research, Vienna, april 1980.
+8. Known properties/peculiarities
+ The considered period runs from march 1968 till november 1972.
+9. Some MATLAB-code to retrieve the data
+ !guzip erie.dat.Z
+ load erie.dat
+ U=erie(:,1:20);
+ Y=erie(:,21:28);
+ U_erie=U(:,1:5);
+ U_erie_n10=U(:,6:10);
+ U_erie_n20=U(:,11:15);
+ U_erie_n30=U(:,16:20);
+ Y_erie=Y(:,1:2);
+ Y_erie_n10=Y(:,3:4);
+ Y_erie_n20=Y(:,5:6);
+ Y_erie_n30=Y(:,7:8);
+%}
+
+clear all, close all, clc
+
+% DaISy code is wrong,
+% first column is sample number
+load erie.dat
+U=erie(:,2:21);
+Y=erie(:,22:29);
+U_erie=U(:,1:5);
+U_erie_n10=U(:,6:10);
+U_erie_n20=U(:,11:15);
+U_erie_n30=U(:,16:20);
+Y_erie=Y(:,1:2);
+Y_erie_n10=Y(:,3:4);
+Y_erie_n20=Y(:,5:6);
+Y_erie_n30=Y(:,7:8);
+
+Y = {Y_erie; Y_erie_n10; Y_erie_n20; Y_erie_n30};
+U = {U_erie; U_erie_n10; U_erie_n20; U_erie_n30};
+
+dat = iddata (Y, U, [], 'inname', {'a. water temperature';
+ 'b. water conductivity';
+ 'c. water alkalinity';
+ 'd. NO3';
+ 'e. total hardness'}, \
+ 'outname', {'a. dissolved oxygen';
+ 'b. algae'})
+
+[sys, x0] = ident (dat, 5, 4) % s=5, n=4
+sys2 = arx (dat, 4, 4)
+
+x0=x0{1};
+
+[y, t] = lsim (sys, U_erie, [], x0);
+[y2, t2] = lsim (sys2(:, 1:5), U_erie);
+
+err = norm (Y_erie - y, 1) / norm (Y_erie, 1)
+err2 = norm (Y_erie - y2, 1) / norm (Y_erie, 1)
+
+
+figure (1)
+p = columns (Y_erie);
+for k = 1 : p
+ subplot (2, 1, k)
+ plot (t, Y_erie(:,k), t, y(:,k), t, y2(:,k))
+endfor
+
+legend ('y measured', 'y MOEN4', 'y ARX', 'location', 'southeast')
+
+
Modified: trunk/octave-forge/extra/control-devel/inst/arx.m
===================================================================
--- trunk/octave-forge/extra/control-devel/inst/arx.m 2012-05-18 21:59:47 UTC (rev 10463)
+++ trunk/octave-forge/extra/control-devel/inst/arx.m 2012-05-19 08:39:04 UTC (rev 10464)
@@ -76,7 +76,7 @@
den = cell (p, m+1);
for i = 1 : p # for every output
- Phi = cell (1, ex);
+ Phi = cell (ex, 1);
for e = 1 : ex # for every experiment
## avoid warning: toeplitz: column wins anti-diagonal conflict
## therefore set first row element equal to y(1)
@@ -123,7 +123,7 @@
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 + ...)
tmp = cellfun (@(Phi, Y) Phi.' * Y(n(i)+1:end, i), phi, y, "uniformoutput", false);
PhiTY = plus (tmp{:});
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <par...@us...> - 2012-07-20 20:08:03
|
Revision: 10757
http://octave.svn.sourceforge.net/octave/?rev=10757&view=rev
Author: paramaniac
Date: 2012-07-20 20:07:57 +0000 (Fri, 20 Jul 2012)
Log Message:
-----------
control-devel: minor changes
Modified Paths:
--------------
trunk/octave-forge/extra/control-devel/devel/Destillation.m
trunk/octave-forge/extra/control-devel/devel/HeatingSystem.m
trunk/octave-forge/extra/control-devel/devel/LakeErieARX.m
trunk/octave-forge/extra/control-devel/inst/arx.m
Modified: trunk/octave-forge/extra/control-devel/devel/Destillation.m
===================================================================
--- trunk/octave-forge/extra/control-devel/devel/Destillation.m 2012-07-20 19:37:05 UTC (rev 10756)
+++ trunk/octave-forge/extra/control-devel/devel/Destillation.m 2012-07-20 20:07:57 UTC (rev 10757)
@@ -85,4 +85,12 @@
plot (t, Y_dest, 'b', t, y, 'r')
legend ('y measured', 'y simulated', 'location', 'southeast')
+figure (2)
+p = columns (Y_dest);
+for k = 1 : 3
+ subplot (3, 1, k)
+ plot (t, Y_dest(:,k), 'b', t, y(:,k), 'r')
+ xlim ([0, 90])
+endfor
+legend ('y measured', 'y simulated', 'location', 'southeast')
Modified: trunk/octave-forge/extra/control-devel/devel/HeatingSystem.m
===================================================================
--- trunk/octave-forge/extra/control-devel/devel/HeatingSystem.m 2012-07-20 19:37:05 UTC (rev 10756)
+++ trunk/octave-forge/extra/control-devel/devel/HeatingSystem.m 2012-07-20 20:07:57 UTC (rev 10757)
@@ -93,7 +93,7 @@
xlim ([t(1), t(end)])
xlabel ('Time [s]')
ylabel ('Temperature [Degree Celsius]')
-legend ('measured', 'simulated subspace', 'simulated ARX', 'location', 'southeast')
+legend ('measurement DaISy', 'simulation MOEN4', 'simulation ARX', 'location', 'northeast')
Modified: trunk/octave-forge/extra/control-devel/devel/LakeErieARX.m
===================================================================
--- trunk/octave-forge/extra/control-devel/devel/LakeErieARX.m 2012-07-20 19:37:05 UTC (rev 10756)
+++ trunk/octave-forge/extra/control-devel/devel/LakeErieARX.m 2012-07-20 20:07:57 UTC (rev 10757)
@@ -25,7 +25,7 @@
d. NO3
e. total hardness
6. Outputs:
- a. dissolved oxigen
+ a. dissolved oxygen
b. algae
7. References:
R.P. Guidorzi, M.P. Losito, T. Muratori, On the last eigenvalue
@@ -76,7 +76,7 @@
'outname', {'a. dissolved oxygen';
'b. algae'})
-[sys, x0] = moen4 (dat, 's', 5, 'n', 4) % s=5, n=4
+[sys, x0, info] = moen4 (dat, 's', 5, 'n', 4) % s=5, n=4
% sys2 = arx (dat, 4, 4)
[sys2, x02] = arx (dat, 4)
@@ -97,8 +97,59 @@
for k = 1 : p
subplot (2, 1, k)
plot (t, Y_erie(:,k), t, y(:,k), t, y2(:,k))
+ grid on
endfor
-legend ('y measured', 'y MOEN4', 'y ARX', 'location', 'southeast')
+subplot (2, 1, 1)
+title ('DaISy: Lake Erie [96-005]')
+ylabel ('Dissolved Oxygen [n.s.]')
+xlim ([0, 56])
+subplot (2, 1, 2)
+ylabel ('Algae [n.s.]')
+xlabel ('Time [months]')
+xlim ([0, 56])
+legend ('measurement DaISy', 'simulation MOEN4', 'simulation ARX', 'location', 'northeast')
+
+
+
+
+
+l = lqe (sys, info.Q, 100*info.Ry)
+
+
+
+[a, b, c, d] = ssdata (sys);
+
+sys2 = ss ([a-l*c], [b-l*d, l], c, [d, zeros(2)], -1)
+
+[sys, ~, info] = moen4 (dat, 's', 5, 'n', 4, 'noise', 'k')
+
+[y, t] = lsim (sys, [U_erie, Y_erie], [], x0);
+[y2, t2] = lsim (sys2, [U_erie, Y_erie], [], x0);
+
+errkp = norm (Y_erie - y, 1) / norm (Y_erie, 1)
+err2kp = norm (Y_erie - y2, 1) / norm (Y_erie, 1)
+
+figure (2)
+p = columns (Y_erie);
+for k = 1 : p
+ subplot (2, 1, k)
+ plot (t, Y_erie(:,k), t, y(:,k), t, y2(:,k))
+ grid on
+endfor
+
+subplot (2, 1, 1)
+title ('DaISy: Lake Erie [96-005]')
+ylabel ('Dissolved Oxygen [n.s.]')
+xlim ([0, 56])
+
+subplot (2, 1, 2)
+ylabel ('Algae [n.s.]')
+xlabel ('Time [months]')
+xlim ([0, 56])
+
+legend ('measurement DaISy', 'Kalman Predictor', 'Kalman Predictor weak', 'location', 'northeast')
+
+
Modified: trunk/octave-forge/extra/control-devel/inst/arx.m
===================================================================
--- trunk/octave-forge/extra/control-devel/inst/arx.m 2012-07-20 19:37:05 UTC (rev 10756)
+++ trunk/octave-forge/extra/control-devel/inst/arx.m 2012-07-20 20:07:57 UTC (rev 10757)
@@ -91,11 +91,11 @@
error ("arx: first argument must be a time-domain iddata dataset");
endif
- if (is_real_scalar (varargin{1})) # arx (dat, n, ...)
+ if (is_real_scalar (varargin{1})) # arx (dat, n, ...)
varargin = horzcat (varargin(2:end), {"na"}, varargin(1), {"nb"}, varargin(1));
endif
- if (isstruct (varargin{1})) # arx (dat, opt, ...), arx (dat, n, opt, ...)
+ if (isstruct (varargin{1})) # arx (dat, opt, ...), arx (dat, n, opt, ...)
varargin = horzcat (__opt2cell__ (varargin{1}), varargin(2:end));
endif
@@ -232,7 +232,7 @@
## 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 + ...)
+ ## Theta = (Phi1' Phi1 + Phi2' Phi2 + ...) \ (Phi1' Y1 + Phi2' Y2 + ...)
## covariance matrix C = (Phi1' Phi + Phi2' Phi2 + ...)
tmp = cellfun (@(Phi) Phi.' * Phi, phi, "uniformoutput", false);
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