Mercurial > forge

changeset 10278:5a5f3499e34a octave-forge

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
control-devel: fix dimension conflict for multi-experiment ARX identification
author paramaniac
date Sat, 19 May 2012 08:39:04 +0000
parents 115f015abea5
children df4726eaff1e
files extra/control-devel/devel/DestillationMEarx.m extra/control-devel/devel/LakeErieARX.m extra/control-devel/inst/arx.m
diffstat 3 files changed, 194 insertions(+), 2 deletions(-) [+]
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/extra/control-devel/devel/DestillationMEarx.m	Sat May 19 08:39:04 2012 +0000
@@ -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
+	Peter.Vanoverschee@esat.kuleuven.ac.be
+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')
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/extra/control-devel/devel/LakeErieARX.m	Sat May 19 08:39:04 2012 +0000
@@ -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
+	Peter.Vanoverschee@esat.kuleuven.ac.be
+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')
+
+
--- a/extra/control-devel/inst/arx.m	Fri May 18 21:59:47 2012 +0000
+++ b/extra/control-devel/inst/arx.m	Sat May 19 08:39:04 2012 +0000
@@ -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{:});