From: <par...@us...> - 2012-03-07 19:09:31
|
Revision: 9774 http://octave.svn.sourceforge.net/octave/?rev=9774&view=rev Author: paramaniac Date: 2012-03-07 19:09:25 +0000 (Wed, 07 Mar 2012) Log Message: ----------- control-devel: touch up cat (2) Modified Paths: -------------- trunk/octave-forge/extra/control-devel/devel/test_iddata.m trunk/octave-forge/extra/control-devel/inst/@iddata/cat.m Modified: trunk/octave-forge/extra/control-devel/devel/test_iddata.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/test_iddata.m 2012-03-07 18:46:36 UTC (rev 9773) +++ trunk/octave-forge/extra/control-devel/devel/test_iddata.m 2012-03-07 19:09:25 UTC (rev 9774) @@ -41,7 +41,9 @@ cat (3, d, e) +%cat (1, b, 4) + un = iddata ({(1:10).', (21:30).'}, {(41:50).', (61:70).'}, [], "expname", strseq ("alpha", 1:2)); vn = iddata ({(11:20).', (31:40).'}, {(51:60).', (71:80).'}, [], "expname", strseq ("beta", 1:2)); n = [un; vn] @@ -49,4 +51,5 @@ cat (1, un, un, vn, vn, vn) -cat (1, b) +%dat = iddata (ones (100, 3)); +%dat2 = cat (1, dat, zeros (4, 3), dat) Modified: trunk/octave-forge/extra/control-devel/inst/@iddata/cat.m =================================================================== --- trunk/octave-forge/extra/control-devel/inst/@iddata/cat.m 2012-03-07 18:46:36 UTC (rev 9773) +++ trunk/octave-forge/extra/control-devel/inst/@iddata/cat.m 2012-03-07 19:09:25 UTC (rev 9774) @@ -37,8 +37,13 @@ endif ## store all datasets in a single struct 'tmp' + ## tmp is not a valid iddata set anymore, + ## but it doesn't matter, we want just a struct tmp = cellfun (@iddata, varargin); [n, p, m, e] = cellfun (@size, varargin, "uniformoutput", false); + + ## TODO: dat = iddata (ones (100, 3)); + ## dat = cat (1, dat, zeros (4, 3), dat) ## default values for metadata ## some of them are overwritten in the switch statement below @@ -152,6 +157,8 @@ mat2str (vertcat (n{:}), 10)); endif + ## TODO: check sampling times + endfunction This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2012-03-13 10:30:05
|
Revision: 9858 http://octave.svn.sourceforge.net/octave/?rev=9858&view=rev Author: paramaniac Date: 2012-03-13 10:29:54 +0000 (Tue, 13 Mar 2012) Log Message: ----------- control-devel: quicksave draft code (2) Modified Paths: -------------- trunk/octave-forge/extra/control-devel/devel/data_ib01ad.m trunk/octave-forge/extra/control-devel/src/slib01ad.cc Modified: trunk/octave-forge/extra/control-devel/devel/data_ib01ad.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/data_ib01ad.m 2012-03-13 09:11:51 UTC (rev 9857) +++ trunk/octave-forge/extra/control-devel/devel/data_ib01ad.m 2012-03-13 10:29:54 UTC (rev 9858) @@ -1,5 +1,6 @@ % IB01AD EXAMPLE PROGRAM DATA % 15 0 1 1 1000 0.0 -1.0 M C N O N N +% nobr, n, m, l, nsmp, rcond, tol, meth, alg, jobd, batch,conct,ctrl U = [ 6.41 Modified: trunk/octave-forge/extra/control-devel/src/slib01ad.cc =================================================================== --- trunk/octave-forge/extra/control-devel/src/slib01ad.cc 2012-03-13 09:11:51 UTC (rev 9857) +++ trunk/octave-forge/extra/control-devel/src/slib01ad.cc 2012-03-13 10:29:54 UTC (rev 9858) @@ -59,88 +59,116 @@ int nargin = args.length (); octave_value_list retval; - if (nargin != 13) + if (nargin != 11) { print_usage (); } else { // arguments in - char dico; + char meth; + char alg; char jobd; - char jobmr; - char jobcf; - char ordsel; + char batch; + char conct; + char ctrl; - Matrix a = args(0).matrix_value (); - Matrix b = args(1).matrix_value (); - Matrix c = args(2).matrix_value (); - Matrix d = args(3).matrix_value (); + Matrix y = args(0).matrix_value (); + Matrix u = args(1).matrix_value (); + int nobr = args(2).int_value (); - const int idico = args(4).int_value (); - int ncr = args(5).int_value (); - const int iordsel = args(6).int_value (); - const int ijobd = args(7).int_value (); - const int ijobmr = args(8).int_value (); - - Matrix f = args(9).matrix_value (); - Matrix g = args(10).matrix_value (); + const int imeth = args(3).int_value (); + const int ialg = args(4).int_value (); + const int ijobd = args(5).int_value (); + const int ibatch = args(6).int_value (); + const int iconct = args(7).int_value (); + const int ictrl = args(8).int_value (); + + double rcond = args(9).double_value (); + double tol = args(10).double_value (); + - const int ijobcf = args(11).int_value (); - double tol = args(12).double_value (); - - if (idico == 0) - dico = 'C'; + if (imeth == 0) + meth = 'M'; else - dico = 'D'; + meth = 'N'; - if (iordsel == 0) - ordsel = 'F'; - else - ordsel = 'A'; - + switch (ialg) + { + case 0: + alg = 'C'; + break; + case 1: + alg = 'F'; + break; + case 2: + alg = 'Q'; + break; + default: + error ("slib01ad: argument 'alg' invalid"); + } + if (ijobd == 0) - jobd = 'Z'; + jobd = 'M'; else - jobd = 'D'; - - if (ijobcf == 0) - jobcf = 'L'; - else - jobcf = 'R'; - - switch (ijobmr) + jobd = 'N'; + + switch (ibatch) { case 0: - jobmr = 'B'; + batch = 'F'; break; case 1: - jobmr = 'F'; + batch = 'I'; break; + case 2: + batch = 'L'; + break; + case 3: + batch = 'O'; + break; default: - error ("slib01ad: argument jobmr invalid"); + error ("slib01ad: argument 'batch' invalid"); } + if (iconct == 0) + conct = 'C'; + else + conct = 'N'; - int n = a.rows (); // n: number of states - int m = b.columns (); // m: number of inputs - int p = c.rows (); // p: number of outputs + if (ictrl == 0) + ctrl = 'C'; + else + ctrl = 'N'; - int lda = max (1, n); - int ldb = max (1, n); - int ldc = max (1, p); - int ldd; + + int m = u.columns (); // m: number of inputs + int l = y.columns (); // l: number of outputs + int nsmp = y.rows (); // nsmp: number of samples + // y.rows == u.rows is checked by iddata class + // TODO: check minimal nsmp size - if (jobd == 'Z') - ldd = 1; - else - ldd = max (1, p); + int ldu; + + if (m == 0) + ldu = 1; + else // m > 0 + ldu = nsmp; - int ldf = max (1, m); - int ldg = max (1, n); + int ldy = nsmp; // arguments out - ColumnVector hsv (n); + int ldr; + + if (meth == 'M' && jobd == 'M') + ldr = max (2*(m+l)*nobr, 3*m*nobr); + else if (meth == 'N' || (meth == 'M' && jobd == 'N')) + ldr = 2*(m+l)*nobr; + else + error ("slib01ad: could not handle 'ldr' case"); + + Matrix r (ldr, 2*(m+l)*nobr); + ColumnVector sv (l*nobr); // workspace int liwork; @@ -220,7 +248,7 @@ "coefficient matrix"}; - error_msg ("ident", info, 6, err_msg); + error_msg ("ident", info, 2, err_msg); warning_msg ("ident", iwarn, 5, warn_msg); This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2012-03-16 18:08:58
|
Revision: 9922 http://octave.svn.sourceforge.net/octave/?rev=9922&view=rev Author: paramaniac Date: 2012-03-16 18:08:52 +0000 (Fri, 16 Mar 2012) Log Message: ----------- control-devel: finish oct-file for n4sid & moesp identification Modified Paths: -------------- trunk/octave-forge/extra/control-devel/src/devel_slicot_functions.cc trunk/octave-forge/extra/control-devel/src/slident.cc Added Paths: ----------- trunk/octave-forge/extra/control-devel/devel/test_slident.m Added: trunk/octave-forge/extra/control-devel/devel/test_slident.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/test_slident.m (rev 0) +++ trunk/octave-forge/extra/control-devel/devel/test_slident.m 2012-03-16 18:08:52 UTC (rev 9922) @@ -0,0 +1,2090 @@ +% IB01AD EXAMPLE PROGRAM DATA +% 15 0 1 1 1000 0.0 -1.0 M C N O N N +% nobr, n, m, l, nsmp, rcond, tol, meth, alg, jobd, batch,conct,ctrl + +% IB01BD EXAMPLE PROGRAM DATA +% 15 0 1 1 1000 0.0 -1.0 C C N O N N A K + + +U = [ + 6.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 6.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 6.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 6.41 + 3.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 3.41 + 6.41 + 3.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 6.41 + 3.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 6.41 + 3.41 + 3.41 + 3.41 + 3.41 + 6.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 3.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 6.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 6.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 6.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 6.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 6.41 + 3.41 + 3.41 + 3.41 + 6.41 + 3.41 + 3.41 + 3.41 + 6.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 6.41 + 6.41 + 3.41 + 6.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 3.41 + 3.41 + 6.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 6.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 3.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 + 6.41 +]; + +Y = [ + 4.766099 + 4.763659 + 4.839359 + 5.002979 + 5.017629 + 5.056699 + 5.154379 + 5.361949 + 5.425439 + 5.569519 + 5.681849 + 5.742899 + 5.803949 + 5.918729 + 5.821049 + 5.447419 + 5.061589 + 4.629349 + 4.267939 + 4.011519 + 3.850349 + 3.711159 + 3.569519 + 3.518239 + 3.652549 + 3.818609 + 3.862559 + 4.011519 + 4.353409 + 4.705049 + 5.083559 + 5.344859 + 5.274039 + 5.127519 + 4.761219 + 4.451089 + 4.221539 + 4.045709 + 3.874769 + 3.730689 + 3.662319 + 3.576849 + 3.542659 + 3.479169 + 3.454749 + 3.359509 + 3.298459 + 3.225199 + 3.200779 + 3.225199 + 3.227639 + 3.274039 + 3.457189 + 3.867449 + 4.321659 + 4.492599 + 4.431549 + 4.243519 + 4.050599 + 3.857679 + 3.730689 + 3.791739 + 3.921169 + 3.955359 + 3.847909 + 3.725809 + 3.611039 + 3.716039 + 4.092109 + 4.480389 + 4.814939 + 5.054259 + 5.303339 + 5.486489 + 5.672089 + 5.779529 + 5.799069 + 5.664759 + 5.291129 + 4.880879 + 4.558529 + 4.184909 + 3.889419 + 3.708719 + 3.623249 + 3.569519 + 3.718479 + 4.033499 + 4.412009 + 4.629349 + 4.558529 + 4.394919 + 4.180019 + 4.197119 + 4.431549 + 4.714819 + 4.961459 + 5.300899 + 5.567079 + 5.681849 + 5.545099 + 5.188569 + 4.883319 + 4.600049 + 4.270379 + 4.038389 + 3.838139 + 3.711159 + 3.591499 + 3.535329 + 3.486489 + 3.476729 + 3.425439 + 3.381489 + 3.369279 + 3.364389 + 3.347299 + 3.381489 + 3.420559 + 3.413229 + 3.452309 + 3.635459 + 4.038389 + 4.375379 + 4.727029 + 5.056699 + 5.298459 + 5.532889 + 5.466959 + 5.195899 + 4.885759 + 4.763659 + 4.875989 + 5.042049 + 5.283809 + 5.491379 + 5.596379 + 5.672089 + 5.772209 + 5.830819 + 5.933379 + 5.899189 + 5.935819 + 5.894309 + 5.918729 + 5.994429 + 5.957799 + 6.031059 + 6.062809 + 6.040829 + 6.096999 + 6.123859 + 6.162929 + 6.040829 + 5.845469 + 5.772209 + 5.799069 + 5.923609 + 5.928499 + 6.001759 + 6.001759 + 6.060369 + 5.882099 + 5.510909 + 5.322879 + 5.371719 + 5.454749 + 5.437649 + 5.159269 + 4.902859 + 4.587839 + 4.502369 + 4.595159 + 4.824709 + 5.064029 + 5.271599 + 5.466959 + 5.615919 + 5.528009 + 5.254499 + 4.883319 + 4.517019 + 4.197119 + 4.001759 + 3.806399 + 3.904079 + 3.923609 + 3.869889 + 3.806399 + 3.720929 + 3.818609 + 4.140949 + 4.529229 + 4.805179 + 5.086009 + 5.339969 + 5.532889 + 5.576849 + 5.667199 + 5.791739 + 5.850349 + 5.923609 + 5.921169 + 5.977339 + 5.740459 + 5.388809 + 5.000539 + 4.849129 + 4.944369 + 5.173919 + 5.369279 + 5.447419 + 5.603709 + 5.730689 + 5.850349 + 5.979779 + 5.991989 + 6.084789 + 5.940709 + 5.803949 + 5.791739 + 5.603709 + 5.264269 + 4.946809 + 4.619579 + 4.514579 + 4.433989 + 4.285029 + 4.121419 + 3.945589 + 3.984659 + 4.219099 + 4.546319 + 4.873549 + 5.154379 + 5.388809 + 5.613479 + 5.835699 + 5.884539 + 5.955359 + 5.762439 + 5.459629 + 5.061589 + 4.707499 + 4.458409 + 4.267939 + 4.053039 + 3.943149 + 3.825929 + 3.967569 + 4.280149 + 4.480389 + 4.492599 + 4.390039 + 4.197119 + 4.111649 + 3.982219 + 3.867449 + 3.767319 + 3.872329 + 4.236189 + 4.663539 + 4.971229 + 5.066469 + 4.902859 + 4.675749 + 4.392479 + 4.099439 + 4.114089 + 4.326539 + 4.643999 + 4.971229 + 5.159269 + 5.388809 + 5.576849 + 5.652549 + 5.803949 + 5.913839 + 5.886979 + 5.799069 + 5.730689 + 5.762439 + 5.813719 + 5.821049 + 5.928499 + 6.013969 + 5.764879 + 5.413229 + 5.098219 + 4.678189 + 4.372939 + 4.392479 + 4.590279 + 4.919949 + 5.017629 + 4.858899 + 4.675749 + 4.619579 + 4.834479 + 5.090889 + 5.376599 + 5.681849 + 5.823489 + 5.952919 + 6.062809 + 6.089669 + 6.075019 + 6.026179 + 5.994429 + 6.077459 + 5.857679 + 5.701389 + 5.730689 + 5.784419 + 5.823489 + 5.894309 + 5.762439 + 5.415679 + 4.961459 + 4.595159 + 4.331429 + 4.297239 + 4.582949 + 4.861339 + 5.173919 + 5.166589 + 4.919949 + 4.607369 + 4.370499 + 4.182469 + 4.038389 + 4.145839 + 4.431549 + 4.556089 + 4.480389 + 4.375379 + 4.370499 + 4.558529 + 4.858899 + 4.895529 + 4.741679 + 4.744129 + 4.875989 + 5.105539 + 5.239849 + 5.518239 + 5.652549 + 5.723369 + 5.855239 + 5.962679 + 5.984659 + 5.984659 + 6.055479 + 6.062809 + 6.055479 + 6.070129 + 5.784419 + 5.440099 + 5.056699 + 4.941929 + 5.010299 + 5.134849 + 5.313109 + 5.479169 + 5.623249 + 5.562199 + 5.330209 + 5.010299 + 4.665979 + 4.414459 + 4.201999 + 4.048159 + 4.079899 + 4.189789 + 4.131179 + 4.004199 + 3.916289 + 3.960239 + 4.199559 + 4.624469 + 4.883319 + 5.137289 + 5.379049 + 5.623249 + 5.762439 + 5.833259 + 5.686739 + 5.366839 + 5.225199 + 5.239849 + 5.354629 + 5.508469 + 5.596379 + 5.752669 + 5.874769 + 5.906519 + 5.894309 + 5.742899 + 5.447419 + 5.024959 + 4.883319 + 4.885759 + 4.893089 + 4.714819 + 4.451089 + 4.233749 + 4.043269 + 3.864999 + 3.757559 + 3.669639 + 3.593939 + 3.547539 + 3.506029 + 3.454749 + 3.398579 + 3.361949 + 3.339969 + 3.374159 + 3.520679 + 3.713599 + 3.757559 + 3.779529 + 3.696509 + 3.777089 + 3.886979 + 3.904079 + 3.850349 + 3.965129 + 4.282589 + 4.521899 + 4.714819 + 4.971229 + 5.220319 + 5.532889 + 5.652549 + 5.781979 + 5.955359 + 6.035939 + 6.118969 + 6.133629 + 6.153159 + 6.192229 + 6.143389 + 6.167809 + 5.991989 + 5.652549 + 5.459629 + 5.437649 + 5.339969 + 5.098219 + 4.785639 + 4.492599 + 4.236189 + 4.067689 + 3.933379 + 3.823489 + 3.730689 + 3.611039 + 3.564639 + 3.549989 + 3.557309 + 3.513359 + 3.515799 + 3.694059 + 4.072579 + 4.480389 + 4.705049 + 4.612259 + 4.385149 + 4.201999 + 4.026179 + 3.904079 + 3.774649 + 3.691619 + 3.845469 + 4.201999 + 4.585399 + 4.902859 + 5.256949 + 5.510909 + 5.640339 + 5.843029 + 5.974889 + 5.935819 + 5.821049 + 5.528009 + 5.171479 + 4.810059 + 4.453529 + 4.380269 + 4.565859 + 4.805179 + 5.125079 + 5.354629 + 5.589059 + 5.764879 + 5.923609 + 5.940709 + 5.857679 + 5.694059 + 5.486489 + 5.149499 + 4.844249 + 4.541439 + 4.267939 + 4.060369 + 3.960239 + 3.789299 + 3.642779 + 3.525569 + 3.498699 + 3.454749 + 3.408349 + 3.379049 + 3.376599 + 3.361949 + 3.359509 + 3.369279 + 3.398579 + 3.579289 + 3.948029 + 4.412009 + 4.585399 + 4.514579 + 4.343639 + 4.155599 + 3.984659 + 4.043269 + 4.307009 + 4.421779 + 4.353409 + 4.223979 + 4.053039 + 3.940709 + 3.838139 + 3.730689 + 3.652549 + 3.611039 + 3.564639 + 3.496259 + 3.462069 + 3.454749 + 3.425439 + 3.379049 + 3.432769 + 3.623249 + 3.974889 + 4.380269 + 4.714819 + 5.073799 + 5.369279 + 5.603709 + 5.745349 + 5.652549 + 5.401019 + 5.015189 + 4.709939 + 4.416899 + 4.236189 + 4.236189 + 4.248399 + 4.221539 + 4.297239 + 4.590279 + 4.893089 + 5.134849 + 5.427889 + 5.379049 + 5.364389 + 5.452309 + 5.567079 + 5.672089 + 5.769769 + 5.830819 + 5.923609 + 5.965129 + 6.057919 + 6.050599 + 6.072579 + 6.111649 + 6.070129 + 5.896749 + 5.755109 + 5.718479 + 5.821049 + 6.001759 + 6.001759 + 5.901629 + 5.557309 + 5.173919 + 4.800289 + 4.431549 + 4.194679 + 4.006639 + 3.850349 + 3.747789 + 3.642779 + 3.591499 + 3.569519 + 3.528009 + 3.537779 + 3.554869 + 3.493819 + 3.447419 + 3.440099 + 3.408349 + 3.410789 + 3.452309 + 3.681849 + 4.060369 + 4.441319 + 4.854019 + 5.154379 + 5.425439 + 5.596379 + 5.586619 + 5.354629 + 5.027399 + 4.863779 + 4.761219 + 4.570739 + 4.368059 + 4.397359 + 4.573189 + 4.841809 + 5.203219 + 5.452309 + 5.652549 + 5.855239 + 5.906519 + 5.952919 + 5.828369 + 5.791739 + 5.799069 + 5.813719 + 5.877209 + 5.955359 + 5.781979 + 5.518239 + 5.127519 + 4.763659 + 4.492599 + 4.233749 + 4.011519 + 3.855239 + 3.691619 + 3.635459 + 3.818609 + 4.155599 + 4.590279 + 4.988329 + 5.076239 + 4.907739 + 4.648889 + 4.377829 + 4.216649 + 4.287469 + 4.590279 + 4.846689 + 5.139729 + 5.388809 + 5.689179 + 5.884539 + 6.043269 + 6.170259 + 6.211769 + 6.250839 + 6.209329 + 6.013969 + 5.701389 + 5.469399 + 5.479169 + 5.557309 + 5.728249 + 5.882099 + 5.984659 + 5.901629 + 5.581729 + 5.371719 + 5.418119 + 5.510909 + 5.667199 + 5.791739 + 5.698949 + 5.484049 + 5.154379 + 4.980999 + 5.061589 + 5.195899 + 5.359509 + 5.615919 + 5.762439 + 5.857679 + 5.948029 + 5.835699 + 5.706269 + 5.498699 + 5.188569 + 5.117749 + 5.191009 + 5.315549 + 5.532889 + 5.444979 + 5.396139 + 5.274039 + 5.027399 + 4.744129 + 4.668419 + 4.651329 + 4.514579 + 4.267939 + 4.260609 + 4.263049 + 4.189789 + 4.277699 + 4.600049 + 4.932159 + 5.283809 + 5.528009 + 5.740459 + 5.874769 + 5.955359 + 5.991989 + 5.845469 + 5.528009 + 5.061589 + 4.734359 + 4.534109 + 4.534109 + 4.697729 + 4.744129 + 4.619579 + 4.643999 + 4.832039 + 5.132399 + 5.410789 + 5.625689 + 5.603709 + 5.315549 + 4.961459 + 4.619579 + 4.358289 + 4.155599 + 4.033499 + 3.886979 + 3.772209 + 3.640339 + 3.532889 + 3.435209 + 3.427889 + 3.422999 + 3.398579 + 3.603709 + 4.023729 + 4.451089 + 4.792969 + 4.902859 + 4.780759 + 4.590279 + 4.336309 + 4.145839 + 4.216649 + 4.433989 + 4.714819 + 5.098219 + 5.359509 + 5.569519 + 5.772209 + 5.921169 + 6.055479 + 5.962679 + 5.642779 + 5.435209 + 5.388809 + 5.537779 + 5.681849 + 5.701389 + 5.615919 + 5.667199 + 5.740459 + 5.803949 + 5.882099 + 5.950469 + 6.072579 + 6.148279 + 6.116529 + 6.177579 + 6.201999 + 6.206889 + 5.991989 + 5.564639 + 5.178799 + 4.998089 + 5.051819 + 5.232529 + 5.484049 + 5.686739 + 5.899189 + 5.869889 + 5.977339 + 6.053039 + 6.079899 + 6.128739 + 6.079899 + 6.167809 + 6.194679 + 6.236189 + 6.053039 + 5.652549 + 5.274039 + 4.858899 + 4.534109 + 4.455969 + 4.619579 + 4.866229 + 5.117749 + 5.166589 + 5.056699 + 5.002979 + 5.098219 + 5.325319 + 5.567079 + 5.466959 + 5.252059 + 4.946809 + 4.880879 + 4.980999 + 5.225199 + 5.459629 + 5.723369 + 5.791739 + 5.906519 + 5.991989 + 5.835699 + 5.528009 + 5.142169 + 4.775869 + 4.490159 + 4.236189 + 4.023729 + 3.886979 + 3.752669 + 3.681849 + 3.806399 + 4.145839 + 4.600049 + 5.002979 + 5.303339 + 5.552429 + 5.615919 + 5.523119 + 5.611039 + 5.713599 + 5.845469 + 5.899189 + 5.994429 + 6.092109 + 6.092109 + 6.143389 + 6.153159 + 6.233749 + 6.187349 + 6.013969 + 5.835699 + 5.774649 + 5.686739 + 5.537779 + 5.327759 + 5.054259 + 4.700169 + 4.394919 + 4.180019 + 4.043269 + 3.877209 + 3.752669 + 3.728249 + 3.869889 + 4.206889 + 4.355849 + 4.426669 + 4.453529 + 4.521899 + 4.392479 + 4.155599 + 3.965129 + 3.877209 + 3.970009 + 4.258169 + 4.421779 + 4.336309 + 4.299679 + 4.392479 + 4.675749 + 4.761219 + 4.658659 + 4.490159 + 4.307009 + 4.126299 + 3.972449 + 4.077459 + 4.372939 + 4.741679 + 5.088449 + 5.186129 + 5.037169 + 4.785639 + 4.563419 + 4.534109 + 4.705049 + 4.741679 + 4.648889 + 4.431549 + 4.238629 + 4.065249 + 3.943149 + 3.811279 + 3.691619 + 3.652549 + 3.825929 + 4.223979 + 4.424219 + 4.429109 + 4.319219 + 4.138509 + 3.965129 + 3.886979 + 3.801509 + 3.701389 + 3.640339 + 3.767319 + 4.150719 + 4.648889 + 4.990769 + 5.088449 + 5.022509 + 4.783199 + 4.685519 + 4.665979 + 4.707499 + 4.912619 + 5.195899 + 5.415679 + 5.623249 + 5.740459 + 5.899189 + 5.928499 + 6.050599 + 6.153159 + 5.965129 + 5.586619 + 5.381489 + 5.371719 + 5.486489 + 5.567079 + 5.821049 + 5.913839 + 5.994429 + 6.011519 + 5.999309 + 6.018849 + 5.821049 + 5.728249 + 5.740459 + 5.764879 + 5.882099 + 5.926049 + 5.750229 + 5.415679 + 4.995649 + 4.861339 + 4.902859 + 5.103099 + 5.364389 + 5.596379 + 5.752669 + 5.845469 + 5.928499 + 6.006639 + 5.840579 + 5.518239 + 5.173919 + 4.739239 + 4.458409 + 4.426669 + 4.602489 + 4.822269 + 5.183689 + 5.430329 + 5.652549 + 5.821049 + 5.706269 + 5.369279 + 5.027399 + 4.705049 + 4.414459 + 4.145839 + 3.965129 + 4.033499 + 4.372939 + 4.683079 +]; + +dat = iddata (Y, U) + +% IB01AD EXAMPLE PROGRAM DATA +% 15 0 1 1 1000 0.0 -1.0 M C N O N N +% nobr, n, m, l, nsmp, rcond, tol, meth, alg, jobd, batch,conct,ctrl + +% IB01BD EXAMPLE PROGRAM DATA +% 15 0 1 1 1000 0.0 -1.0 C C N O N N A K + +nobr = 15; +meth = 2; +alg = 0; +jobd = 1; +batch = 3; +conct = 1; +ctrl = 1; +rcond = 0.0; +tol = -1.0; + +[a, b, c, d, q, ry, s, k] = slident (Y, U, nobr, meth, alg, jobd, batch, conct, ctrl, rcond, tol); + +n +sv + +%{ + IB01AD EXAMPLE PROGRAM RESULTS + + The order of the system is 4 + The singular values are + 69.8841 14.9963 3.6675 1.9677 0.3000 0.2078 0.1651 0.1373 + 0.1133 0.1059 0.0856 0.0784 0.0733 0.0678 0.0571 + + + IB01BD EXAMPLE PROGRAM RESULTS + + + The system state matrix A is + 0.8924 0.3887 0.1285 0.1716 + -0.0837 0.6186 -0.6273 -0.4582 + 0.0052 0.1307 0.6685 -0.6755 + 0.0055 0.0734 -0.2148 0.4788 + + The system output matrix C is + -0.4442 0.6663 0.3961 0.4102 + + The system input matrix B is + -0.2142 + -0.1968 + 0.0525 + 0.0361 + + The system input-output matrix D is + -0.0041 + + The Kalman gain matrix K is + -1.9513 + -0.1867 + 0.6348 + -0.3486 + + The state covariance matrix Q is + 0.0052 0.0005 -0.0017 0.0009 + 0.0005 0.0000 -0.0002 0.0001 + -0.0017 -0.0002 0.0006 -0.0003 + 0.0009 0.0001 -0.0003 0.0002 + + The output covariance matrix Ry is + 0.0012 + + The state-output cross-covariance matrix S is + -0.0025 + -0.0002 + 0.0008 + -0.0005 + +%} Modified: trunk/octave-forge/extra/control-devel/src/devel_slicot_functions.cc =================================================================== --- trunk/octave-forge/extra/control-devel/src/devel_slicot_functions.cc 2012-03-16 15:03:53 UTC (rev 9921) +++ trunk/octave-forge/extra/control-devel/src/devel_slicot_functions.cc 2012-03-16 18:08:52 UTC (rev 9922) @@ -1,2 +1,2 @@ -#include "slib01ad.cc" // preprocess the input-output data +// #include "slib01ad.cc" // preprocess the input-output data #include "slident.cc" // system identification Modified: trunk/octave-forge/extra/control-devel/src/slident.cc =================================================================== --- trunk/octave-forge/extra/control-devel/src/slident.cc 2012-03-16 15:03:53 UTC (rev 9921) +++ trunk/octave-forge/extra/control-devel/src/slident.cc 2012-03-16 18:08:52 UTC (rev 9922) @@ -79,7 +79,7 @@ int nargin = args.length (); octave_value_list retval; - if (nargin != 12) + if (nargin != 11) { print_usage (); } @@ -92,6 +92,8 @@ char batch; char conct; char ctrl; + char metha; + char jobda; Matrix y = args(0).matrix_value (); Matrix u = args(1).matrix_value (); @@ -106,13 +108,23 @@ double rcond = args(9).double_value (); double tol = args(10).double_value (); - double tolb = args(11).double_value (); - + double tolb = args(9).double_value (); // tolb = rcond - if (imeth == 0) - meth = 'M'; - else - meth = 'N'; + + switch (imeth) + { + case 0: + meth = 'M'; + metha = 'M'; + case 1: + meth = 'N'; + metha = 'N'; + case 3: + meth = 'C'; + metha = 'N'; // no typo here + default: + error ("slib01ad: argument 'meth' invalid"); + } switch (ialg) { @@ -129,7 +141,9 @@ error ("slib01ad: argument 'alg' invalid"); } - if (ijobd == 0) + if (meth == 'C') + jobd = 'N'; + else if (ijobd == 0) jobd = 'M'; else jobd = 'N'; @@ -182,9 +196,9 @@ int n; int ldr; - if (meth == 'M' && jobd == 'M') + if (metha == 'M' && jobd == 'M') ldr = max (2*(m+l)*nobr, 3*m*nobr); - else if (meth == 'N' || (meth == 'M' && jobd == 'N')) + else if (metha == 'N' || (metha == 'M' && jobd == 'N')) ldr = 2*(m+l)*nobr; else error ("slib01ad: could not handle 'ldr' case"); @@ -195,7 +209,7 @@ // workspace int liwork; - if (meth == 'N') // if METH = 'N' + if (metha == 'N') // if METH = 'N' liwork = (m+l)*nobr; else if (alg == 'F') // if METH = 'M' and ALG = 'F' liwork = m+l; @@ -205,8 +219,6 @@ // TODO: Handle 'k' for DWORK int ldwork; - - ldwork = 0; if (alg == 'C') { @@ -217,7 +229,7 @@ else // (conct == 'N') ldwork = 1; } - else if (meth == 'M') // && (batch == 'L' || batch == 'O') + else if (metha == 'M') // && (batch == 'L' || batch == 'O') { if (conct == 'C' && batch == 'L') ldwork = max ((4*nobr-2)*(m+l), 5*l*nobr); @@ -233,44 +245,11 @@ } else if (alg == 'F') { -/* -For the second LDWORK case, code and documentation don't match: -doc line 276: BATCH = 'F', 'I' -code line 586: BATCH = 'F', 'I', 'O' -The third case with BATCH = 'O' is never reached. - - -IB01AD.f Lines 273-279: -C LDWORK >= (M+L)*2*NOBR*(M+L+3), if ALG = 'F', -C BATCH <> 'O' and CONCT = 'C'; -C LDWORK >= (M+L)*2*NOBR*(M+L+1), if ALG = 'F', -C BATCH = 'F', 'I' and CONCT = 'N'; -C LDWORK >= (M+L)*4*NOBR*(M+L+1)+(M+L)*2*NOBR, if ALG = 'F', -C BATCH = 'L' and CONCT = 'N', or -C BATCH = 'O'; - - -IB01AD.f Lines 499-500: - ONEBCH = LSAME( BATCH, 'O' ) - FIRST = LSAME( BATCH, 'F' ) .OR. ONEBCH - - -IB01AD.f Lines 583-591: - ELSE IF ( FQRALG ) THEN - IF ( .NOT.ONEBCH .AND. CONNEC ) THEN - MINWRK = NR*( M + L + 3 ) - ELSE IF ( FIRST .OR. INTERM ) THEN // (batch = F || O) || batch = I - MINWRK = NR*( M + L + 1 ) ^ - ELSE | - MINWRK = 2*NR*( M + L + 1 ) + NR ??? - END IF - ELSE -*/ if (batch != 'O' && conct == 'C') ldwork = (m+l)*2*nobr*(m+l+3); - else if (batch == 'F' || batch == 'O' || batch == 'I') // && conct == 'N' + else if (batch == 'F' || batch == 'I') // && conct == 'N' ldwork = (m+l)*2*nobr*(m+l+1); - else // (batch == 'L' && conct == 'N') + else // (batch == 'L' || '0' && conct == 'N') ldwork = (m+l)*4*nobr*(m+l+1)+(m+l)*2*nobr; } else // (alg == 'Q') @@ -283,7 +262,7 @@ } else if (ldr >= ns && batch == 'O') { - if (meth == 'M') + if (metha == 'M') ldwork = max (4*(m+l)*nobr, 5*l*nobr); else // (meth == 'N') ldwork = 5*(m+l)*nobr + 1; @@ -321,7 +300,7 @@ // SLICOT routine IB01AD F77_XFCN (ib01ad, IB01AD, - (meth, alg, jobd, + (metha, alg, jobd, batch, conct, ctrl, nobr, m, l, nsmp, This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2012-03-16 18:23:57
|
Revision: 9923 http://octave.svn.sourceforge.net/octave/?rev=9923&view=rev Author: paramaniac Date: 2012-03-16 18:23:51 +0000 (Fri, 16 Mar 2012) Log Message: ----------- control-devel: fix silly mistakes, identification works and computes slicot test case correctly Modified Paths: -------------- trunk/octave-forge/extra/control-devel/devel/test_slident.m trunk/octave-forge/extra/control-devel/src/slident.cc Modified: trunk/octave-forge/extra/control-devel/devel/test_slident.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/test_slident.m 2012-03-16 18:08:52 UTC (rev 9922) +++ trunk/octave-forge/extra/control-devel/devel/test_slident.m 2012-03-16 18:23:51 UTC (rev 9923) @@ -2031,10 +2031,10 @@ rcond = 0.0; tol = -1.0; -[a, b, c, d, q, ry, s, k] = slident (Y, U, nobr, meth, alg, jobd, batch, conct, ctrl, rcond, tol); +[a, b, c, d, q, ry, s, k] = slident (Y, U, nobr, meth, alg, jobd, batch, conct, ctrl, rcond, tol) -n -sv +%n +%sv %{ IB01AD EXAMPLE PROGRAM RESULTS Modified: trunk/octave-forge/extra/control-devel/src/slident.cc =================================================================== --- trunk/octave-forge/extra/control-devel/src/slident.cc 2012-03-16 18:08:52 UTC (rev 9922) +++ trunk/octave-forge/extra/control-devel/src/slident.cc 2012-03-16 18:23:51 UTC (rev 9923) @@ -116,12 +116,15 @@ case 0: meth = 'M'; metha = 'M'; + break; case 1: meth = 'N'; metha = 'N'; - case 3: + break; + case 2: meth = 'C'; metha = 'N'; // no typo here + break; default: error ("slib01ad: argument 'meth' invalid"); } This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2012-03-17 08:50:54
|
Revision: 9931 http://octave.svn.sourceforge.net/octave/?rev=9931&view=rev Author: paramaniac Date: 2012-03-17 08:50:48 +0000 (Sat, 17 Mar 2012) Log Message: ----------- control-devel: integrate Slicot IB01CD Modified Paths: -------------- trunk/octave-forge/extra/control-devel/src/slident.cc Removed Paths: ------------- trunk/octave-forge/extra/control-devel/devel/makefile_ident.m Deleted: trunk/octave-forge/extra/control-devel/devel/makefile_ident.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/makefile_ident.m 2012-03-17 07:20:50 UTC (rev 9930) +++ trunk/octave-forge/extra/control-devel/devel/makefile_ident.m 2012-03-17 08:50:48 UTC (rev 9931) @@ -1,54 +0,0 @@ -## ============================================================================== -## Developer Makefile for OCT-files -## ============================================================================== -## USAGE: * fetch control-devel from Octave-Forge by svn -## * add control-devel/inst, control-devel/src and control-devel/devel -## to your Octave path (by an .octaverc file) -## * run makefile_ident -## ============================================================================== - -homedir = pwd (); -develdir = fileparts (which ("makefile_ident")); -srcdir = [develdir, "/../src"]; -cd (srcdir); - -## preprocess the input-output data -mkoctfile IB01AD.f IB01MD.f IB01ND.f IB01OD.f IB01MY.f \ - MB04OD.f MB03UD.f MB04ID.f MA02AD.f MB03OD.f \ - MB04IY.f IB01OY.f MA02ED.f MA02FD.f MB04OY.f \ - "$(mkoctfile -p LAPACK_LIBS)" \ - "$(mkoctfile -p BLAS_LIBS)" \ - "$(mkoctfile -p FLIBS)" - -## estimating system matrices, Kalman gain, and covariances -mkoctfile IB01BD.f IB01PD.f MA02AD.f SB02MT.f SB02RD.f \ - SB02ND.f MB02UD.f MA02ED.f IB01PY.f MB03OD.f \ - MB02QY.f IB01PX.f SB02MS.f SB02RU.f SB02SD.f \ - MB01RU.f SB02QD.f SB02MV.f SB02MW.f SB02MR.f \ - MB02PD.f MB01SD.f MB04KD.f MB03UD.f MB04OD.f \ - MB04OY.f MB01VD.f select.f MB01UD.f SB03SY.f \ - MB01RX.f SB03MX.f SB03SX.f MB01RY.f SB03QY.f \ - SB03QX.f SB03MY.f SB04PX.f SB03MV.f SB03MW.f \ - "$(mkoctfile -p LAPACK_LIBS)" \ - "$(mkoctfile -p BLAS_LIBS)" - -## estimating the initial state -mkoctfile IB01CD.f TB01WD.f IB01RD.f IB01QD.f select.f \ - MB01TD.f MA02AD.f MB04OD.f MB04OY.f MB02UD.f \ - MB03UD.f MB01SD.f \ - "$(mkoctfile -p LAPACK_LIBS)" \ - "$(mkoctfile -p BLAS_LIBS)" \ - "$(mkoctfile -p FLIBS)" - -## fit state-space model to frequency response data -mkoctfile slsb10yd.cc \ - SB10YD.f DG01MD.f AB04MD.f SB10ZP.f AB07ND.f \ - MC01PD.f TD04AD.f TD03AY.f TB01PD.f TB01XD.f \ - AB07MD.f TB01UD.f TB01ID.f MB01PD.f MB03OY.f \ - MB01QD.f \ - "$(mkoctfile -p LAPACK_LIBS)" \ - "$(mkoctfile -p BLAS_LIBS)" - -system ("rm *.o"); -cd (homedir); - Modified: trunk/octave-forge/extra/control-devel/src/slident.cc =================================================================== --- trunk/octave-forge/extra/control-devel/src/slident.cc 2012-03-17 07:20:50 UTC (rev 9930) +++ trunk/octave-forge/extra/control-devel/src/slident.cc 2012-03-17 08:50:48 UTC (rev 9931) @@ -556,7 +556,75 @@ //////////////////////////////////////////////////////////////////////////////////// // SLICOT IB01CD - estimating the initial state // //////////////////////////////////////////////////////////////////////////////////// +/* + // arguments in + char jobx0 = 'X'; + char comuse = 'U'; + char jobbd = 'D'; + + // arguments out + int ldv = max (1, n); + + ColumnVector x0 (n); + Matrix (ldv, n); + + // workspace + int liwork_c = n; // if JOBX0 = 'X' and COMUSE <> 'C' + int ldwork_c; + + ldwork_c = ldw1 + n*( n + m + l ) + max (5*n, ldw1, min (ldw2, ldw3)) + + + // I can't find a definition for parameter 't' + +C The length of the array DWORK. +C LDWORK >= 2, if JOBX0 = 'N' and COMUSE <> 'C', or +C if max( N, M ) = 0. +C Otherwise, +C LDWORK >= LDW1 + N*( N + M + L ) + +C max( 5*N, LDW1, min( LDW2, LDW3 ) ), +C where, if COMUSE = 'C', then +C LDW1 = 2, if M = 0 or JOB = 'B', +C LDW1 = 3, if M > 0 and JOB = 'D', +C LDWa = t*L*(r + 1) + max( N + max( d, f ), 6*r ), +C LDW2 = LDWa, if M = 0 or JOB = 'B', +C LDW2 = max( LDWa, t*L*(r + 1) + 2*M*M + 6*M ), +C if M > 0 and JOB = 'D', +C LDWb = (b + r)*(r + 1) + +C max( q*(r + 1) + N*N*M + c + max( d, f ), 6*r ), +C LDW3 = LDWb, if M = 0 or JOB = 'B', +C LDW3 = max( LDWb, (b + r)*(r + 1) + 2*M*M + 6*M ), +C if M > 0 and JOB = 'D', +C r = N*M + a, +C a = 0, if JOBX0 = 'N', +C a = N, if JOBX0 = 'X'; +C b = 0, if JOB = 'B', +C b = L*M, if JOB = 'D'; +C c = 0, if JOBX0 = 'N', +C c = L*N, if JOBX0 = 'X'; +C d = 0, if JOBX0 = 'N', +C d = 2*N*N + N, if JOBX0 = 'X'; +C f = 2*r, if JOB = 'B' or M = 0, +C f = M + max( 2*r, M ), if JOB = 'D' and M > 0; +C q = b + r*L; +C and, if JOBX0 = 'X' and COMUSE <> 'C', then +C LDW1 = 2, +C LDW2 = t*L*(N + 1) + 2*N + max( 2*N*N, 4*N ), +C LDW3 = N*(N + 1) + 2*N + max( q*(N + 1) + 2*N*N + L*N, +C 4*N ), +C q = N*L. +C For good performance, LDWORK should be larger. + + + OCTAVE_LOCAL_BUFFER (int, iwork_c, liwork_c); + OCTAVE_LOCAL_BUFFER (double, dwork_c, ldwork_c); + + // error indicators + int iwarn_c = 0; + int info_c = 0; + + // SLICOT routine IB01CD F77_XFCN (ib01cd, IB01CD, (jobx0, comuse, job***, @@ -603,6 +671,7 @@ error_msg ("ident", info_c, 2, err_msg_c); warning_msg ("ident", iwarn_c, 6, warn_msg_c); +*/ // return values retval(0) = a; @@ -614,6 +683,8 @@ retval(5) = ry; retval(6) = s; retval(7) = k; + + //retval(8) = x0; //retval(0) = octave_value (n); //retval(1) = r; //retval(2) = sv; This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2012-03-17 16:00:11
|
Revision: 9936 http://octave.svn.sourceforge.net/octave/?rev=9936&view=rev Author: paramaniac Date: 2012-03-17 16:00:05 +0000 (Sat, 17 Mar 2012) Log Message: ----------- control-devel: check results with assert Modified Paths: -------------- trunk/octave-forge/extra/control-devel/devel/test_slident.m trunk/octave-forge/extra/control-devel/src/slident.cc Modified: trunk/octave-forge/extra/control-devel/devel/test_slident.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/test_slident.m 2012-03-17 12:41:13 UTC (rev 9935) +++ trunk/octave-forge/extra/control-devel/devel/test_slident.m 2012-03-17 16:00:05 UTC (rev 9936) @@ -2033,6 +2033,48 @@ [a, b, c, d, q, ry, s, k] = slident (Y, U, nobr, meth, alg, jobd, batch, conct, ctrl, rcond, tol) + +ae = [ 0.8924 0.3887 0.1285 0.1716 + -0.0837 0.6186 -0.6273 -0.4582 + 0.0052 0.1307 0.6685 -0.6755 + 0.0055 0.0734 -0.2148 0.4788 ]; + +ce = [ -0.4442 0.6663 0.3961 0.4102 ]; + +be = [ -0.2142 + -0.1968 + 0.0525 + 0.0361 ]; + +de = [ -0.0041 ]; + +ke = [ -1.9513 + -0.1867 + 0.6348 + -0.3486 ]; + +qe = [ 0.0052 0.0005 -0.0017 0.0009 + 0.0005 0.0000 -0.0002 0.0001 + -0.0017 -0.0002 0.0006 -0.0003 + 0.0009 0.0001 -0.0003 0.0002 ]; + +rye = [ 0.0012 ]; + +se = [ -0.0025 + -0.0002 + 0.0008 + -0.0005 ]; + +assert (a, ae, 1e-4); +assert (b, be, 1e-4); +assert (c, ce, 1e-4); +assert (d, de, 1e-4); +assert (k, ke, 1e-4); +assert (q, qe, 1e-4); +assert (ry, rye, 1e-4); +assert (s, se, 1e-4); + + %n %sv Modified: trunk/octave-forge/extra/control-devel/src/slident.cc =================================================================== --- trunk/octave-forge/extra/control-devel/src/slident.cc 2012-03-17 12:41:13 UTC (rev 9935) +++ trunk/octave-forge/extra/control-devel/src/slident.cc 2012-03-17 16:00:05 UTC (rev 9936) @@ -556,6 +556,11 @@ //////////////////////////////////////////////////////////////////////////////////// // SLICOT IB01CD - estimating the initial state // //////////////////////////////////////////////////////////////////////////////////// + +// TODO: use only one iwork and dwork for all three slicot routines +// ldwork = max (ldwork_a, ldwork_b, ldwork_c) + + /* // arguments in char jobx0 = 'X'; This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2012-03-19 16:16:58
|
Revision: 9970 http://octave.svn.sourceforge.net/octave/?rev=9970&view=rev Author: paramaniac Date: 2012-03-19 16:16:49 +0000 (Mon, 19 Mar 2012) Log Message: ----------- control-devel: calculate initial state vector Modified Paths: -------------- trunk/octave-forge/extra/control-devel/devel/test_slident.m trunk/octave-forge/extra/control-devel/src/slident.cc Modified: trunk/octave-forge/extra/control-devel/devel/test_slident.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/test_slident.m 2012-03-19 15:35:53 UTC (rev 9969) +++ trunk/octave-forge/extra/control-devel/devel/test_slident.m 2012-03-19 16:16:49 UTC (rev 9970) @@ -2031,7 +2031,7 @@ rcond = 0.0; tol = -1.0; -[a, b, c, d, q, ry, s, k] = slident (Y, U, nobr, meth, alg, jobd, batch, conct, ctrl, rcond, tol) +[a, b, c, d, q, ry, s, k, x0] = slident (Y, U, nobr, meth, alg, jobd, batch, conct, ctrl, rcond, tol) ae = [ 0.8924 0.3887 0.1285 0.1716 @@ -2077,23 +2077,27 @@ -figure (1) -plot (Y) +% figure (1) +% plot (Y) P = ss (a, b, c, d, 1); -figure (2) -lsim (P, U) -%lsim (P, U, [], x0) % initial values dependent on realization, IB01BD != IB01CD +% figure (2) +% lsim (P, U) +% lsim (P, U, [], x0) % initial values dependent on realization, IB01BD != IB01CD -[y, t] = lsim (P, U); +% [y, t] = lsim (P, U); +[y, t] = lsim (P, U, [], x0); + figure (3) plot (t, Y, 'b', t, y, 'r') legend ('y measured', 'y simulated', 'location', 'southeast') -axis tight +%axis tight +% print -depsc2 ident.eps + %n %sv Modified: trunk/octave-forge/extra/control-devel/src/slident.cc =================================================================== --- trunk/octave-forge/extra/control-devel/src/slident.cc 2012-03-19 15:35:53 UTC (rev 9969) +++ trunk/octave-forge/extra/control-devel/src/slident.cc 2012-03-19 16:16:49 UTC (rev 9970) @@ -561,7 +561,6 @@ // ldwork = max (ldwork_a, ldwork_b, ldwork_c) -/* // arguments in char jobx0 = 'X'; char comuse = 'U'; @@ -571,57 +570,19 @@ int ldv = max (1, n); ColumnVector x0 (n); - Matrix (ldv, n); + Matrix v (ldv, n); // workspace int liwork_c = n; // if JOBX0 = 'X' and COMUSE <> 'C' int ldwork_c; - - ldwork_c = ldw1 + n*( n + m + l ) + max (5*n, ldw1, min (ldw2, ldw3)) - - - // I can't find a definition for parameter 't' - -C The length of the array DWORK. -C LDWORK >= 2, if JOBX0 = 'N' and COMUSE <> 'C', or -C if max( N, M ) = 0. -C Otherwise, -C LDWORK >= LDW1 + N*( N + M + L ) + -C max( 5*N, LDW1, min( LDW2, LDW3 ) ), -C where, if COMUSE = 'C', then -C LDW1 = 2, if M = 0 or JOB = 'B', -C LDW1 = 3, if M > 0 and JOB = 'D', -C LDWa = t*L*(r + 1) + max( N + max( d, f ), 6*r ), -C LDW2 = LDWa, if M = 0 or JOB = 'B', -C LDW2 = max( LDWa, t*L*(r + 1) + 2*M*M + 6*M ), -C if M > 0 and JOB = 'D', -C LDWb = (b + r)*(r + 1) + -C max( q*(r + 1) + N*N*M + c + max( d, f ), 6*r ), -C LDW3 = LDWb, if M = 0 or JOB = 'B', -C LDW3 = max( LDWb, (b + r)*(r + 1) + 2*M*M + 6*M ), -C if M > 0 and JOB = 'D', -C r = N*M + a, -C a = 0, if JOBX0 = 'N', -C a = N, if JOBX0 = 'X'; -C b = 0, if JOB = 'B', -C b = L*M, if JOB = 'D'; -C c = 0, if JOBX0 = 'N', -C c = L*N, if JOBX0 = 'X'; -C d = 0, if JOBX0 = 'N', -C d = 2*N*N + N, if JOBX0 = 'X'; -C f = 2*r, if JOB = 'B' or M = 0, -C f = M + max( 2*r, M ), if JOB = 'D' and M > 0; -C q = b + r*L; -C and, if JOBX0 = 'X' and COMUSE <> 'C', then -C LDW1 = 2, -C LDW2 = t*L*(N + 1) + 2*N + max( 2*N*N, 4*N ), -C LDW3 = N*(N + 1) + 2*N + max( q*(N + 1) + 2*N*N + L*N, -C 4*N ), -C q = N*L. -C For good performance, LDWORK should be larger. + int t = nsmp; + + int ldw1_c = 2; + int ldw2_c = t*l*(n + 1) + 2*n + max (2*n*n, 4*n); + int ldw3_c = n*(n + 1) + 2*n + max (n*l*(n + 1) + 2*n*n + l*n, 4*n); + ldwork_c = ldw1_c + n*( n + m + l ) + max (5*n, ldw1_c, min (ldw2_c, ldw3_c)); - OCTAVE_LOCAL_BUFFER (int, iwork_c, liwork_c); OCTAVE_LOCAL_BUFFER (double, dwork_c, ldwork_c); @@ -632,7 +593,7 @@ // SLICOT routine IB01CD F77_XFCN (ib01cd, IB01CD, - (jobx0, comuse, job***, + (jobx0, comuse, jobbd, n, m, l, nsmp, a.fortran_vec (), lda, @@ -676,7 +637,7 @@ error_msg ("ident", info_c, 2, err_msg_c); warning_msg ("ident", iwarn_c, 6, warn_msg_c); -*/ + // return values retval(0) = a; @@ -689,7 +650,7 @@ retval(6) = s; retval(7) = k; - //retval(8) = x0; + retval(8) = x0; //retval(0) = octave_value (n); //retval(1) = r; //retval(2) = sv; This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2012-03-20 10:36:56
|
Revision: 9976 http://octave.svn.sourceforge.net/octave/?rev=9976&view=rev Author: paramaniac Date: 2012-03-20 10:36:50 +0000 (Tue, 20 Mar 2012) Log Message: ----------- control-devel: first steps to fix crash with DaISy data Modified Paths: -------------- trunk/octave-forge/extra/control-devel/devel/destillation.m trunk/octave-forge/extra/control-devel/src/slident.cc Modified: trunk/octave-forge/extra/control-devel/devel/destillation.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/destillation.m 2012-03-20 09:37:44 UTC (rev 9975) +++ trunk/octave-forge/extra/control-devel/devel/destillation.m 2012-03-20 10:36:50 UTC (rev 9976) @@ -63,7 +63,16 @@ Y_dest_n30=Y(:,10:12); -dat = iddata (Y_dest, Y_dest) +dat = iddata (Y_dest, U_dest) -[sys, x0] = ident (dat, 15) % nobr? 5, 90, ? +[sys, x0] = ident (dat, 5) % nobr? 5, 90, ? + +[y, t] = lsim (sys, U_dest, [], x0); + +figure (1) +plot (t, Y_dest, 'b') +%plot (t, Y_dest, 'b', t, y, 'r') +legend ('y measured', 'y simulated', 'location', 'southeast') + + Modified: trunk/octave-forge/extra/control-devel/src/slident.cc =================================================================== --- trunk/octave-forge/extra/control-devel/src/slident.cc 2012-03-20 09:37:44 UTC (rev 9975) +++ trunk/octave-forge/extra/control-devel/src/slident.cc 2012-03-20 10:36:50 UTC (rev 9976) @@ -206,6 +206,17 @@ // y.rows == u.rows is checked by iddata class // TODO: check minimal nsmp size + if (batch == 'O') + { + if (nsmp < 2*(m+l+1)*nobr - 1) + error ("slident: require NSMP >= 2*(M+L+1)*NOBR - 1"); + } + else + { + if (nsmp < 2*nobr) + error ("slident: require NSMP >= 2*NOBR"); + } + int ldu; if (m == 0) @@ -299,6 +310,11 @@ ldwork = 6*(m+l)*nobr; } } +//////////////////////////////////////////////////////////////////// +// TO BE REMOVED !!! +//////////////////////////////////////////////////////////////////// +ldwork = 100000; +//////////////////////////////////////////////////////////////////// /* IB01AD.f Lines 291-195: This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2012-03-20 13:11:44
|
Revision: 9977 http://octave.svn.sourceforge.net/octave/?rev=9977&view=rev Author: paramaniac Date: 2012-03-20 13:11:38 +0000 (Tue, 20 Mar 2012) Log Message: ----------- control-devel: first steps to fix crash with DaISy data (2) Modified Paths: -------------- trunk/octave-forge/extra/control-devel/devel/destillation.m trunk/octave-forge/extra/control-devel/devel/ident.m trunk/octave-forge/extra/control-devel/src/slident.cc Modified: trunk/octave-forge/extra/control-devel/devel/destillation.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/destillation.m 2012-03-20 10:36:50 UTC (rev 9976) +++ trunk/octave-forge/extra/control-devel/devel/destillation.m 2012-03-20 13:11:38 UTC (rev 9977) @@ -65,10 +65,11 @@ dat = iddata (Y_dest, U_dest) -[sys, x0] = ident (dat, 5) % nobr? 5, 90, ? +[sys, x0] = ident (dat) %, 5) % nobr? 5, 90, ? [y, t] = lsim (sys, U_dest, [], x0); +%[y, t] = lsim (sys, U_dest); figure (1) plot (t, Y_dest, 'b') Modified: trunk/octave-forge/extra/control-devel/devel/ident.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/ident.m 2012-03-20 10:36:50 UTC (rev 9976) +++ trunk/octave-forge/extra/control-devel/devel/ident.m 2012-03-20 13:11:38 UTC (rev 9977) @@ -6,11 +6,18 @@ jobd = 1; batch = 3; conct = 1; - ctrl = 1; + ctrl = 0; %1; rcond = 0.0; tol = -1.0; + + [n, l, m, e] = size (dat); + + nsmp = n(1) + nobr = fix ((nsmp+1)/(2*(m+l+1))) + % nsmp >= 2*(m+l+1)*nobr - 1 + % nobr <= (nsmp+1)/(2*(m+l+1)) - [a, b, c, d, q, ry, s, k, x0] = slident (dat.y{1}, dat.u{1}, nobr, meth, alg, jobd, batch, conct, ctrl, rcond, tol) + [a, b, c, d, q, ry, s, k, x0] = slident (dat.y{1}, dat.u{1}, nobr, meth, alg, jobd, batch, conct, ctrl, rcond, tol); sys = ss (a, b, c, d, -1); Modified: trunk/octave-forge/extra/control-devel/src/slident.cc =================================================================== --- trunk/octave-forge/extra/control-devel/src/slident.cc 2012-03-20 10:36:50 UTC (rev 9976) +++ trunk/octave-forge/extra/control-devel/src/slident.cc 2012-03-20 13:11:38 UTC (rev 9977) @@ -313,7 +313,7 @@ //////////////////////////////////////////////////////////////////// // TO BE REMOVED !!! //////////////////////////////////////////////////////////////////// -ldwork = 100000; +ldwork = 1000000; //////////////////////////////////////////////////////////////////// /* This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2012-03-21 10:21:37
|
Revision: 9989 http://octave.svn.sourceforge.net/octave/?rev=9989&view=rev Author: paramaniac Date: 2012-03-21 10:21:26 +0000 (Wed, 21 Mar 2012) Log Message: ----------- control-devel: first steps towards an interface for ident function (2) Modified Paths: -------------- trunk/octave-forge/extra/control-devel/devel/destillation.m trunk/octave-forge/extra/control-devel/devel/glass_furnace.m trunk/octave-forge/extra/control-devel/devel/ident.m trunk/octave-forge/extra/control-devel/devel/test_slident.m trunk/octave-forge/extra/control-devel/src/slident.cc Modified: trunk/octave-forge/extra/control-devel/devel/destillation.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/destillation.m 2012-03-21 09:48:12 UTC (rev 9988) +++ trunk/octave-forge/extra/control-devel/devel/destillation.m 2012-03-21 10:21:26 UTC (rev 9989) @@ -65,7 +65,7 @@ dat = iddata (Y_dest, U_dest) -[sys, x0] = ident (dat, 5) % s=5, n=4 +[sys, x0] = ident (dat, 5, 4) % s=5, n=4 [y, t] = lsim (sys, U_dest, [], x0); Modified: trunk/octave-forge/extra/control-devel/devel/glass_furnace.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/glass_furnace.m 2012-03-21 09:48:12 UTC (rev 9988) +++ trunk/octave-forge/extra/control-devel/devel/glass_furnace.m 2012-03-21 10:21:26 UTC (rev 9989) @@ -47,7 +47,7 @@ dat = iddata (Y, U) -[sys, x0] = ident (dat, 10) % s=10, n=5 +[sys, x0] = ident (dat, 10, 5) % s=10, n=5 [y, t] = lsim (sys, U, [], x0); Modified: trunk/octave-forge/extra/control-devel/devel/ident.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/ident.m 2012-03-21 09:48:12 UTC (rev 9988) +++ trunk/octave-forge/extra/control-devel/devel/ident.m 2012-03-21 10:21:26 UTC (rev 9989) @@ -18,6 +18,7 @@ nsmp = ns(1); nobr = fix ((nsmp+1)/(2*(m+l+1))); ctrl = 0; # confirm system order estimate + n = 0; % nsmp >= 2*(m+l+1)*nobr - 1 % nobr <= (nsmp+1)/(2*(m+l+1)) elseif (isempty (s)) @@ -25,11 +26,11 @@ nsmp = ns(1); nobr = fix ((nsmp+1)/(2*(m+l+1))); nobr = min (nobr, s); - %% ctrl = 1; # no confirmation - ctrl = 0; % setting of n not yet possible + ctrl = 1; # no confirmation elseif (isempty (n)) nobr = s; ctrl = 0; # confirm system order estimate + n = 0; else # s & n non-empty nsmp = ns(1); nobr = fix ((nsmp+1)/(2*(m+l+1))); @@ -37,6 +38,7 @@ error ("ident: s > nobr"); endif nobr = s; + ctrl = 1; ## TODO: specify n for IB01BD endif @@ -45,7 +47,7 @@ % nsmp >= 2*(m+l+1)*nobr - 1 % nobr <= (nsmp+1)/(2*(m+l+1)) %nobr = 10 - [a, b, c, d, q, ry, s, k, x0] = slident (dat.y{1}, dat.u{1}, nobr, meth, alg, jobd, batch, conct, ctrl, rcond, tol); + [a, b, c, d, q, ry, s, k, x0] = slident (dat.y{1}, dat.u{1}, nobr, n, meth, alg, jobd, batch, conct, ctrl, rcond, tol); sys = ss (a, b, c, d, -1); Modified: trunk/octave-forge/extra/control-devel/devel/test_slident.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/test_slident.m 2012-03-21 09:48:12 UTC (rev 9988) +++ trunk/octave-forge/extra/control-devel/devel/test_slident.m 2012-03-21 10:21:26 UTC (rev 9989) @@ -2030,8 +2030,9 @@ ctrl = 1; rcond = 0.0; tol = -1.0; +n = 0; -[a, b, c, d, q, ry, s, k, x0] = slident (Y, U, nobr, meth, alg, jobd, batch, conct, ctrl, rcond, tol) +[a, b, c, d, q, ry, s, k, x0] = slident (Y, U, nobr, n, meth, alg, jobd, batch, conct, ctrl, rcond, tol) ae = [ 0.8924 0.3887 0.1285 0.1716 @@ -2090,6 +2091,7 @@ % [y, t] = lsim (P, U); [y, t] = lsim (P, U, [], x0); +err = norm (Y - y, 1) / norm (Y, 1) figure (3) plot (t, Y, 'b', t, y, 'r') Modified: trunk/octave-forge/extra/control-devel/src/slident.cc =================================================================== --- trunk/octave-forge/extra/control-devel/src/slident.cc 2012-03-21 09:48:12 UTC (rev 9988) +++ trunk/octave-forge/extra/control-devel/src/slident.cc 2012-03-21 10:21:26 UTC (rev 9989) @@ -95,7 +95,7 @@ int nargin = args.length (); octave_value_list retval; - if (nargin != 11) + if (nargin != 12) { print_usage (); } @@ -118,17 +118,18 @@ Matrix y = args(0).matrix_value (); Matrix u = args(1).matrix_value (); int nobr = args(2).int_value (); + int nuser = args(3).int_value (); - const int imeth = args(3).int_value (); - const int ialg = args(4).int_value (); - const int ijobd = args(5).int_value (); - const int ibatch = args(6).int_value (); - const int iconct = args(7).int_value (); - const int ictrl = args(8).int_value (); + const int imeth = args(4).int_value (); + const int ialg = args(5).int_value (); + const int ijobd = args(6).int_value (); + const int ibatch = args(7).int_value (); + const int iconct = args(8).int_value (); + const int ictrl = args(9).int_value (); - double rcond = args(9).double_value (); - double tol = args(10).double_value (); - double tolb = args(9).double_value (); // tolb = rcond + double rcond = args(10).double_value (); + double tol = args(11).double_value (); + double tolb = args(10).double_value (); // tolb = rcond switch (imeth) @@ -395,6 +396,14 @@ int rs = 2*(m+l)*nobr; r.resize (rs, rs); + if (nuser > 0) + { + if (nuser < nobr) + n = nuser; + else + error ("ident: 'nuser' invalid"); + } + //////////////////////////////////////////////////////////////////////////////////// // SLICOT IB01BD - estimating system matrices, Kalman gain, and covariances // //////////////////////////////////////////////////////////////////////////////////// This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2012-03-25 14:42:00
|
Revision: 10050 http://octave.svn.sourceforge.net/octave/?rev=10050&view=rev Author: paramaniac Date: 2012-03-25 14:41:48 +0000 (Sun, 25 Mar 2012) Log Message: ----------- control-devel: change ldwork size, CamelCase for examples Modified Paths: -------------- trunk/octave-forge/extra/control-devel/src/slident.cc Added Paths: ----------- trunk/octave-forge/extra/control-devel/devel/Destillation.m trunk/octave-forge/extra/control-devel/devel/GlassFurnace.m trunk/octave-forge/extra/control-devel/devel/PowerPlant.m Removed Paths: ------------- trunk/octave-forge/extra/control-devel/devel/destillation.m trunk/octave-forge/extra/control-devel/devel/glass_furnace.m trunk/octave-forge/extra/control-devel/devel/power_plant.m Copied: trunk/octave-forge/extra/control-devel/devel/Destillation.m (from rev 10049, trunk/octave-forge/extra/control-devel/devel/destillation.m) =================================================================== --- trunk/octave-forge/extra/control-devel/devel/Destillation.m (rev 0) +++ trunk/octave-forge/extra/control-devel/devel/Destillation.m 2012-03-25 14:41:48 UTC (rev 10050) @@ -0,0 +1,81 @@ +%{ +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); +%} + +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); + + +dat = iddata (Y_dest, U_dest) + +[sys, x0] = ident (dat, 5, 4) % s=5, n=4 + + +[y, t] = lsim (sys, U_dest, [], x0); +%[y, t] = lsim (sys, U_dest); + +err = norm (Y_dest - y, 1) / norm (Y_dest, 1) + +figure (1) +plot (t, Y_dest, 'b') +%plot (t, Y_dest, 'b', t, y, 'r') +legend ('y measured', 'y simulated', 'location', 'southeast') + + Copied: trunk/octave-forge/extra/control-devel/devel/GlassFurnace.m (from rev 10049, trunk/octave-forge/extra/control-devel/devel/glass_furnace.m) =================================================================== --- trunk/octave-forge/extra/control-devel/devel/GlassFurnace.m (rev 0) +++ trunk/octave-forge/extra/control-devel/devel/GlassFurnace.m 2012-03-25 14:41:48 UTC (rev 10050) @@ -0,0 +1,66 @@ +%{ +This file describes the data in the glassfurnace.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 glassfurnace (Philips) +3. Sampling time + +4. Number of samples: + 1247 samples +5. Inputs: + a. heating input + b. cooling input + c. heating input +6. Outputs: + a. 6 outputs from temperature sensors in a cross section of the + furnace +7. References: + a. Van Overschee P., De Moor B., N4SID : Subspace Algorithms for + the Identification of Combined Deterministic-Stochastic Systems, + Automatica, Special Issue on Statistical Signal Processing and Control, + Vol. 30, No. 1, 1994, pp. 75-93 + b. Van Overschee P., "Subspace identification : Theory, + Implementation, Application" , Ph.D. Thesis, K.U.Leuven, February 1995. +8. Known properties/peculiarities + +9. Some MATLAB-code to retrieve the data + !gunzip glassfurnace.dat.Z + load glassfurnace.dat + T=glassfurnace(:,1); + U=glassfurnace(:,2:4); + Y=glassfurnace(:,5:10); + +%} + +clear all, close all, clc + +load glassfurnace.dat +T=glassfurnace(:,1); +U=glassfurnace(:,2:4); +Y=glassfurnace(:,5:10); + + +dat = iddata (Y, U) + +[sys, x0] = ident (dat, 10, 5) % s=10, n=5 + + +[y, t] = lsim (sys, U, [], x0); + +err = norm (Y - y, 1) / norm (Y, 1) + +figure (1) +p = columns (Y); +for k = 1 : p + subplot (3, 2, k) + plot (t, Y(:,k), 'b', t, y(:,k), 'r') +endfor +%title ('DaISy: Glass Furnace') +%legend ('y measured', 'y simulated', 'location', 'southeast') + + Copied: trunk/octave-forge/extra/control-devel/devel/PowerPlant.m (from rev 10049, trunk/octave-forge/extra/control-devel/devel/power_plant.m) =================================================================== --- trunk/octave-forge/extra/control-devel/devel/PowerPlant.m (rev 0) +++ trunk/octave-forge/extra/control-devel/devel/PowerPlant.m 2012-03-25 14:41:48 UTC (rev 10050) @@ -0,0 +1,80 @@ +%{ +This file describes the data in the powerplant.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 power plant (Pont-sur-Sambre (France)) of 120 MW +3. Sampling time + 1228.8 sec +4. Number of samples: + 200 samples +5. Inputs: + 1. gas flow + 2. turbine valves opening + 3. super heater spray flow + 4. gas dampers + 5. air flow +6. Outputs: + 1. steam pressure + 2. main stem temperature + 3. reheat steam temperature +7. References: + a. R.P. Guidorzi, P. Rossi, Identification of a power plant from normal + operating records. Automatic control theory and applications (Canada, + Vol 2, pp 63-67, sept 1974. + b. Moonen M., De Moor B., Vandenberghe L., Vandewalle J., On- and + off-line identification of linear state-space models, International + Journal of Control, Vol. 49, Jan. 1989, pp.219-232 +8. Known properties/peculiarities + +9. Some MATLAB-code to retrieve the data + !gunzip powerplant.dat.Z + load powerplant.dat + U=powerplant(:,1:5); + Y=powerplant(:,6:8); + Yr=powerplant(:,9:11); + +%} + +clear all, close all, clc + +load powerplant.dat +U=powerplant(:,1:5); +Y=powerplant(:,6:8); +Yr=powerplant(:,9:11); + +inname = {'gas flow', + 'turbine valves opening', + 'super heater spray flow', + 'gas dampers', + 'air flow'}; + +outname = {'steam pressure', + 'main steam temperature', + 'reheat steam temperature'}; + +tsam = 1228.8; + +dat = iddata (Y, U, tsam, 'outname', outname, 'inname', inname) + +[sys, x0] = ident (dat, 10, 8) % s=10, n=8 + + +[y, t] = lsim (sys, U, [], x0); + +err = norm (Y - y, 1) / norm (Y, 1) + +figure (1) +p = columns (Y); +for k = 1 : p + subplot (3, 1, k) + plot (t, Y(:,k), 'b', t, y(:,k), 'r') +endfor +%title ('DaISy: Power Plant') +%legend ('y measured', 'y simulated', 'location', 'southeast') + + Deleted: trunk/octave-forge/extra/control-devel/devel/destillation.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/destillation.m 2012-03-25 09:01:46 UTC (rev 10049) +++ trunk/octave-forge/extra/control-devel/devel/destillation.m 2012-03-25 14:41:48 UTC (rev 10050) @@ -1,81 +0,0 @@ -%{ -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); -%} - -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); - - -dat = iddata (Y_dest, U_dest) - -[sys, x0] = ident (dat, 5, 4) % s=5, n=4 - - -[y, t] = lsim (sys, U_dest, [], x0); -%[y, t] = lsim (sys, U_dest); - -err = norm (Y_dest - y, 1) / norm (Y_dest, 1) - -figure (1) -plot (t, Y_dest, 'b') -%plot (t, Y_dest, 'b', t, y, 'r') -legend ('y measured', 'y simulated', 'location', 'southeast') - - Deleted: trunk/octave-forge/extra/control-devel/devel/glass_furnace.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/glass_furnace.m 2012-03-25 09:01:46 UTC (rev 10049) +++ trunk/octave-forge/extra/control-devel/devel/glass_furnace.m 2012-03-25 14:41:48 UTC (rev 10050) @@ -1,66 +0,0 @@ -%{ -This file describes the data in the glassfurnace.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 glassfurnace (Philips) -3. Sampling time - -4. Number of samples: - 1247 samples -5. Inputs: - a. heating input - b. cooling input - c. heating input -6. Outputs: - a. 6 outputs from temperature sensors in a cross section of the - furnace -7. References: - a. Van Overschee P., De Moor B., N4SID : Subspace Algorithms for - the Identification of Combined Deterministic-Stochastic Systems, - Automatica, Special Issue on Statistical Signal Processing and Control, - Vol. 30, No. 1, 1994, pp. 75-93 - b. Van Overschee P., "Subspace identification : Theory, - Implementation, Application" , Ph.D. Thesis, K.U.Leuven, February 1995. -8. Known properties/peculiarities - -9. Some MATLAB-code to retrieve the data - !gunzip glassfurnace.dat.Z - load glassfurnace.dat - T=glassfurnace(:,1); - U=glassfurnace(:,2:4); - Y=glassfurnace(:,5:10); - -%} - -clear all, close all, clc - -load glassfurnace.dat -T=glassfurnace(:,1); -U=glassfurnace(:,2:4); -Y=glassfurnace(:,5:10); - - -dat = iddata (Y, U) - -[sys, x0] = ident (dat, 10, 5) % s=10, n=5 - - -[y, t] = lsim (sys, U, [], x0); - -err = norm (Y - y, 1) / norm (Y, 1) - -figure (1) -p = columns (Y); -for k = 1 : p - subplot (3, 2, k) - plot (t, Y(:,k), 'b', t, y(:,k), 'r') -endfor -%title ('DaISy: Glass Furnace') -%legend ('y measured', 'y simulated', 'location', 'southeast') - - Deleted: trunk/octave-forge/extra/control-devel/devel/power_plant.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/power_plant.m 2012-03-25 09:01:46 UTC (rev 10049) +++ trunk/octave-forge/extra/control-devel/devel/power_plant.m 2012-03-25 14:41:48 UTC (rev 10050) @@ -1,80 +0,0 @@ -%{ -This file describes the data in the powerplant.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 power plant (Pont-sur-Sambre (France)) of 120 MW -3. Sampling time - 1228.8 sec -4. Number of samples: - 200 samples -5. Inputs: - 1. gas flow - 2. turbine valves opening - 3. super heater spray flow - 4. gas dampers - 5. air flow -6. Outputs: - 1. steam pressure - 2. main stem temperature - 3. reheat steam temperature -7. References: - a. R.P. Guidorzi, P. Rossi, Identification of a power plant from normal - operating records. Automatic control theory and applications (Canada, - Vol 2, pp 63-67, sept 1974. - b. Moonen M., De Moor B., Vandenberghe L., Vandewalle J., On- and - off-line identification of linear state-space models, International - Journal of Control, Vol. 49, Jan. 1989, pp.219-232 -8. Known properties/peculiarities - -9. Some MATLAB-code to retrieve the data - !gunzip powerplant.dat.Z - load powerplant.dat - U=powerplant(:,1:5); - Y=powerplant(:,6:8); - Yr=powerplant(:,9:11); - -%} - -clear all, close all, clc - -load powerplant.dat -U=powerplant(:,1:5); -Y=powerplant(:,6:8); -Yr=powerplant(:,9:11); - -inname = {'gas flow', - 'turbine valves opening', - 'super heater spray flow', - 'gas dampers', - 'air flow'}; - -outname = {'steam pressure', - 'main steam temperature', - 'reheat steam temperature'}; - -tsam = 1228.8; - -dat = iddata (Y, U, tsam, 'outname', outname, 'inname', inname) - -[sys, x0] = ident (dat, 10, 8) % s=10, n=8 - - -[y, t] = lsim (sys, U, [], x0); - -err = norm (Y - y, 1) / norm (Y, 1) - -figure (1) -p = columns (Y); -for k = 1 : p - subplot (3, 1, k) - plot (t, Y(:,k), 'b', t, y(:,k), 'r') -endfor -%title ('DaISy: Power Plant') -%legend ('y measured', 'y simulated', 'location', 'southeast') - - Modified: trunk/octave-forge/extra/control-devel/src/slident.cc =================================================================== --- trunk/octave-forge/extra/control-devel/src/slident.cc 2012-03-25 09:01:46 UTC (rev 10049) +++ trunk/octave-forge/extra/control-devel/src/slident.cc 2012-03-25 14:41:48 UTC (rev 10050) @@ -314,7 +314,8 @@ //////////////////////////////////////////////////////////////////// // TO BE REMOVED !!! //////////////////////////////////////////////////////////////////// -ldwork = 1000000; +ldwork *= 2; +//ldwork = 1000000; //////////////////////////////////////////////////////////////////// /* This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2012-03-29 15:07:27
|
Revision: 10093 http://octave.svn.sourceforge.net/octave/?rev=10093&view=rev Author: paramaniac Date: 2012-03-29 15:07:17 +0000 (Thu, 29 Mar 2012) Log Message: ----------- control-devel: debugging stuff for ib01ad Modified Paths: -------------- trunk/octave-forge/extra/control-devel/inst/@iddata/set.m trunk/octave-forge/extra/control-devel/src/devel_slicot_functions.cc Added Paths: ----------- trunk/octave-forge/extra/control-devel/devel/PowerPlantTest.m trunk/octave-forge/extra/control-devel/devel/identtest.m trunk/octave-forge/extra/control-devel/src/slitest.cc Added: trunk/octave-forge/extra/control-devel/devel/PowerPlantTest.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/PowerPlantTest.m (rev 0) +++ trunk/octave-forge/extra/control-devel/devel/PowerPlantTest.m 2012-03-29 15:07:17 UTC (rev 10093) @@ -0,0 +1,83 @@ +%{ +This file describes the data in the powerplant.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 power plant (Pont-sur-Sambre (France)) of 120 MW +3. Sampling time + 1228.8 sec +4. Number of samples: + 200 samples +5. Inputs: + 1. gas flow + 2. turbine valves opening + 3. super heater spray flow + 4. gas dampers + 5. air flow +6. Outputs: + 1. steam pressure + 2. main stem temperature + 3. reheat steam temperature +7. References: + a. R.P. Guidorzi, P. Rossi, Identification of a power plant from normal + operating records. Automatic control theory and applications (Canada, + Vol 2, pp 63-67, sept 1974. + b. Moonen M., De Moor B., Vandenberghe L., Vandewalle J., On- and + off-line identification of linear state-space models, International + Journal of Control, Vol. 49, Jan. 1989, pp.219-232 +8. Known properties/peculiarities + +9. Some MATLAB-code to retrieve the data + !gunzip powerplant.dat.Z + load powerplant.dat + U=powerplant(:,1:5); + Y=powerplant(:,6:8); + Yr=powerplant(:,9:11); + +%} + +clear all, close all, clc + +load powerplant.dat +U=powerplant(:,1:5); +Y=powerplant(:,6:8); +Yr=powerplant(:,9:11); + +inname = {'gas flow', + 'turbine valves opening', + 'super heater spray flow', + 'gas dampers', + 'air flow'}; + +outname = {'steam pressure', + 'main steam temperature', + 'reheat steam temperature'}; + +tsam = 1228.8; + +dat = iddata (Y, U, tsam, 'outname', outname, 'inname', inname); + + +% ldwork = [401, 802, 1203, 1604] +% warning: implicit conversion from real matrix to real scalar + +ldwork = [802, 1203, 1604] + +r = arrayfun (@(x) identtest (dat, 10, 8, x), ldwork, 'uniformoutput', false); + + % s=10, n=8 + +l = length (ldwork); + +err = cell (l-1, 1); + +for k = 2 : l + err(k-1) = norm (abs(r{1}) - abs(r{k}), 1); + % err(k-1) = norm (abs(abs(r{1}) - abs(r{k})), 1); +endfor + +err Added: trunk/octave-forge/extra/control-devel/devel/identtest.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/identtest.m (rev 0) +++ trunk/octave-forge/extra/control-devel/devel/identtest.m 2012-03-29 15:07:17 UTC (rev 10093) @@ -0,0 +1,52 @@ +function [r, sv, n] = identtest (dat, s = [], n = [], ldwork) + + + + %nobr = 15; + meth = 2; + alg = 0; + jobd = 1; + batch = 3; + conct = 1; + ctrl = 0; %1; + rcond = 0.0; + tol = -1.0; % 0; + + [ns, l, m, e] = size (dat); + + if (isempty (s) && isempty (n)) + nsmp = ns(1); + nobr = fix ((nsmp+1)/(2*(m+l+1))); + ctrl = 0; # confirm system order estimate + n = 0; + % nsmp >= 2*(m+l+1)*nobr - 1 + % nobr <= (nsmp+1)/(2*(m+l+1)) + elseif (isempty (s)) + s = min (2*n, n+10); + nsmp = ns(1); + nobr = fix ((nsmp+1)/(2*(m+l+1))); + nobr = min (nobr, s); + ctrl = 1; # no confirmation + elseif (isempty (n)) + nobr = s; + ctrl = 0; # confirm system order estimate + n = 0; + else # s & n non-empty + nsmp = ns(1); + nobr = fix ((nsmp+1)/(2*(m+l+1))); + if (s > nobr) + error ("ident: s > nobr"); + endif + nobr = s; + ctrl = 1; + ## TODO: specify n for IB01BD + endif + + %nsmp = ns(1) + %nobr = fix ((nsmp+1)/(2*(m+l+1))) + % nsmp >= 2*(m+l+1)*nobr - 1 + % nobr <= (nsmp+1)/(2*(m+l+1)) +%nobr = 10 + [r, sv, n] = slitest (dat.y{1}, dat.u{1}, nobr, n, meth, alg, jobd, batch, conct, ctrl, rcond, tol, ldwork); + +endfunction Modified: trunk/octave-forge/extra/control-devel/inst/@iddata/set.m =================================================================== --- trunk/octave-forge/extra/control-devel/inst/@iddata/set.m 2012-03-29 14:33:05 UTC (rev 10092) +++ trunk/octave-forge/extra/control-devel/inst/@iddata/set.m 2012-03-29 15:07:17 UTC (rev 10093) @@ -84,8 +84,6 @@ dat.outunit = __adjust_labels__ (val, p); case {"inunit", "inputunit"} dat.inunit = __adjust_labels__ (val, m); - case {"tsam", "ts"} - dat.tsam; case {"timeunit"} dat.timeunit case {"expname", "experimentname"} Modified: trunk/octave-forge/extra/control-devel/src/devel_slicot_functions.cc =================================================================== --- trunk/octave-forge/extra/control-devel/src/devel_slicot_functions.cc 2012-03-29 14:33:05 UTC (rev 10092) +++ trunk/octave-forge/extra/control-devel/src/devel_slicot_functions.cc 2012-03-29 15:07:17 UTC (rev 10093) @@ -1,2 +1,3 @@ // #include "slib01ad.cc" // preprocess the input-output data #include "slident.cc" // system identification +#include "slitest.cc" // debug system identification Added: trunk/octave-forge/extra/control-devel/src/slitest.cc =================================================================== --- trunk/octave-forge/extra/control-devel/src/slitest.cc (rev 0) +++ trunk/octave-forge/extra/control-devel/src/slitest.cc 2012-03-29 15:07:17 UTC (rev 10093) @@ -0,0 +1,418 @@ +/* + +Copyright (C) 2012 Lukas F. Reichlin + +This file is part of LTI Syncope. + +LTI Syncope is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +LTI Syncope is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with LTI Syncope. If not, see <http://www.gnu.org/licenses/>. + +SLICOT system identification +Uses SLICOT IB01AD, IB01BD and IB01CD by courtesy of NICONET e.V. +<http://www.slicot.org> + +Author: Lukas Reichlin <luk...@gm...> +Created: March 2012 +Version: 0.1 + +*/ + +#include <octave/oct.h> +#include <f77-fcn.h> +#include "common.h" + +extern "C" +{ + int F77_FUNC (ib01ad, IB01AD) + (char& METH, char& ALG, char& JOBD, + char& BATCH, char& CONCT, char& CTRL, + int& NOBR, int& M, int& L, + int& NSMP, + double* U, int& LDU, + double* Y, int& LDY, + int& N, + double* R, int& LDR, + double* SV, + double& RCOND, double& TOL, + int* IWORK, + double* DWORK, int& LDWORK, + int& IWARN, int& INFO); + + int F77_FUNC (ib01bd, IB01BD) + (char& METH, char& JOB, char& JOBCK, + int& NOBR, int& N, int& M, int& L, + int& NSMPL, + double* R, int& LDR, + double* A, int& LDA, + double* C, int& LDC, + double* B, int& LDB, + double* D, int& LDD, + double* Q, int& LDQ, + double* RY, int& LDRY, + double* S, int& LDS, + double* K, int& LDK, + double& TOL, + int* IWORK, + double* DWORK, int& LDWORK, + bool* BWORK, + int& IWARN, int& INFO); + + int F77_FUNC (ib01cd, IB01CD) + (char& JOBX0, char& COMUSE, char& JOB, + int& N, int& M, int& L, + int& NSMP, + double* A, int& LDA, + double* B, int& LDB, + double* C, int& LDC, + double* D, int& LDD, + double* U, int& LDU, + double* Y, int& LDY, + double* X0, + double* V, int& LDV, + double& TOL, + int* IWORK, + double* DWORK, int& LDWORK, + int& IWARN, int& INFO); +} + +// PKG_ADD: autoload ("slitest", "devel_slicot_functions.oct"); +DEFUN_DLD (slitest, args, nargout, + "-*- texinfo -*-\n\ +Slicot IB01AD Release 5.0\n\ +No argument checking.\n\ +For internal use only.") +{ + int nargin = args.length (); + octave_value_list retval; + + if (nargin != 13) + { + print_usage (); + } + else + { +//////////////////////////////////////////////////////////////////////////////////// +// SLICOT IB01AD - preprocess the input-output data // +//////////////////////////////////////////////////////////////////////////////////// + + // arguments in + char meth; + char alg; + char jobd; + char batch; + char conct; + char ctrl; + char metha; + char jobda; + + Matrix y = args(0).matrix_value (); + Matrix u = args(1).matrix_value (); + int nobr = args(2).int_value (); + int nuser = args(3).int_value (); + + const int imeth = args(4).int_value (); + const int ialg = args(5).int_value (); + const int ijobd = args(6).int_value (); + const int ibatch = args(7).int_value (); + const int iconct = args(8).int_value (); + const int ictrl = args(9).int_value (); + + double rcond = args(10).double_value (); + double tol = args(11).double_value (); + double tolb = args(10).double_value (); // tolb = rcond + + int ldwork = args(12).int_value (); + + + switch (imeth) + { + case 0: + meth = 'M'; + metha = 'M'; + break; + case 1: + meth = 'N'; + metha = 'N'; + break; + case 2: + meth = 'C'; + metha = 'N'; // no typo here + break; + default: + error ("slib01ad: argument 'meth' invalid"); + } + + switch (ialg) + { + case 0: + alg = 'C'; + break; + case 1: + alg = 'F'; + break; + case 2: + alg = 'Q'; + break; + default: + error ("slib01ad: argument 'alg' invalid"); + } + + if (meth == 'C') + jobd = 'N'; + else if (ijobd == 0) + jobd = 'M'; + else + jobd = 'N'; + + switch (ibatch) + { + case 0: + batch = 'F'; + break; + case 1: + batch = 'I'; + break; + case 2: + batch = 'L'; + break; + case 3: + batch = 'O'; + break; + default: + error ("slib01ad: argument 'batch' invalid"); + } + + if (iconct == 0) + conct = 'C'; + else + conct = 'N'; + + if (ictrl == 0) + ctrl = 'C'; + else + ctrl = 'N'; + + + int m = u.columns (); // m: number of inputs + int l = y.columns (); // l: number of outputs + int nsmp = y.rows (); // nsmp: number of samples + // y.rows == u.rows is checked by iddata class + // TODO: check minimal nsmp size + + if (batch == 'O') + { + if (nsmp < 2*(m+l+1)*nobr - 1) + error ("slident: require NSMP >= 2*(M+L+1)*NOBR - 1"); + } + else + { + if (nsmp < 2*nobr) + error ("slident: require NSMP >= 2*NOBR"); + } + + int ldu; + + if (m == 0) + ldu = 1; + else // m > 0 + ldu = nsmp; + + int ldy = nsmp; + + // arguments out + int n; + int ldr; + + if (metha == 'M' && jobd == 'M') + ldr = max (2*(m+l)*nobr, 3*m*nobr); + else if (metha == 'N' || (metha == 'M' && jobd == 'N')) + ldr = 2*(m+l)*nobr; + else + error ("slib01ad: could not handle 'ldr' case"); + + Matrix r (ldr, 2*(m+l)*nobr); + ColumnVector sv (l*nobr); + + // workspace + int liwork; + + if (metha == 'N') // if METH = 'N' + liwork = (m+l)*nobr; + else if (alg == 'F') // if METH = 'M' and ALG = 'F' + liwork = m+l; + else // if METH = 'M' and ALG = 'C' or 'Q' + liwork = 0; + + // TODO: Handle 'k' for DWORK +/* + int ldwork; + + if (alg == 'C') + { + if (batch == 'F' || batch == 'I') + { + if (conct == 'C') + ldwork = (4*nobr-2)*(m+l); + else // (conct == 'N') + ldwork = 1; + } + else if (metha == 'M') // && (batch == 'L' || batch == 'O') + { + if (conct == 'C' && batch == 'L') + ldwork = max ((4*nobr-2)*(m+l), 5*l*nobr); + else if (jobd == 'M') + ldwork = max ((2*m-1)*nobr, (m+l)*nobr, 5*l*nobr); + else // (jobd == 'N') + ldwork = 5*l*nobr; + } + else // meth == 'N' && (batch == 'L' || batch == 'O') + { + ldwork = 5*(m+l)*nobr + 1; + } + } + else if (alg == 'F') + { + if (batch != 'O' && conct == 'C') + ldwork = (m+l)*2*nobr*(m+l+3); + else if (batch == 'F' || batch == 'I') // && conct == 'N' + ldwork = (m+l)*2*nobr*(m+l+1); + else // (batch == 'L' || '0' && conct == 'N') + ldwork = (m+l)*4*nobr*(m+l+1)+(m+l)*2*nobr; + } + else // (alg == 'Q') + { + int ns = nsmp - 2*nobr + 1; + + if (ldr >= ns && batch == 'F') + { + ldwork = 4*(m+l)*nobr; + } + else if (ldr >= ns && batch == 'O') + { + if (metha == 'M') + ldwork = max (4*(m+l)*nobr, 5*l*nobr); + else // (meth == 'N') + ldwork = 5*(m+l)*nobr + 1; + } + else if (conct == 'C' && (batch == 'I' || batch == 'L')) + { + ldwork = 4*(nobr+1)*(m+l)*nobr; + } + else // if ALG = 'Q', (BATCH = 'F' or 'O', and LDR < NS), or (BATCH = 'I' or 'L' and CONCT = 'N') + { + ldwork = 6*(m+l)*nobr; + } + } +//////////////////////////////////////////////////////////////////// +// TO BE REMOVED !!! +//////////////////////////////////////////////////////////////////// +ldwork *= 2; +//ldwork = 1000000; +//////////////////////////////////////////////////////////////////// + +*/ + +/* +IB01AD.f Lines 291-195: +c the workspace used for alg = 'q' is +c ldrwrk*2*(m+l)*nobr + 4*(m+l)*nobr, +c where ldrwrk = ldwork/(2*(m+l)*nobr) - 2; recommended +c value ldrwrk = ns, assuming a large enough cache size. +c for good performance, ldwork should be larger. + +somehow ldrwrk and ldwork must have been mixed up here + +*/ + + + OCTAVE_LOCAL_BUFFER (int, iwork, liwork); + OCTAVE_LOCAL_BUFFER (double, dwork, ldwork); + + // error indicators + int iwarn = 0; + int info = 0; + + + // SLICOT routine IB01AD + F77_XFCN (ib01ad, IB01AD, + (metha, alg, jobd, + batch, conct, ctrl, + nobr, m, l, + nsmp, + u.fortran_vec (), ldu, + y.fortran_vec (), ldy, + n, + r.fortran_vec (), ldr, + sv.fortran_vec (), + rcond, tol, + iwork, + dwork, ldwork, + iwarn, info)); + + + if (f77_exception_encountered) + error ("ident: exception in SLICOT subroutine IB01AD"); + + static const char* err_msg[] = { + "0: OK", + "1: a fast algorithm was requested (ALG = 'C', or 'F') " + "in sequential data processing, but it failed; the " + "routine can be repeatedly called again using the " + "standard QR algorithm", + "2: the singular value decomposition (SVD) algorithm did " + "not converge"}; + + static const char* warn_msg[] = { + "0: OK", + "1: the number of 100 cycles in sequential data " + "processing has been exhausted without signaling " + "that the last block of data was get; the cycle " + "counter was reinitialized", + "2: a fast algorithm was requested (ALG = 'C' or 'F'), " + "but it failed, and the QR algorithm was then used " + "(non-sequential data processing)", + "3: all singular values were exactly zero, hence N = 0 " + "(both input and output were identically zero)", + "4: the least squares problems with coefficient matrix " + "U_f, used for computing the weighted oblique " + "projection (for METH = 'N'), have a rank-deficient " + "coefficient matrix", + "5: the least squares problem with coefficient matrix " + "r_1 [6], used for computing the weighted oblique " + "projection (for METH = 'N'), has a rank-deficient " + "coefficient matrix"}; + + + error_msg ("ident", info, 2, err_msg); + warning_msg ("ident", iwarn, 5, warn_msg); + + + // resize + int rs = 2*(m+l)*nobr; + r.resize (rs, rs); + + if (nuser > 0) + { + if (nuser < nobr) + n = nuser; + else + error ("ident: 'nuser' invalid"); + } + + retval(0) = r; + retval(1) = sv; + retval(2) = octave_value (n); + } + + return retval; +} This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2012-03-29 15:33:52
|
Revision: 10097 http://octave.svn.sourceforge.net/octave/?rev=10097&view=rev Author: paramaniac Date: 2012-03-29 15:33:41 +0000 (Thu, 29 Mar 2012) Log Message: ----------- control-devel: update debugging stuff for ib01ad Modified Paths: -------------- trunk/octave-forge/extra/control-devel/devel/identtest.m trunk/octave-forge/extra/control-devel/src/slident.cc trunk/octave-forge/extra/control-devel/src/slitest.cc Modified: trunk/octave-forge/extra/control-devel/devel/identtest.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/identtest.m 2012-03-29 15:21:47 UTC (rev 10096) +++ trunk/octave-forge/extra/control-devel/devel/identtest.m 2012-03-29 15:33:41 UTC (rev 10097) @@ -1,52 +1,5 @@ function [r, sv, n] = identtest (dat, s = [], n = [], ldwork) + [r, sv, n] = slitest (dat.y{1}, dat.u{1}, s, ldwork); - - %nobr = 15; - meth = 2; - alg = 0; - jobd = 1; - batch = 3; - conct = 1; - ctrl = 0; %1; - rcond = 0.0; - tol = -1.0; % 0; - - [ns, l, m, e] = size (dat); - - if (isempty (s) && isempty (n)) - nsmp = ns(1); - nobr = fix ((nsmp+1)/(2*(m+l+1))); - ctrl = 0; # confirm system order estimate - n = 0; - % nsmp >= 2*(m+l+1)*nobr - 1 - % nobr <= (nsmp+1)/(2*(m+l+1)) - elseif (isempty (s)) - s = min (2*n, n+10); - nsmp = ns(1); - nobr = fix ((nsmp+1)/(2*(m+l+1))); - nobr = min (nobr, s); - ctrl = 1; # no confirmation - elseif (isempty (n)) - nobr = s; - ctrl = 0; # confirm system order estimate - n = 0; - else # s & n non-empty - nsmp = ns(1); - nobr = fix ((nsmp+1)/(2*(m+l+1))); - if (s > nobr) - error ("ident: s > nobr"); - endif - nobr = s; - ctrl = 1; - ## TODO: specify n for IB01BD - endif - - %nsmp = ns(1) - %nobr = fix ((nsmp+1)/(2*(m+l+1))) - % nsmp >= 2*(m+l+1)*nobr - 1 - % nobr <= (nsmp+1)/(2*(m+l+1)) -%nobr = 10 - [r, sv, n] = slitest (dat.y{1}, dat.u{1}, nobr, n, meth, alg, jobd, batch, conct, ctrl, rcond, tol, ldwork); - endfunction Modified: trunk/octave-forge/extra/control-devel/src/slident.cc =================================================================== --- trunk/octave-forge/extra/control-devel/src/slident.cc 2012-03-29 15:21:47 UTC (rev 10096) +++ trunk/octave-forge/extra/control-devel/src/slident.cc 2012-03-29 15:33:41 UTC (rev 10097) @@ -113,7 +113,7 @@ char conct; char ctrl; char metha; - char jobda; + char jobda; // ??? unused Matrix y = args(0).matrix_value (); Matrix u = args(1).matrix_value (); Modified: trunk/octave-forge/extra/control-devel/src/slitest.cc =================================================================== --- trunk/octave-forge/extra/control-devel/src/slitest.cc 2012-03-29 15:21:47 UTC (rev 10096) +++ trunk/octave-forge/extra/control-devel/src/slitest.cc 2012-03-29 15:33:41 UTC (rev 10097) @@ -95,7 +95,7 @@ int nargin = args.length (); octave_value_list retval; - if (nargin != 13) + if (nargin != 4) { print_usage (); } @@ -106,127 +106,36 @@ //////////////////////////////////////////////////////////////////////////////////// // arguments in - char meth; - char alg; - char jobd; - char batch; - char conct; - char ctrl; - char metha; - char jobda; + char meth = 'N'; // <--- not used, use metha + char alg = 'C'; + char jobd = 'N'; + char batch = 'O'; + char conct = 'N'; + char ctrl = 'N'; + + char metha = 'N'; + + // ??? char jobda ; + Matrix y = args(0).matrix_value (); Matrix u = args(1).matrix_value (); int nobr = args(2).int_value (); - int nuser = args(3).int_value (); - const int imeth = args(4).int_value (); - const int ialg = args(5).int_value (); - const int ijobd = args(6).int_value (); - const int ibatch = args(7).int_value (); - const int iconct = args(8).int_value (); - const int ictrl = args(9).int_value (); + double rcond = 0.0; + double tol = -1.0; - double rcond = args(10).double_value (); - double tol = args(11).double_value (); - double tolb = args(10).double_value (); // tolb = rcond - - int ldwork = args(12).int_value (); + int ldwork = args(3).int_value (); - - switch (imeth) - { - case 0: - meth = 'M'; - metha = 'M'; - break; - case 1: - meth = 'N'; - metha = 'N'; - break; - case 2: - meth = 'C'; - metha = 'N'; // no typo here - break; - default: - error ("slib01ad: argument 'meth' invalid"); - } - switch (ialg) - { - case 0: - alg = 'C'; - break; - case 1: - alg = 'F'; - break; - case 2: - alg = 'Q'; - break; - default: - error ("slib01ad: argument 'alg' invalid"); - } - - if (meth == 'C') - jobd = 'N'; - else if (ijobd == 0) - jobd = 'M'; - else - jobd = 'N'; - - switch (ibatch) - { - case 0: - batch = 'F'; - break; - case 1: - batch = 'I'; - break; - case 2: - batch = 'L'; - break; - case 3: - batch = 'O'; - break; - default: - error ("slib01ad: argument 'batch' invalid"); - } - if (iconct == 0) - conct = 'C'; - else - conct = 'N'; - - if (ictrl == 0) - ctrl = 'C'; - else - ctrl = 'N'; - - int m = u.columns (); // m: number of inputs int l = y.columns (); // l: number of outputs int nsmp = y.rows (); // nsmp: number of samples // y.rows == u.rows is checked by iddata class // TODO: check minimal nsmp size - - if (batch == 'O') - { - if (nsmp < 2*(m+l+1)*nobr - 1) - error ("slident: require NSMP >= 2*(M+L+1)*NOBR - 1"); - } - else - { - if (nsmp < 2*nobr) - error ("slident: require NSMP >= 2*NOBR"); - } - - int ldu; - - if (m == 0) - ldu = 1; - else // m > 0 - ldu = nsmp; - + + int ldu = nsmp; int ldy = nsmp; // arguments out @@ -322,19 +231,7 @@ */ -/* -IB01AD.f Lines 291-195: -c the workspace used for alg = 'q' is -c ldrwrk*2*(m+l)*nobr + 4*(m+l)*nobr, -c where ldrwrk = ldwork/(2*(m+l)*nobr) - 2; recommended -c value ldrwrk = ns, assuming a large enough cache size. -c for good performance, ldwork should be larger. -somehow ldrwrk and ldwork must have been mixed up here - -*/ - - OCTAVE_LOCAL_BUFFER (int, iwork, liwork); OCTAVE_LOCAL_BUFFER (double, dwork, ldwork); @@ -401,13 +298,6 @@ int rs = 2*(m+l)*nobr; r.resize (rs, rs); - if (nuser > 0) - { - if (nuser < nobr) - n = nuser; - else - error ("ident: 'nuser' invalid"); - } retval(0) = r; retval(1) = sv; This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2012-03-29 15:45:36
|
Revision: 10098 http://octave.svn.sourceforge.net/octave/?rev=10098&view=rev Author: paramaniac Date: 2012-03-29 15:45:27 +0000 (Thu, 29 Mar 2012) Log Message: ----------- control-devel: update debugging stuff for ib01ad (2) Modified Paths: -------------- trunk/octave-forge/extra/control-devel/devel/PowerPlantTest.m trunk/octave-forge/extra/control-devel/src/slitest.cc Modified: trunk/octave-forge/extra/control-devel/devel/PowerPlantTest.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/PowerPlantTest.m 2012-03-29 15:33:41 UTC (rev 10097) +++ trunk/octave-forge/extra/control-devel/devel/PowerPlantTest.m 2012-03-29 15:45:27 UTC (rev 10098) @@ -77,6 +77,7 @@ for k = 2 : l err(k-1) = norm (abs(r{1}) - abs(r{k}), 1); + % err(k-1) = norm (r{1} - r{k}, 1); % err(k-1) = norm (abs(abs(r{1}) - abs(r{k})), 1); endfor Modified: trunk/octave-forge/extra/control-devel/src/slitest.cc =================================================================== --- trunk/octave-forge/extra/control-devel/src/slitest.cc 2012-03-29 15:33:41 UTC (rev 10097) +++ trunk/octave-forge/extra/control-devel/src/slitest.cc 2012-03-29 15:45:27 UTC (rev 10098) @@ -106,15 +106,15 @@ //////////////////////////////////////////////////////////////////////////////////// // arguments in - char meth = 'N'; // <--- not used, use metha - char alg = 'C'; - char jobd = 'N'; + char meth = 'M'; // <--- not used, use metha + char alg = 'Q'; + char jobd = 'M'; char batch = 'O'; char conct = 'N'; char ctrl = 'N'; - char metha = 'N'; + char metha = 'M'; // ??? char jobda ; This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2012-04-02 11:36:23
|
Revision: 10130 http://octave.svn.sourceforge.net/octave/?rev=10130&view=rev Author: paramaniac Date: 2012-04-02 11:36:12 +0000 (Mon, 02 Apr 2012) Log Message: ----------- control-devel: improve plot for multi-experiment iddata Modified Paths: -------------- trunk/octave-forge/extra/control-devel/devel/PowerPlantTest.m trunk/octave-forge/extra/control-devel/inst/@iddata/plot.m Modified: trunk/octave-forge/extra/control-devel/devel/PowerPlantTest.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/PowerPlantTest.m 2012-04-02 11:34:53 UTC (rev 10129) +++ trunk/octave-forge/extra/control-devel/devel/PowerPlantTest.m 2012-04-02 11:36:12 UTC (rev 10130) @@ -65,7 +65,7 @@ % ldwork = [401, 802, 1203, 1604] % warning: implicit conversion from real matrix to real scalar -ldwork = [802, 1203, 1604] +ldwork = [802, 1203, 1604, 3000, 10000] r = arrayfun (@(x) identtest (dat, 10, 8, x), ldwork, 'uniformoutput', false); Modified: trunk/octave-forge/extra/control-devel/inst/@iddata/plot.m =================================================================== --- trunk/octave-forge/extra/control-devel/inst/@iddata/plot.m 2012-04-02 11:34:53 UTC (rev 10129) +++ trunk/octave-forge/extra/control-devel/inst/@iddata/plot.m 2012-04-02 11:36:12 UTC (rev 10130) @@ -26,21 +26,32 @@ function plot (dat) - [n, p, m, e] = size (dat) + [n, p, m, e] = size (dat); + expname = __labels__ (dat.expname, "exp"); if (m == 0) # time series for k = 1 : e + if (k > 1) + pause + endif plot (dat.y{k}) - hold on + title (expname{k}) + % hold on endfor else # inputs present for k = 1 : e + if (k > 1) + pause + endif subplot (2, 1, 1) plot (dat.y{k}) - hold on + title (expname{k}) + legend (__labels__ (dat.outname, "y"){:}) + % hold on subplot (2, 1, 2) stairs (dat.u{k}) - hold on + legend (__labels__ (dat.inname, "u"){:}) + % hold on endfor endif This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2012-04-04 09:19:51
|
Revision: 10143 http://octave.svn.sourceforge.net/octave/?rev=10143&view=rev Author: paramaniac Date: 2012-04-04 09:19:42 +0000 (Wed, 04 Apr 2012) Log Message: ----------- control-devel: touch up fft method (3) Modified Paths: -------------- trunk/octave-forge/extra/control-devel/inst/@iddata/fft.m Added Paths: ----------- trunk/octave-forge/extra/control-devel/devel/fixtest.m Added: trunk/octave-forge/extra/control-devel/devel/fixtest.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/fixtest.m (rev 0) +++ trunk/octave-forge/extra/control-devel/devel/fixtest.m 2012-04-04 09:19:42 UTC (rev 10143) @@ -0,0 +1,34 @@ +% Extract FFT +for N = 1:500 + + if (rem (N, 2)) % odd + n1 = (N+1)/2; + else % even + n1 = N/2+1; + endif + + n2 = fix (N/2) + 1; + + if (n1 != n2) + warning ("FFT %d: n1=%d, n2=%d", N, n1, n2); + endif + +endfor + + +% Frequency Vector +for N = 1:500 + + if (rem (N, 2)) % odd + n1 = (N-1)/2; + else % even + n1 = N/2; + endif + + n2 = fix (N/2); + + if (n1 != n2) + warning ("W %d: n1=%d, n2=%d", N, n1, n2); + endif + +endfor Modified: trunk/octave-forge/extra/control-devel/inst/@iddata/fft.m =================================================================== --- trunk/octave-forge/extra/control-devel/inst/@iddata/fft.m 2012-04-04 09:09:45 UTC (rev 10142) +++ trunk/octave-forge/extra/control-devel/inst/@iddata/fft.m 2012-04-04 09:19:42 UTC (rev 10143) @@ -53,6 +53,8 @@ dat.y = cellfun (@(y, n) fft (y, n)(1:fix(n/2)+1, :) / sqrt (n), dat.y, n, "uniformoutput", false); dat.u = cellfun (@(u, n) fft (u, n)(1:fix(n/2)+1, :) / sqrt (n), dat.u, n, "uniformoutput", false); + % w = (0:fix(n/2)) * (2*pi/tsam/n) + dat.timedomain = false; % dat.y = cellfun (@(y) fft (y, n), dat.y, "uniformoutput", false); This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2012-04-11 07:20:55
|
Revision: 10184 http://octave.svn.sourceforge.net/octave/?rev=10184&view=rev Author: paramaniac Date: 2012-04-11 07:20:44 +0000 (Wed, 11 Apr 2012) Log Message: ----------- control-devel: fix corner case in detrend method (set experiments with only 1 sample to zero) Modified Paths: -------------- trunk/octave-forge/extra/control-devel/inst/@iddata/detrend.m Added Paths: ----------- trunk/octave-forge/extra/control-devel/devel/test_frd2iddata.m Added: trunk/octave-forge/extra/control-devel/devel/test_frd2iddata.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/test_frd2iddata.m (rev 0) +++ trunk/octave-forge/extra/control-devel/devel/test_frd2iddata.m 2012-04-11 07:20:44 UTC (rev 10184) @@ -0,0 +1,15 @@ +sys = ss (-2,3,4,5) + +H = idfrd (sys) + +H.frequency +H.responsedata + + +dat = iddata (H) +dat.y +H.responsedata + +dat.u % alles 1! +dat.frequency +H.frequency \ No newline at end of file Modified: trunk/octave-forge/extra/control-devel/inst/@iddata/detrend.m =================================================================== --- trunk/octave-forge/extra/control-devel/inst/@iddata/detrend.m 2012-04-11 05:23:52 UTC (rev 10183) +++ trunk/octave-forge/extra/control-devel/inst/@iddata/detrend.m 2012-04-11 07:20:44 UTC (rev 10184) @@ -38,9 +38,18 @@ error ("iddata: detrend: second argument must be a positve integer"); endif + [n, p, m] = size (dat); + dat.y = cellfun (@(y) detrend (y, ord), dat.y, "uniformoutput", false); dat.u = cellfun (@(u) detrend (u, ord), dat.u, "uniformoutput", false); + ## if a MIMO experiment has only 1 sample, detrend works + ## row-wisely instead of column-wisely + ## therefore we set these experiments to zero + idx = (n == 1); + dat.y(idx) = zeros (1, p); + dat.u(idx) = zeros (1, m); + endfunction This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2012-04-17 13:34:02
|
Revision: 10263 http://octave.svn.sourceforge.net/octave/?rev=10263&view=rev Author: paramaniac Date: 2012-04-17 13:33:51 +0000 (Tue, 17 Apr 2012) Log Message: ----------- control-devel: add ifft method to iddata class Modified Paths: -------------- trunk/octave-forge/extra/control-devel/devel/PowerPlantTest.m trunk/octave-forge/extra/control-devel/inst/@iddata/fft.m trunk/octave-forge/extra/control-devel/src/slitest.cc Added Paths: ----------- trunk/octave-forge/extra/control-devel/devel/PowerPlantFFT.m trunk/octave-forge/extra/control-devel/inst/@iddata/ifft.m Added: trunk/octave-forge/extra/control-devel/devel/PowerPlantFFT.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/PowerPlantFFT.m (rev 0) +++ trunk/octave-forge/extra/control-devel/devel/PowerPlantFFT.m 2012-04-17 13:33:51 UTC (rev 10263) @@ -0,0 +1,68 @@ +%{ +This file describes the data in the powerplant.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 power plant (Pont-sur-Sambre (France)) of 120 MW +3. Sampling time + 1228.8 sec +4. Number of samples: + 200 samples +5. Inputs: + 1. gas flow + 2. turbine valves opening + 3. super heater spray flow + 4. gas dampers + 5. air flow +6. Outputs: + 1. steam pressure + 2. main stem temperature + 3. reheat steam temperature +7. References: + a. R.P. Guidorzi, P. Rossi, Identification of a power plant from normal + operating records. Automatic control theory and applications (Canada, + Vol 2, pp 63-67, sept 1974. + b. Moonen M., De Moor B., Vandenberghe L., Vandewalle J., On- and + off-line identification of linear state-space models, International + Journal of Control, Vol. 49, Jan. 1989, pp.219-232 +8. Known properties/peculiarities + +9. Some MATLAB-code to retrieve the data + !gunzip powerplant.dat.Z + load powerplant.dat + U=powerplant(:,1:5); + Y=powerplant(:,6:8); + Yr=powerplant(:,9:11); + +%} + +clear all, close all, clc + +load powerplant.dat +U=powerplant(:,1:5); +Y=powerplant(:,6:8); +Yr=powerplant(:,9:11); + +inname = {'gas flow', + 'turbine valves opening', + 'super heater spray flow', + 'gas dampers', + 'air flow'}; + +outname = {'steam pressure', + 'main steam temperature', + 'reheat steam temperature'}; + +tsam = 1228.8; + +dat = iddata (Y, U, tsam, 'outputname', outname, 'inputname', inname) + + +a = fft (dat) + +b = ifft (a) + Modified: trunk/octave-forge/extra/control-devel/devel/PowerPlantTest.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/PowerPlantTest.m 2012-04-17 12:03:41 UTC (rev 10262) +++ trunk/octave-forge/extra/control-devel/devel/PowerPlantTest.m 2012-04-17 13:33:51 UTC (rev 10263) @@ -65,7 +65,7 @@ % ldwork = [401, 802, 1203, 1604] % warning: implicit conversion from real matrix to real scalar -ldwork = [802, 1203, 1604, 3000, 10000] +ldwork = [802, 1203, 1604, 3000, 1e4, 1e5] r = arrayfun (@(x) identtest (dat, 10, 8, x), ldwork, 'uniformoutput', false); Modified: trunk/octave-forge/extra/control-devel/inst/@iddata/fft.m =================================================================== --- trunk/octave-forge/extra/control-devel/inst/@iddata/fft.m 2012-04-17 12:03:41 UTC (rev 10262) +++ trunk/octave-forge/extra/control-devel/inst/@iddata/fft.m 2012-04-17 13:33:51 UTC (rev 10263) @@ -33,6 +33,10 @@ if (nargin > 2) # no need to test nargin == 0, this is handled by built-in fft print_usage (); endif + + if (! dat.timedomain) + return; + endif [x, ~, ~, e] = size (dat); Added: trunk/octave-forge/extra/control-devel/inst/@iddata/ifft.m =================================================================== --- trunk/octave-forge/extra/control-devel/inst/@iddata/ifft.m (rev 0) +++ trunk/octave-forge/extra/control-devel/inst/@iddata/ifft.m 2012-04-17 13:33:51 UTC (rev 10263) @@ -0,0 +1,67 @@ +## Copyright (C) 2012 Lukas F. Reichlin +## +## This file is part of LTI Syncope. +## +## LTI Syncope is free software: you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation, either version 3 of the License, or +## (at your option) any later version. +## +## LTI Syncope is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with LTI Syncope. If not, see <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{dat} =} ifft (@var{dat}) +## @deftypefnx {Function File} {@var{dat} =} ifft (@var{dat}, @var{ord}) +## Detrend outputs and inputs of dataset @var{dat} by +## removing the best fit of a polynomial of order @var{ord}. +## If @var{ord} is not specified, default value 0 is taken. +## This corresponds to removing a constant. +## @end deftypefn + +## Author: Lukas Reichlin <luk...@gm...> +## Created: April 2012 +## Version: 0.1 + +function dat = ifft (dat) + + if (nargin > 1) # no need to test nargin == 0, this is handled by built-in fft + print_usage (); + endif + + if (dat.timedomain) + return; + endif + + [x, ~, ~, e] = size (dat); + + x = x(:); + n = num2cell (x); + %nconj = num2cell (x - ! rem (x, 2)); + nconj = num2cell (x - rem (x, 2)); + + + dat.y = cellfun (@(y, n, nconj) real (ifft ([y; conj(y(nconj:-1:2, :))], [], 1)) * sqrt (n+nconj), dat.y, n, nconj, "uniformoutput", false); + dat.u = cellfun (@(u, n, nconj) real (ifft ([u; conj(u(nconj:-1:2, :))], [], 1)) * sqrt (n+nconj), dat.u, n, nconj, "uniformoutput", false); + + ## ifft (x, n, dim=1) because x could be a row vector (n=1) + + % dat.w = cellfun (@(n, tsam) (0:fix(n/2)).' * (2*pi/abs(tsam)/n), n, dat.tsam, "uniformoutput", false); + ## abs(tsam) because of -1 for undefined sampling times + dat.timedomain = true; + +endfunction + + +%!shared DATD, Y, U +%! Y = 1:10; +%! U = 20:-2:1; +%! DAT = iddata (Y, U); +%! DATD = fft (DAT); +%!assert (DATD.y{1}, Y, 1e-10); +%!assert (DATD.u{1}, U, 1e-10); Modified: trunk/octave-forge/extra/control-devel/src/slitest.cc =================================================================== --- trunk/octave-forge/extra/control-devel/src/slitest.cc 2012-04-17 12:03:41 UTC (rev 10262) +++ trunk/octave-forge/extra/control-devel/src/slitest.cc 2012-04-17 13:33:51 UTC (rev 10263) @@ -134,6 +134,8 @@ int nsmp = y.rows (); // nsmp: number of samples // y.rows == u.rows is checked by iddata class // TODO: check minimal nsmp size + if (nsmp < 2*nobr) + error ("slitest: nsmp < 2*nobr"); int ldu = nsmp; int ldy = nsmp; This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2012-04-18 08:53:19
|
Revision: 10275 http://octave.svn.sourceforge.net/octave/?rev=10275&view=rev Author: paramaniac Date: 2012-04-18 08:53:07 +0000 (Wed, 18 Apr 2012) Log Message: ----------- control-devel: use larger ldwork as specified under further comments in slicot ib01ad Modified Paths: -------------- trunk/octave-forge/extra/control-devel/devel/ident.m trunk/octave-forge/extra/control-devel/src/slident.cc Modified: trunk/octave-forge/extra/control-devel/devel/ident.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/ident.m 2012-04-18 07:15:07 UTC (rev 10274) +++ trunk/octave-forge/extra/control-devel/devel/ident.m 2012-04-18 08:53:07 UTC (rev 10275) @@ -32,6 +32,7 @@ ctrl = 0; # confirm system order estimate n = 0; else # s & n non-empty + disp ("======== hallo ===============") nsmp = ns(1); nobr = fix ((nsmp+1)/(2*(m+l+1))); if (s > nobr) Modified: trunk/octave-forge/extra/control-devel/src/slident.cc =================================================================== --- trunk/octave-forge/extra/control-devel/src/slident.cc 2012-04-18 07:15:07 UTC (rev 10274) +++ trunk/octave-forge/extra/control-devel/src/slident.cc 2012-04-18 08:53:07 UTC (rev 10275) @@ -254,6 +254,7 @@ // TODO: Handle 'k' for DWORK int ldwork; + int ns = nsmp - 2*nobr + 1; if (alg == 'C') { @@ -289,7 +290,7 @@ } else // (alg == 'Q') { - int ns = nsmp - 2*nobr + 1; + // int ns = nsmp - 2*nobr + 1; if (ldr >= ns && batch == 'F') { @@ -311,14 +312,26 @@ ldwork = 6*(m+l)*nobr; } } -//////////////////////////////////////////////////////////////////// -// TO BE REMOVED !!! -//////////////////////////////////////////////////////////////////// -ldwork *= 2; -//ldwork = 1000000; -//////////////////////////////////////////////////////////////////// /* +IB01AD.f Lines 438-445 +C FURTHER COMMENTS +C +C For ALG = 'Q', BATCH = 'O' and LDR < NS, or BATCH <> 'O', the +C calculations could be rather inefficient if only minimal workspace +C (see argument LDWORK) is provided. It is advisable to provide as +C much workspace as possible. Almost optimal efficiency can be +C obtained for LDWORK = (NS+2)*(2*(M+L)*NOBR), assuming that the +C cache size is large enough to accommodate R, U, Y, and DWORK. +*/ + +// warning ("==================== ldwork before: %d =====================", ldwork); +// ldwork = (ns+2)*(2*(m+l)*nobr); +ldwork = max (ldwork, (ns+2)*(2*(m+l)*nobr)); +// warning ("==================== ldwork after: %d =====================", ldwork); + + +/* IB01AD.f Lines 291-195: c the workspace used for alg = 'q' is c ldrwrk*2*(m+l)*nobr + 4*(m+l)*nobr, This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2012-04-28 08:57:02
|
Revision: 10333 http://octave.svn.sourceforge.net/octave/?rev=10333&view=rev Author: paramaniac Date: 2012-04-28 08:56:56 +0000 (Sat, 28 Apr 2012) Log Message: ----------- control-devel: minor changes for debugging Modified Paths: -------------- trunk/octave-forge/extra/control-devel/devel/Destillation.m trunk/octave-forge/extra/control-devel/devel/ident.m trunk/octave-forge/extra/control-devel/src/slident.cc Modified: trunk/octave-forge/extra/control-devel/devel/Destillation.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/Destillation.m 2012-04-27 14:16:16 UTC (rev 10332) +++ trunk/octave-forge/extra/control-devel/devel/Destillation.m 2012-04-28 08:56:56 UTC (rev 10333) @@ -50,6 +50,8 @@ Y_dest_n30=Y(:,10:12); %} +clear all, close all, clc + load destill.dat U=destill(:,1:20); Y=destill(:,21:32); Modified: trunk/octave-forge/extra/control-devel/devel/ident.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/ident.m 2012-04-27 14:16:16 UTC (rev 10332) +++ trunk/octave-forge/extra/control-devel/devel/ident.m 2012-04-28 08:56:56 UTC (rev 10333) @@ -3,8 +3,8 @@ %nobr = 15; - meth = 2; - alg = 0; + meth = 2; % 2 % geht: meth/alg 1/1, + alg = 0; % 0 % geht nicht: meth/alg 0/1 jobd = 1; batch = 3; conct = 1; Modified: trunk/octave-forge/extra/control-devel/src/slident.cc =================================================================== --- trunk/octave-forge/extra/control-devel/src/slident.cc 2012-04-27 14:16:16 UTC (rev 10332) +++ trunk/octave-forge/extra/control-devel/src/slident.cc 2012-04-28 08:56:56 UTC (rev 10333) @@ -328,6 +328,7 @@ // warning ("==================== ldwork before: %d =====================", ldwork); // ldwork = (ns+2)*(2*(m+l)*nobr); ldwork = max (ldwork, (ns+2)*(2*(m+l)*nobr)); +// ldwork *= 3; // warning ("==================== ldwork after: %d =====================", ldwork); @@ -413,7 +414,10 @@ if (nuser > 0) { if (nuser < nobr) + { n = nuser; + // warning ("ident: nuser (%d) < nobr (%d), n = nuser", nuser, nobr); + } else error ("ident: 'nuser' invalid"); } @@ -430,6 +434,9 @@ int nsmpl = nsmp; + if (nsmpl < 2*(m+l)*nobr) + error ("slident: nsmpl (%d) < 2*(m+l)*nobr (%d)", nsmpl, nobr); + // arguments out int lda = max (1, n); int ldc = max (1, l); @@ -463,7 +470,7 @@ int ldw1; int ldw2; int ldw3; - +/* if (meth == 'M') { int ldw1a = max (2*(l*nobr-l)*n+2*n, (l*nobr-l)*n+n*n+7*n); @@ -507,9 +514,57 @@ ldw2 = l*nobr*n+m*nobr*(n+l)*(m*(n+l)+1)+ max ((n+l)*(n+l), 4*m*(n+l)+1); } +*/ + + int ldw1ax = max (2*(l*nobr-l)*n+2*n, (l*nobr-l)*n+n*n+7*n); + int ldw1bx = max (2*(l*nobr-l)*n+n*n+7*n, + (l*nobr-l)*n+n+6*m*nobr, + (l*nobr-l)*n+n+max (l+m*nobr, l*nobr + max (3*l*nobr+1, m))); + int ldw1x = max (ldw1ax, ldw1bx); + int aw; + + if (m == 0 || job == 'C') + aw = n + n*n; + else + aw = 0; + + int ldw2x = l*nobr*n + max ((l*nobr-l)*n+aw+2*n+max(5*n,(2*m+l)*nobr+l), 4*(m*nobr+n)+1, m*nobr+2*n+l ); + + + + int ldw1y = l*nobr*n + max ((l*nobr-l)*n+2*n+(2*m+l)*nobr+l, + 2*(l*nobr-l)*n+n*n+8*n, + n+4*(m*nobr+n)+1, + m*nobr+3*n+l); + int ldw2y; + if (m == 0 || job == 'C') + int ldw2y = 0; + else + int ldw2y = l*nobr*n+m*nobr*(n+l)*(m*(n+l)+1)+ max ((n+l)*(n+l), 4*m*(n+l)+1); + + + int ldw1az = max (2*(l*nobr-l)*n+2*n, (l*nobr-l)*n+n*n+7*n); + int ldw1bz = l*nobr*n + max ((l*nobr-l)*n+2*n+(2*m+l)*nobr+l, + 2*(l*nobr-l)*n+n*n+8*n, + n+4*(m*nobr+n)+1, + m*nobr+3*n+l); + + int ldw1z = max (ldw1az, ldw1bz); + + int ldw2z = l*nobr*n+m*nobr*(n+l)*(m*(n+l)+1)+ max ((n+l)*(n+l), 4*m*(n+l)+1); + + + ldw1 = max (ldw1x, ldw1y, ldw1z); + ldw2 = max (ldw2x, ldw2y, ldw2z); + + + ldw3 = max(4*n*n + 2*n*l + l*l + max (3*l, n*l), 14*n*n + 12*n + 5); ldwork_b = max (ldw1, ldw2, ldw3); + + // + ldwork_b *= 3; OCTAVE_LOCAL_BUFFER (int, iwork_b, liwork_b); OCTAVE_LOCAL_BUFFER (double, dwork_b, ldwork_b); @@ -690,6 +745,7 @@ retval(7) = k; retval(8) = x0; + //retval(8) = ColumnVector (n); //retval(0) = octave_value (n); //retval(1) = r; //retval(2) = sv; This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2012-05-05 13:04:06
|
Revision: 10365 http://octave.svn.sourceforge.net/octave/?rev=10365&view=rev Author: paramaniac Date: 2012-05-05 13:03:58 +0000 (Sat, 05 May 2012) Log Message: ----------- control-devel: add some debug code Modified Paths: -------------- trunk/octave-forge/extra/control-devel/devel/GlassFurnace.m trunk/octave-forge/extra/control-devel/devel/PowerPlant.m trunk/octave-forge/extra/control-devel/src/devel_slicot_functions.cc Added Paths: ----------- trunk/octave-forge/extra/control-devel/devel/CDplayer_a.m trunk/octave-forge/extra/control-devel/devel/Destillation_a.m trunk/octave-forge/extra/control-devel/devel/GlassFurnace_a.m trunk/octave-forge/extra/control-devel/devel/PowerPlant_a.m trunk/octave-forge/extra/control-devel/devel/ident_a.m trunk/octave-forge/extra/control-devel/src/slident_a.cc trunk/octave-forge/extra/control-devel/src/slident_b.cc trunk/octave-forge/extra/control-devel/src/slident_c.cc Added: trunk/octave-forge/extra/control-devel/devel/CDplayer_a.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/CDplayer_a.m (rev 0) +++ trunk/octave-forge/extra/control-devel/devel/CDplayer_a.m 2012-05-05 13:03:58 UTC (rev 10365) @@ -0,0 +1,67 @@ +%{ +Contributed by: + Favoreel + KULeuven + Departement Electrotechniek ESAT/SISTA +Kardinaal Mercierlaan 94 +B-3001 Leuven +Belgium + wou...@es... +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 + +%} + +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_a (dat, 15, 8) % s=15, n=8 + + +[y, t] = lsim (sys, U, [], x0); + +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 + + Added: trunk/octave-forge/extra/control-devel/devel/Destillation_a.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/Destillation_a.m (rev 0) +++ trunk/octave-forge/extra/control-devel/devel/Destillation_a.m 2012-05-05 13:03:58 UTC (rev 10365) @@ -0,0 +1,83 @@ +%{ +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); +%} + + close all, clc + +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); + + +dat = iddata (Y_dest, U_dest) + +[sys, x0] = ident_a (dat, 5, 4) % s=5, n=4 + + +[y, t] = lsim (sys, U_dest, [], x0); +%[y, t] = lsim (sys, U_dest); + +err = norm (Y_dest - y, 1) / norm (Y_dest, 1) + +figure (1) +%plot (t, Y_dest, 'b') +plot (t, Y_dest, 'b', t, y, 'r') +legend ('y measured', 'y simulated', 'location', 'southeast') + + Modified: trunk/octave-forge/extra/control-devel/devel/GlassFurnace.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/GlassFurnace.m 2012-05-05 06:48:41 UTC (rev 10364) +++ trunk/octave-forge/extra/control-devel/devel/GlassFurnace.m 2012-05-05 13:03:58 UTC (rev 10365) @@ -37,8 +37,9 @@ %} -clear all, close all, clc +close all, clc + load glassfurnace.dat T=glassfurnace(:,1); U=glassfurnace(:,2:4); Added: trunk/octave-forge/extra/control-devel/devel/GlassFurnace_a.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/GlassFurnace_a.m (rev 0) +++ trunk/octave-forge/extra/control-devel/devel/GlassFurnace_a.m 2012-05-05 13:03:58 UTC (rev 10365) @@ -0,0 +1,66 @@ +%{ +This file describes the data in the glassfurnace.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 glassfurnace (Philips) +3. Sampling time + +4. Number of samples: + 1247 samples +5. Inputs: + a. heating input + b. cooling input + c. heating input +6. Outputs: + a. 6 outputs from temperature sensors in a cross section of the + furnace +7. References: + a. Van Overschee P., De Moor B., N4SID : Subspace Algorithms for + the Identification of Combined Deterministic-Stochastic Systems, + Automatica, Special Issue on Statistical Signal Processing and Control, + Vol. 30, No. 1, 1994, pp. 75-93 + b. Van Overschee P., "Subspace identification : Theory, + Implementation, Application" , Ph.D. Thesis, K.U.Leuven, February 1995. +8. Known properties/peculiarities + +9. Some MATLAB-code to retrieve the data + !gunzip glassfurnace.dat.Z + load glassfurnace.dat + T=glassfurnace(:,1); + U=glassfurnace(:,2:4); + Y=glassfurnace(:,5:10); + +%} + + close all, clc + +load glassfurnace.dat +T=glassfurnace(:,1); +U=glassfurnace(:,2:4); +Y=glassfurnace(:,5:10); + + +dat = iddata (Y, U) + +[sys, x0] = ident_a (dat, 10, 5) % s=10, n=5 + + +[y, t] = lsim (sys, U, [], x0); + +err = norm (Y - y, 1) / norm (Y, 1) + +figure (1) +p = columns (Y); +for k = 1 : p + subplot (3, 2, k) + plot (t, Y(:,k), 'b', t, y(:,k), 'r') +endfor +%title ('DaISy: Glass Furnace') +%legend ('y measured', 'y simulated', 'location', 'southeast') + + Modified: trunk/octave-forge/extra/control-devel/devel/PowerPlant.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/PowerPlant.m 2012-05-05 06:48:41 UTC (rev 10364) +++ trunk/octave-forge/extra/control-devel/devel/PowerPlant.m 2012-05-05 13:03:58 UTC (rev 10365) @@ -40,7 +40,7 @@ %} -clear all, close all, clc +close all, clc load powerplant.dat U=powerplant(:,1:5); Added: trunk/octave-forge/extra/control-devel/devel/PowerPlant_a.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/PowerPlant_a.m (rev 0) +++ trunk/octave-forge/extra/control-devel/devel/PowerPlant_a.m 2012-05-05 13:03:58 UTC (rev 10365) @@ -0,0 +1,81 @@ +%{ +This file describes the data in the powerplant.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 power plant (Pont-sur-Sambre (France)) of 120 MW +3. Sampling time + 1228.8 sec +4. Number of samples: + 200 samples +5. Inputs: + 1. gas flow + 2. turbine valves opening + 3. super heater spray flow + 4. gas dampers + 5. air flow +6. Outputs: + 1. steam pressure + 2. main stem temperature + 3. reheat steam temperature +7. References: + a. R.P. Guidorzi, P. Rossi, Identification of a power plant from normal + operating records. Automatic control theory and applications (Canada, + Vol 2, pp 63-67, sept 1974. + b. Moonen M., De Moor B., Vandenberghe L., Vandewalle J., On- and + off-line identification of linear state-space models, International + Journal of Control, Vol. 49, Jan. 1989, pp.219-232 +8. Known properties/peculiarities + +9. Some MATLAB-code to retrieve the data + !gunzip powerplant.dat.Z + load powerplant.dat + U=powerplant(:,1:5); + Y=powerplant(:,6:8); + Yr=powerplant(:,9:11); + +%} + + close all, clc + +load powerplant.dat +U=powerplant(:,1:5); +Y=powerplant(:,6:8); +Yr=powerplant(:,9:11); + +inname = {'gas flow', + 'turbine valves opening', + 'super heater spray flow', + 'gas dampers', + 'air flow'}; + +outname = {'steam pressure', + 'main steam temperature', + 'reheat steam temperature'}; + +tsam = 1228.8; + +dat = iddata (Y, U, tsam, 'outname', outname, 'inname', inname) + +[sys, x0] = ident_a (dat, 10, 8) % s=10, n=8 + + +[y, t] = lsim (sys, U, [], x0); + +err = norm (Y - y, 1) / norm (Y, 1) + +figure (1) +p = columns (Y); +for k = 1 : p + subplot (3, 1, k) + plot (t, Y(:,k), 'b', t, y(:,k), 'r') +endfor +%title ('DaISy: Power Plant') +%legend ('y measured', 'y simulated', 'location', 'southeast') + +st = isstable (sys) + Added: trunk/octave-forge/extra/control-devel/devel/ident_a.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/ident_a.m (rev 0) +++ trunk/octave-forge/extra/control-devel/devel/ident_a.m 2012-05-05 13:03:58 UTC (rev 10365) @@ -0,0 +1,61 @@ +function [sys, x0] = ident_a (dat, s = [], nn = []) + +n=nn; + + %nobr = 15; + meth = 2; % 2 % geht: meth/alg 1/1, + alg = 0; % 0 % geht nicht: meth/alg 0/1 + jobd = 1; + batch = 3; + conct = 1; + ctrl = 0; %1; + rcond = 0.0; + tol = -1.0; % 0; + + [ns, l, m, e] = size (dat); + + if (isempty (s) && isempty (n)) + nsmp = ns(1); + nobr = fix ((nsmp+1)/(2*(m+l+1))); + ctrl = 0; # confirm system order estimate + n = 0; + % nsmp >= 2*(m+l+1)*nobr - 1 + % nobr <= (nsmp+1)/(2*(m+l+1)) + elseif (isempty (s)) + s = min (2*n, n+10); + nsmp = ns(1); + nobr = fix ((nsmp+1)/(2*(m+l+1))); + nobr = min (nobr, s); + ctrl = 1; # no confirmation + elseif (isempty (n)) + nobr = s; + ctrl = 0; # confirm system order estimate + n = 0; + else # s & n non-empty + nsmp = ns(1); + nobr = fix ((nsmp+1)/(2*(m+l+1))); + if (s > nobr) + error ("ident: s > nobr"); + endif + nobr = s; + ctrl = 1; + ## TODO: specify n for IB01BD + endif + + %nsmp = ns(1) + %nobr = fix ((nsmp+1)/(2*(m+l+1))) + % nsmp >= 2*(m+l+1)*nobr - 1 + % nobr <= (nsmp+1)/(2*(m+l+1)) +%nobr = 10 + [r, sv, n] = slident_a (dat.y{1}, dat.u{1}, nobr, n, meth, alg, jobd, batch, conct, ctrl, rcond, tol); +n +n = nn; + [a, b, c, d, q, ry, s, k] = slident_b (dat.y{1}, dat.u{1}, nobr, n, meth, alg, jobd, batch, conct, ctrl, rcond, tol, \ + r, sv, n); + + x0 = slident_c (dat.y{1}, dat.u{1}, nobr, n, meth, alg, jobd, batch, conct, ctrl, rcond, tol, \ + a, b, c, d); + + sys = ss (a, b, c, d, -1); + +endfunction Modified: trunk/octave-forge/extra/control-devel/src/devel_slicot_functions.cc =================================================================== --- trunk/octave-forge/extra/control-devel/src/devel_slicot_functions.cc 2012-05-05 06:48:41 UTC (rev 10364) +++ trunk/octave-forge/extra/control-devel/src/devel_slicot_functions.cc 2012-05-05 13:03:58 UTC (rev 10365) @@ -1,3 +1,6 @@ // #include "slib01ad.cc" // preprocess the input-output data #include "slident.cc" // system identification #include "slitest.cc" // debug system identification +#include "slident_a.cc" +#include "slident_b.cc" +#include "slident_c.cc" \ No newline at end of file Added: trunk/octave-forge/extra/control-devel/src/slident_a.cc =================================================================== --- trunk/octave-forge/extra/control-devel/src/slident_a.cc (rev 0) +++ trunk/octave-forge/extra/control-devel/src/slident_a.cc 2012-05-05 13:03:58 UTC (rev 10365) @@ -0,0 +1,422 @@ +/* + +Copyright (C) 2012 Lukas F. Reichlin + +This file is part of LTI Syncope. + +LTI Syncope is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +LTI Syncope is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with LTI Syncope. If not, see <http://www.gnu.org/licenses/>. + +SLICOT system identification +Uses SLICOT IB01AD, IB01BD and IB01CD by courtesy of NICONET e.V. +<http://www.slicot.org> + +Author: Lukas Reichlin <luk...@gm...> +Created: March 2012 +Version: 0.1 + +*/ + +#include <octave/oct.h> +#include <f77-fcn.h> +#include "common.h" + +extern "C" +{ + int F77_FUNC (ib01ad, IB01AD) + (char& METH, char& ALG, char& JOBD, + char& BATCH, char& CONCT, char& CTRL, + int& NOBR, int& M, int& L, + int& NSMP, + double* U, int& LDU, + double* Y, int& LDY, + int& N, + double* R, int& LDR, + double* SV, + double& RCOND, double& TOL, + int* IWORK, + double* DWORK, int& LDWORK, + int& IWARN, int& INFO); + + int F77_FUNC (ib01bd, IB01BD) + (char& METH, char& JOB, char& JOBCK, + int& NOBR, int& N, int& M, int& L, + int& NSMPL, + double* R, int& LDR, + double* A, int& LDA, + double* C, int& LDC, + double* B, int& LDB, + double* D, int& LDD, + double* Q, int& LDQ, + double* RY, int& LDRY, + double* S, int& LDS, + double* K, int& LDK, + double& TOL, + int* IWORK, + double* DWORK, int& LDWORK, + bool* BWORK, + int& IWARN, int& INFO); + + int F77_FUNC (ib01cd, IB01CD) + (char& JOBX0, char& COMUSE, char& JOB, + int& N, int& M, int& L, + int& NSMP, + double* A, int& LDA, + double* B, int& LDB, + double* C, int& LDC, + double* D, int& LDD, + double* U, int& LDU, + double* Y, int& LDY, + double* X0, + double* V, int& LDV, + double& TOL, + int* IWORK, + double* DWORK, int& LDWORK, + int& IWARN, int& INFO); +} + +// PKG_ADD: autoload ("slident_a", "devel_slicot_functions.oct"); +DEFUN_DLD (slident_a, args, nargout, + "-*- texinfo -*-\n\ +Slicot IB01AD Release 5.0\n\ +No argument checking.\n\ +For internal use only.") +{ + int nargin = args.length (); + octave_value_list retval; + + if (nargin != 12) + { + print_usage (); + } + else + { +//////////////////////////////////////////////////////////////////////////////////// +// SLICOT IB01AD - preprocess the input-output data // +//////////////////////////////////////////////////////////////////////////////////// + + // arguments in + char meth; + char alg; + char jobd; + char batch; + char conct; + char ctrl; + char metha; + char jobda; // ??? unused + + Matrix y = args(0).matrix_value (); + Matrix u = args(1).matrix_value (); + int nobr = args(2).int_value (); + int nuser = args(3).int_value (); + + const int imeth = args(4).int_value (); + const int ialg = args(5).int_value (); + const int ijobd = args(6).int_value (); + const int ibatch = args(7).int_value (); + const int iconct = args(8).int_value (); + const int ictrl = args(9).int_value (); + + double rcond = args(10).double_value (); + double tol = args(11).double_value (); + double tolb = args(10).double_value (); // tolb = rcond + + + switch (imeth) + { + case 0: + meth = 'M'; + metha = 'M'; + break; + case 1: + meth = 'N'; + metha = 'N'; + break; + case 2: + meth = 'C'; + metha = 'N'; // no typo here + break; + default: + error ("slib01ad: argument 'meth' invalid"); + } + + switch (ialg) + { + case 0: + alg = 'C'; + break; + case 1: + alg = 'F'; + break; + case 2: + alg = 'Q'; + break; + default: + error ("slib01ad: argument 'alg' invalid"); + } + + if (meth == 'C') + jobd = 'N'; + else if (ijobd == 0) + jobd = 'M'; + else + jobd = 'N'; + + switch (ibatch) + { + case 0: + batch = 'F'; + break; + case 1: + batch = 'I'; + break; + case 2: + batch = 'L'; + break; + case 3: + batch = 'O'; + break; + default: + error ("slib01ad: argument 'batch' invalid"); + } + + if (iconct == 0) + conct = 'C'; + else + conct = 'N'; + + if (ictrl == 0) + ctrl = 'C'; + else + ctrl = 'N'; + + + int m = u.columns (); // m: number of inputs + int l = y.columns (); // l: number of outputs + int nsmp = y.rows (); // nsmp: number of samples + // y.rows == u.rows is checked by iddata class + // TODO: check minimal nsmp size + + if (batch == 'O') + { + if (nsmp < 2*(m+l+1)*nobr - 1) + error ("slident: require NSMP >= 2*(M+L+1)*NOBR - 1"); + } + else + { + if (nsmp < 2*nobr) + error ("slident: require NSMP >= 2*NOBR"); + } + + int ldu; + + if (m == 0) + ldu = 1; + else // m > 0 + ldu = nsmp; + + int ldy = nsmp; + + // arguments out + int n; + int ldr; + + if (metha == 'M' && jobd == 'M') + ldr = max (2*(m+l)*nobr, 3*m*nobr); + else if (metha == 'N' || (metha == 'M' && jobd == 'N')) + ldr = 2*(m+l)*nobr; + else + error ("slib01ad: could not handle 'ldr' case"); + + Matrix r (ldr, 2*(m+l)*nobr); + ColumnVector sv (l*nobr); + + // workspace + int liwork; + + if (metha == 'N') // if METH = 'N' + liwork = (m+l)*nobr; + else if (alg == 'F') // if METH = 'M' and ALG = 'F' + liwork = m+l; + else // if METH = 'M' and ALG = 'C' or 'Q' + liwork = 0; + + // TODO: Handle 'k' for DWORK + + int ldwork; + int ns = nsmp - 2*nobr + 1; + + if (alg == 'C') + { + if (batch == 'F' || batch == 'I') + { + if (conct == 'C') + ldwork = (4*nobr-2)*(m+l); + else // (conct == 'N') + ldwork = 1; + } + else if (metha == 'M') // && (batch == 'L' || batch == 'O') + { + if (conct == 'C' && batch == 'L') + ldwork = max ((4*nobr-2)*(m+l), 5*l*nobr); + else if (jobd == 'M') + ldwork = max ((2*m-1)*nobr, (m+l)*nobr, 5*l*nobr); + else // (jobd == 'N') + ldwork = 5*l*nobr; + } + else // meth == 'N' && (batch == 'L' || batch == 'O') + { + ldwork = 5*(m+l)*nobr + 1; + } + } + else if (alg == 'F') + { + if (batch != 'O' && conct == 'C') + ldwork = (m+l)*2*nobr*(m+l+3); + else if (batch == 'F' || batch == 'I') // && conct == 'N' + ldwork = (m+l)*2*nobr*(m+l+1); + else // (batch == 'L' || '0' && conct == 'N') + ldwork = (m+l)*4*nobr*(m+l+1)+(m+l)*2*nobr; + } + else // (alg == 'Q') + { + // int ns = nsmp - 2*nobr + 1; + + if (ldr >= ns && batch == 'F') + { + ldwork = 4*(m+l)*nobr; + } + else if (ldr >= ns && batch == 'O') + { + if (metha == 'M') + ldwork = max (4*(m+l)*nobr, 5*l*nobr); + else // (meth == 'N') + ldwork = 5*(m+l)*nobr + 1; + } + else if (conct == 'C' && (batch == 'I' || batch == 'L')) + { + ldwork = 4*(nobr+1)*(m+l)*nobr; + } + else // if ALG = 'Q', (BATCH = 'F' or 'O', and LDR < NS), or (BATCH = 'I' or 'L' and CONCT = 'N') + { + ldwork = 6*(m+l)*nobr; + } + } + +/* +IB01AD.f Lines 438-445 +C FURTHER COMMENTS +C +C For ALG = 'Q', BATCH = 'O' and LDR < NS, or BATCH <> 'O', the +C calculations could be rather inefficient if only minimal workspace +C (see argument LDWORK) is provided. It is advisable to provide as +C much workspace as possible. Almost optimal efficiency can be +C obtained for LDWORK = (NS+2)*(2*(M+L)*NOBR), assuming that the +C cache size is large enough to accommodate R, U, Y, and DWORK. +*/ + +// warning ("==================== ldwork before: %d =====================", ldwork); +// ldwork = (ns+2)*(2*(m+l)*nobr); +ldwork = max (ldwork, (ns+2)*(2*(m+l)*nobr)); +// ldwork *= 3; +// warning ("==================== ldwork after: %d =====================", ldwork); + + +/* +IB01AD.f Lines 291-195: +c the workspace used for alg = 'q' is +c ldrwrk*2*(m+l)*nobr + 4*(m+l)*nobr, +c where ldrwrk = ldwork/(2*(m+l)*nobr) - 2; recommended +c value ldrwrk = ns, assuming a large enough cache size. +c for good performance, ldwork should be larger. + +somehow ldrwrk and ldwork must have been mixed up here + +*/ + + + OCTAVE_LOCAL_BUFFER (int, iwork, liwork); + OCTAVE_LOCAL_BUFFER (double, dwork, ldwork); + + // error indicators + int iwarn = 0; + int info = 0; + + + // SLICOT routine IB01AD + F77_XFCN (ib01ad, IB01AD, + (metha, alg, jobd, + batch, conct, ctrl, + nobr, m, l, + nsmp, + u.fortran_vec (), ldu, + y.fortran_vec (), ldy, + n, + r.fortran_vec (), ldr, + sv.fortran_vec (), + rcond, tol, + iwork, + dwork, ldwork, + iwarn, info)); + + + if (f77_exception_encountered) + error ("ident: exception in SLICOT subroutine IB01AD"); + + static const char* err_msg[] = { + "0: OK", + "1: a fast algorithm was requested (ALG = 'C', or 'F') " + "in sequential data processing, but it failed; the " + "routine can be repeatedly called again using the " + "standard QR algorithm", + "2: the singular value decomposition (SVD) algorithm did " + "not converge"}; + + static const char* warn_msg[] = { + "0: OK", + "1: the number of 100 cycles in sequential data " + "processing has been exhausted without signaling " + "that the last block of data was get; the cycle " + "counter was reinitialized", + "2: a fast algorithm was requested (ALG = 'C' or 'F'), " + "but it failed, and the QR algorithm was then used " + "(non-sequential data processing)", + "3: all singular values were exactly zero, hence N = 0 " + "(both input and output were identically zero)", + "4: the least squares problems with coefficient matrix " + "U_f, used for computing the weighted oblique " + "projection (for METH = 'N'), have a rank-deficient " + "coefficient matrix", + "5: the least squares problem with coefficient matrix " + "r_1 [6], used for computing the weighted oblique " + "projection (for METH = 'N'), has a rank-deficient " + "coefficient matrix"}; + + + error_msg ("ident", info, 2, err_msg); + warning_msg ("ident", iwarn, 5, warn_msg); + + + // resize + int rs = 2*(m+l)*nobr; + r.resize (rs, rs); + + + // return values + retval(0) = r; + retval(1) = sv; + retval(2) = octave_value (n); + } + + return retval; +} Added: trunk/octave-forge/extra/control-devel/src/slident_b.cc =================================================================== --- trunk/octave-forge/extra/control-devel/src/slident_b.cc (rev 0) +++ trunk/octave-forge/extra/control-devel/src/slident_b.cc 2012-05-05 13:03:58 UTC (rev 10365) @@ -0,0 +1,669 @@ +/* + +Copyright (C) 2012 Lukas F. Reichlin + +This file is part of LTI Syncope. + +LTI Syncope is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +LTI Syncope is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with LTI Syncope. If not, see <http://www.gnu.org/licenses/>. + +SLICOT system identification +Uses SLICOT IB01AD, IB01BD and IB01CD by courtesy of NICONET e.V. +<http://www.slicot.org> + +Author: Lukas Reichlin <luk...@gm...> +Created: March 2012 +Version: 0.1 + +*/ + +#include <octave/oct.h> +#include <f77-fcn.h> +#include "common.h" + +extern "C" +{ + int F77_FUNC (ib01ad, IB01AD) + (char& METH, char& ALG, char& JOBD, + char& BATCH, char& CONCT, char& CTRL, + int& NOBR, int& M, int& L, + int& NSMP, + double* U, int& LDU, + double* Y, int& LDY, + int& N, + double* R, int& LDR, + double* SV, + double& RCOND, double& TOL, + int* IWORK, + double* DWORK, int& LDWORK, + int& IWARN, int& INFO); + + int F77_FUNC (ib01bd, IB01BD) + (char& METH, char& JOB, char& JOBCK, + int& NOBR, int& N, int& M, int& L, + int& NSMPL, + double* R, int& LDR, + double* A, int& LDA, + double* C, int& LDC, + double* B, int& LDB, + double* D, int& LDD, + double* Q, int& LDQ, + double* RY, int& LDRY, + double* S, int& LDS, + double* K, int& LDK, + double& TOL, + int* IWORK, + double* DWORK, int& LDWORK, + bool* BWORK, + int& IWARN, int& INFO); + + int F77_FUNC (ib01cd, IB01CD) + (char& JOBX0, char& COMUSE, char& JOB, + int& N, int& M, int& L, + int& NSMP, + double* A, int& LDA, + double* B, int& LDB, + double* C, int& LDC, + double* D, int& LDD, + double* U, int& LDU, + double* Y, int& LDY, + double* X0, + double* V, int& LDV, + double& TOL, + int* IWORK, + double* DWORK, int& LDWORK, + int& IWARN, int& INFO); +} + +// PKG_ADD: autoload ("slident_b", "devel_slicot_functions.oct"); +DEFUN_DLD (slident_b, args, nargout, + "-*- texinfo -*-\n\ +Slicot IB01AD Release 5.0\n\ +No argument checking.\n\ +For internal use only.") +{ + int nargin = args.length (); + octave_value_list retval; + + if (nargin != 15) + { + print_usage (); + } + else + { +//////////////////////////////////////////////////////////////////////////////////// +// SLICOT IB01AD - preprocess the input-output data // +//////////////////////////////////////////////////////////////////////////////////// + + // arguments in + char meth; + char alg; + char jobd; + char batch; + char conct; + char ctrl; + char metha; + char jobda; // ??? unused + + Matrix y = args(0).matrix_value (); + Matrix u = args(1).matrix_value (); + int nobr = args(2).int_value (); + int nuser = args(3).int_value (); + + const int imeth = args(4).int_value (); + const int ialg = args(5).int_value (); + const int ijobd = args(6).int_value (); + const int ibatch = args(7).int_value (); + const int iconct = args(8).int_value (); + const int ictrl = args(9).int_value (); + + double rcond = args(10).double_value (); + double tol = args(11).double_value (); + double tolb = args(10).double_value (); // tolb = rcond + + Matrix r = args(12).matrix_value (); + Matrix sv = args(13).matrix_value (); + int n = args(14).int_value (); + + + switch (imeth) + { + case 0: + meth = 'M'; + metha = 'M'; + break; + case 1: + meth = 'N'; + metha = 'N'; + break; + case 2: + meth = 'C'; + metha = 'N'; // no typo here + break; + default: + error ("slib01ad: argument 'meth' invalid"); + } + + switch (ialg) + { + case 0: + alg = 'C'; + break; + case 1: + alg = 'F'; + break; + case 2: + alg = 'Q'; + break; + default: + error ("slib01ad: argument 'alg' invalid"); + } + + if (meth == 'C') + jobd = 'N'; + else if (ijobd == 0) + jobd = 'M'; + else + jobd = 'N'; + + switch (ibatch) + { + case 0: + batch = 'F'; + break; + case 1: + batch = 'I'; + break; + case 2: + batch = 'L'; + break; + case 3: + batch = 'O'; + break; + default: + error ("slib01ad: argument 'batch' invalid"); + } + + if (iconct == 0) + conct = 'C'; + else + conct = 'N'; + + if (ictrl == 0) + ctrl = 'C'; + else + ctrl = 'N'; + + + int m = u.columns (); // m: number of inputs + int l = y.columns (); // l: number of outputs + int nsmp = y.rows (); // nsmp: number of samples + // y.rows == u.rows is checked by iddata class + // TODO: check minimal nsmp size + + if (batch == 'O') + { + if (nsmp < 2*(m+l+1)*nobr - 1) + error ("slident: require NSMP >= 2*(M+L+1)*NOBR - 1"); + } + else + { + if (nsmp < 2*nobr) + error ("slident: require NSMP >= 2*NOBR"); + } + + int ldu; + + if (m == 0) + ldu = 1; + else // m > 0 + ldu = nsmp; + + int ldy = nsmp; + + // arguments out + //int n; + int ldr; + + if (metha == 'M' && jobd == 'M') + ldr = max (2*(m+l)*nobr, 3*m*nobr); + else if (metha == 'N' || (metha == 'M' && jobd == 'N')) + ldr = 2*(m+l)*nobr; + else + error ("slib01ad: could not handle 'ldr' case"); + + //Matrix r (ldr, 2*(m+l)*nobr); + //ColumnVector sv (l*nobr); + + // workspace + int liwork; + + if (metha == 'N') // if METH = 'N' + liwork = (m+l)*nobr; + else if (alg == 'F') // if METH = 'M' and ALG = 'F' + liwork = m+l; + else // if METH = 'M' and ALG = 'C' or 'Q' + liwork = 0; + + // TODO: Handle 'k' for DWORK + + int ldwork; + int ns = nsmp - 2*nobr + 1; + + if (alg == 'C') + { + if (batch == 'F' || batch == 'I') + { + if (conct == 'C') + ldwork = (4*nobr-2)*(m+l); + else // (conct == 'N') + ldwork = 1; + } + else if (metha == 'M') // && (batch == 'L' || batch == 'O') + { + if (conct == 'C' && batch == 'L') + ldwork = max ((4*nobr-2)*(m+l), 5*l*nobr); + else if (jobd == 'M') + ldwork = max ((2*m-1)*nobr, (m+l)*nobr, 5*l*nobr); + else // (jobd == 'N') + ldwork = 5*l*nobr; + } + else // meth == 'N' && (batch == 'L' || batch == 'O') + { + ldwork = 5*(m+l)*nobr + 1; + } + } + else if (alg == 'F') + { + if (batch != 'O' && conct == 'C') + ldwork = (m+l)*2*nobr*(m+l+3); + else if (batch == 'F' || batch == 'I') // && conct == 'N' + ldwork = (m+l)*2*nobr*(m+l+1); + else // (batch == 'L' || '0' && conct == 'N') + ldwork = (m+l)*4*nobr*(m+l+1)+(m+l)*2*nobr; + } + else // (alg == 'Q') + { + // int ns = nsmp - 2*nobr + 1; + + if (ldr >= ns && batch == 'F') + { + ldwork = 4*(m+l)*nobr; + } + else if (ldr >= ns && batch == 'O') + { + if (metha == 'M') + ldwork = max (4*(m+l)*nobr, 5*l*nobr); + else // (meth == 'N') + ldwork = 5*(m+l)*nobr + 1; + } + else if (conct == 'C' && (batch == 'I' || batch == 'L')) + { + ldwork = 4*(nobr+1)*(m+l)*nobr; + } + else // if ALG = 'Q', (BATCH = 'F' or 'O', and LDR < NS), or (BATCH = 'I' or 'L' and CONCT = 'N') + { + ldwork = 6*(m+l)*nobr; + } + } + +/* +IB01AD.f Lines 438-445 +C FURTHER COMMENTS +C +C For ALG = 'Q', BATCH = 'O' and LDR < NS, or BATCH <> 'O', the +C calculations could be rather inefficient if only minimal workspace +C (see argument LDWORK) is provided. It is advisable to provide as +C much workspace as possible. Almost optimal efficiency can be +C obtained for LDWORK = (NS+2)*(2*(M+L)*NOBR), assuming that the +C cache size is large enough to accommodate R, U, Y, and DWORK. +*/ + +// warning ("==================== ldwork before: %d =====================", ldwork); +// ldwork = (ns+2)*(2*(m+l)*nobr); +ldwork = max (ldwork, (ns+2)*(2*(m+l)*nobr)); +// ldwork *= 3; +// warning ("==================== ldwork after: %d =====================", ldwork); + + +/* +IB01AD.f Lines 291-195: +c the workspace used for alg = 'q' is +c ldrwrk*2*(m+l)*nobr + 4*(m+l)*nobr, +c where ldrwrk = ldwork/(2*(m+l)*nobr) - 2; recommended +c value ldrwrk = ns, assuming a large enough cache size. +c for good performance, ldwork should be larger. + +somehow ldrwrk and ldwork must have been mixed up here + +*/ + +#if 0 + OCTAVE_LOCAL_BUFFER (int, iwork, liwork); + OCTAVE_LOCAL_BUFFER (double, dwork, ldwork); + + // error indicators + int iwarn = 0; + int info = 0; + + + // SLICOT routine IB01AD + F77_XFCN (ib01ad, IB01AD, + (metha, alg, jobd, + batch, conct, ctrl, + nobr, m, l, + nsmp, + u.fortran_vec (), ldu, + y.fortran_vec (), ldy, + n, + r.fortran_vec (), ldr, + sv.fortran_vec (), + rcond, tol, + iwork, + dwork, ldwork, + iwarn, info)); + + + if (f77_exception_encountered) + error ("ident: exception in SLICOT subroutine IB01AD"); + + static const char* err_msg[] = { + "0: OK", + "1: a fast algorithm was requested (ALG = 'C', or 'F') " + "in sequential data processing, but it failed; the " + "routine can be repeatedly called again using the " + "standard QR algorithm", + "2: the singular value decomposition (SVD) algorithm did " + "not converge"}; + + static const char* warn_msg[] = { + "0: OK", + "1: the number of 100 cycles in sequential data " + "processing has been exhausted without signaling " + "that the last block of data was get; the cycle " + "counter was reinitialized", + "2: a fast algorithm was requested (ALG = 'C' or 'F'), " + "but it failed, and the QR algorithm was then used " + "(non-sequential data processing)", + "3: all singular values were exactly zero, hence N = 0 " + "(both input and output were identically zero)", + "4: the least squares problems with coefficient matrix " + "U_f, used for computing the weighted oblique " + "projection (for METH = 'N'), have a rank-deficient " + "coefficient matrix", + "5: the least squares problem with coefficient matrix " + "r_1 [6], used for computing the weighted oblique " + "projection (for METH = 'N'), has a rank-deficient " + "coefficient matrix"}; + + + error_msg ("ident", info, 2, err_msg); + warning_msg ("ident", iwarn, 5, warn_msg); + + + // resize + int rs = 2*(m+l)*nobr; + r.resize (rs, rs); + + if (nuser > 0) + { + if (nuser < nobr) + { + n = nuser; + // warning ("ident: nuser (%d) < nobr (%d), n = nuser", nuser, nobr); + } + else + error ("ident: 'nuser' invalid"); + } +#endif +//////////////////////////////////////////////////////////////////////////////////// +// SLICOT IB01BD - estimating system matrices, Kalman gain, and covariances // +//////////////////////////////////////////////////////////////////////////////////// + + // arguments in + char job = 'A'; + char jobck = 'K'; + + // TODO: if meth == 'C', which meth should be taken for IB01AD.f, 'M' or 'N'? + + int nsmpl = nsmp; + + if (nsmpl < 2*(m+l)*nobr) + error ("slident: nsmpl (%d) < 2*(m+l)*nobr (%d)", nsmpl, nobr); + + // arguments out + int lda = max (1, n); + int ldc = max (1, l); + int ldb = max (1, n); + int ldd = max (1, l); + int ldq = n; // if JOBCK = 'C' or 'K' + int ldry = l; // if JOBCK = 'C' or 'K' + int lds = n; // if JOBCK = 'C' or 'K' + int ldk = n; // if JOBCK = 'K' + + Matrix a (lda, n); + Matrix c (ldc, n); + Matrix b (ldb, m); + Matrix d (ldd, m); + + Matrix q (ldq, n); + Matrix ry (ldry, l); + Matrix s (lds, l); + Matrix k (ldk, l); + + // workspace + int liwork_b; + int liw1; + int liw2; + + liw1 = max (n, m*nobr+n, l*nobr, m*(n+l)); + liw2 = n*n; // if JOBCK = 'K' + liwork_b = max (liw1, liw2); + + int ldwork_b; + int ldw1; + int ldw2; + int ldw3; +/* + if (meth == 'M') + { + int ldw1a = max (2*(l*nobr-l)*n+2*n, (l*nobr-l)*n+n*n+7*n); + int ldw1b = max (2*(l*nobr-l)*n+n*n+7*n, + (l*nobr-l)*n+n+6*m*nobr, + (l*nobr-l)*n+n+max (l+m*nobr, l*nobr + max (3*l*nobr+1, m))); + ldw1 = max (ldw1a, ldw1b); + + int aw; + + if (m == 0 || job == 'C') + aw = n + n*n; + else + aw = 0; + + ldw2 = l*nobr*n + max ((l*nobr-l)*n+aw+2*n+max(5*n,(2*m+l)*nobr+l), 4*(m*nobr+n)+1, m*nobr+2*n+l ); + } + else if (meth == 'N') + { + ldw1 = l*nobr*n + max ((l*nobr-l)*n+2*n+(2*m+l)*nobr+l, + 2*(l*nobr-l)*n+n*n+8*n, + n+4*(m*nobr+n)+1, + m*nobr+3*n+l); + + if (m == 0 || job == 'C') + ldw2 = 0; + else + ldw2 = l*nobr*n+m*nobr*(n+l)*(m*(n+l)+1)+ max ((n+l)*(n+l), 4*m*(n+l)+1); + + } + else // (meth == 'C') + { + int ldw1a = max (2*(l*nobr-l)*n+2*n, (l*nobr-l)*n+n*n+7*n); + int ldw1b = l*nobr*n + max ((l*nobr-l)*n+2*n+(2*m+l)*nobr+l, + 2*(l*nobr-l)*n+n*n+8*n, + n+4*(m*nobr+n)+1, + m*nobr+3*n+l); + + ldw1 = max (ldw1a, ldw1b); + + ldw2 = l*nobr*n+m*nobr*(n+l)*(m*(n+l)+1)+ max ((n+l)*(n+l), 4*m*(n+l)+1); + + } +*/ + + int ldw1ax = max (2*(l*nobr-l)*n+2*n, (l*nobr-l)*n+n*n+7*n); + int ldw1bx = max (2*(l*nobr-l)*n+n*n+7*n, + (l*nobr-l)*n+n+6*m*nobr, + (l*nobr-l)*n+n+max (l+m*nobr, l*nobr + max (3*l*nobr+1, m))); + int ldw1x = max (ldw1ax, ldw1bx); + + int aw; + + if (m == 0 || job == 'C') + aw = n + n*n; + else + aw = 0; + + int ldw2x = l*nobr*n + max ((l*nobr-l)*n+aw+2*n+max(5*n,(2*m+l)*nobr+l), 4*(m*nobr+n)+1, m*nobr+2*n+l ); + + + + int ldw1y = l*nobr*n + max ((l*nobr-l)*n+2*n+(2*m+l)*nobr+l, + 2*(l*nobr-l)*n+n*n+8*n, + n+4*(m*nobr+n)+1, + m*nobr+3*n+l); + int ldw2y; + if (m == 0 || job == 'C') + int ldw2y = 0; + else + int ldw2y = l*nobr*n+m*nobr*(n+l)*(m*(n+l)+1)+ max ((n+l)*(n+l), 4*m*(n+l)+1); + + + int ldw1az = max (2*(l*nobr-l)*n+2*n, (l*nobr-l)*n+n*n+7*n); + int ldw1bz = l*nobr*n + max ((l*nobr-l)*n+2*n+(2*m+l)*nobr+l, + 2*(l*nobr-l)*n+n*n+8*n, + n+4*(m*nobr+n)+1, + m*nobr+3*n+l); + + int ldw1z = max (ldw1az, ldw1bz); + + int ldw2z = l*nobr*n+m*nobr*(n+l)*(m*(n+l)+1)+ max ((n+l)*(n+l), 4*m*(n+l)+1); + + + ldw1 = max (ldw1x, ldw1y, ldw1z); + ldw2 = max (ldw2x, ldw2y, ldw2z); + + + + ldw3 = max(4*n*n + 2*n*l + l*l + max (3*l, n*l), 14*n*n + 12*n + 5); + ldwork_b = max (ldw1, ldw2, ldw3); + + // + ldwork_b *= 3; + + OCTAVE_LOCAL_BUFFER (int, iwork_b, liwork_b); + OCTAVE_LOCAL_BUFFER (double, dwork_b, ldwork_b); + OCTAVE_LOCAL_BUFFER (bool, bwork, 2*n); + + + // error indicators + int iwarn_b = 0; + int info_b = 0; + + + // SLICOT routine IB01BD + F77_XFCN (ib01bd, IB01BD, + (meth, job, jobck, + nobr, n, m, l, + nsmpl, + r.fortran_vec (), ldr, + a.fortran_vec (), lda, + c.fortran_vec (), ldc, + b.fortran_vec (), ldb, + d.fortran_vec (), ldd, + q.fortran_vec (), ldq, + ry.fortran_vec (), ldry, + s.fortran_vec (), lds, + k.fortran_vec (), ldk, + tolb, + iwork_b, + dwork_b, ldwork_b, + bwork, + iwarn_b, info_b)); + + + if (f77_exception_encountered) + error ("ident: exception in SLICOT subroutine IB01BD"); + + static const char* err_msg_b[] = { + "0: OK", + "1: error message not specified", + "2: the singular value decomposition (SVD) algorithm did " + "not converge", + "3: a singular upper triangular matrix was found", + "4: matrix A is (numerically) singular in discrete-" + "time case", + "5: the Hamiltonian or symplectic matrix H cannot be " + "reduced to real Schur form", + "6: the real Schur form of the Hamiltonian or " + "symplectic matrix H cannot be appropriately ordered", + "7: the Hamiltonian or symplectic matrix H has less " + "than N stable eigenvalues", + "8: the N-th order system of linear algebraic " + "equations, from which the solution matrix X would " + "be obtained, is singular to working precision", + "9: the QR algorithm failed to complete the reduction " + "of the matrix Ac to Schur canonical form, T", + "10: the QR algorithm did not converge"}; + + static const char* warn_msg_b[] = { + "0: OK", + "1: warning message not specified", + "2: warning message not specified", + "3: warning message not specified", + "4: a least squares problem to be solved has a " + "rank-deficient coefficient matrix", + "5: the computed covariance matrices are too small. " + "The problem seems to be a deterministic one; the " + "gain matrix is set to zero"}; + + + error_msg ("ident", info_b, 10, err_msg_b); + warning_msg ("ident", iwarn_b, 5, warn_msg_b); + + // resize + a.resize (n, n); + c.resize (l, n); + b.resize (n, m); + d.resize (l, m); + + q.resize (n, n); + ry.resize (l, l); + s.resize (n, l); + k.resize (n, l); + + + + // return values + retval(0) = a; + retval(1) = b; + retval(2) = c; + retval(3) = d; + + retval(4) = q; + retval(5) = ry; + retval(6) = s; + retval(7) = k; + } + + return retval; +} Added: trunk/octave-forge/extra/control-devel/src/slident_c.cc =================================================================== --- trunk/octave-forge/extra/control-devel/src/slident_c.cc (rev 0) +++ trunk/octave-forge/extra/control-devel/src/slident_c.cc 2012-05-05 13:03:58 UTC (rev 10365) @@ -0,0 +1,473 @@ +/* + +Copyright (C) 2012 Lukas F. Reichlin + +This file is part of LTI Syncope. + +LTI Syncope is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +LTI Syncope is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with LTI Syncope. If not, see <http://www.gnu.org/licenses/>. + +SLICOT system identification +Uses SLICOT IB01AD, IB01BD and IB01CD by courtesy of NICONET e.V. +<http://www.slicot.org> + +Author: Lukas Reichlin <luk...@gm...> +Created: March 2012 +Version: 0.1 + +*/ + +#include <octave/oct.h> +#include <f77-fcn.h> +#include "common.h" + +extern "C" +{ + int F77_FUNC (ib01ad, IB01AD) + (char& METH, char& ALG, char& JOBD, + char& BATCH, char& CONCT, char& CTRL, + int& NOBR, int& M, int& L, + int& NSMP, + double* U, int& LDU, + double* Y, int& LDY, + int& N, + double* R, int& LDR, + double* SV, + double& RCOND, double& TOL, + int* IWORK, + double* DWORK, int& LDWORK, + int& IWARN, int& INFO); + + int F77_FUNC (ib01bd, IB01BD) + (char& METH, char& JOB, char& JOBCK, + int& NOBR, int& N, int& M, int& L, + int& NSMPL, + double* R, int& LDR, + double* A, int& LDA, + double* C, int& LDC, + double* B, int& LDB, + double* D, int& LDD, + double* Q, int& LDQ, + double* RY, int& LDRY, + double* S, int& LDS, + double* K, int& LDK, + double& TOL, + int* IWORK, + double* DWORK, int& LDWORK, + bool* BWORK, + int& IWARN, int& INFO); + + int F77_FUNC (ib01cd, IB01CD) + (char& JOBX0, char& COMUSE, char& JOB, + int& N, int& M, int& L, + int& NSMP, + double* A, int& LDA, + double* B, int& LDB, + double* C, int& LDC, + double* D, int& LDD, + double* U, int& LDU, + double* Y, int& LDY, + double* X0, + double* V, int& LDV, + double& TOL, + int* IWORK, + double* DWORK, int& LDWORK, + int& IWARN, int& INFO); +} + +// PKG_ADD: autoload ("slident_c", "devel_slicot_functions.oct"); +DEFUN_DLD (slident_c, args, nargout, + "-*- texinfo -*-\n\ +Slicot IB01AD Release 5.0\n\ +No argument checking.\n\ +For internal use only.") +{ + int nargin = args.length (); + octave_value_list retval; + + if (nargin != 16) + { + print_usage (); + } + else + { +//////////////////////////////////////////////////////////////////////////////////// +// SLICOT IB01AD - preprocess the input-output data // +//////////////////////////////////////////////////////////////////////////////////// + + // arguments in + char meth; + char alg; + char jobd; + char batch; + char conct; + char ctrl; + char metha; + char jobda; // ??? unused + + Matrix y = args(0).matrix_value (); + Matrix u = args(1).matrix_value (); + int nobr = args(2).int_value (); + int nuser = args(3).int_value (); + + const int imeth = args(4).int_value (); + const int ialg = args(5).int_value (); + const int ijobd = args(6).int_value (); + const int ibatch = args(7).int_value (); + const int iconct = args(8).int_value (); + const int ictrl = args(9).int_value (); + + double rcond = args(10).double_value (); + double tol = args(11).double_value (); + double tolb = args(10).double_value (); // tolb = rcond + + Matrix a = args(12).matrix_value (); + Matrix b = args(13).matrix_value (); + Matrix c = args(14).matrix_value (); + Matrix d = args(15).matrix_value (); + + + switch (imeth) + { + case 0: + meth = 'M'; + metha = 'M'; + break; + case 1: + meth = 'N'; + metha = 'N'; + break; + case 2: + meth = 'C'; + metha = 'N'; // no typo here + break; + default: + error ("slib01ad: argument 'meth' invalid"); + } + + switch (ialg) + { + case 0: + alg = 'C'; + break; + case 1: + alg = 'F'; + break; + case 2: + alg = 'Q'; + break; + default: + error ("slib01ad: argument 'alg' invalid"); + } + + if (meth == 'C') + jobd = 'N'; + else if (ijobd == 0) + jobd = 'M'; + else + jobd = 'N'; + + switch (ibatch) + { + case 0: + batch = 'F'; + break; + case 1: + batch = 'I'; + break; + case 2: + batch = 'L'; + break; + case 3: + batch = 'O'; + break; + default: + error ("slib01ad: argument 'batch' invalid"); + } + + if (iconct == 0) + conct = 'C'; + else + conct = 'N'; + + if (ictrl == 0) + ctrl = 'C'; + else + ctrl = 'N'; + + + int m = u.columns (); // m: number of inputs + int l = y.columns (); // l: number of outputs + int nsmp = y.rows (); // nsmp: number of samples + // y.rows == u.rows is checked by iddata class + // TODO: check minimal nsmp size + + if (batch == 'O') + { + if (nsmp < 2*(m+l+1)*nobr - 1) + error ("slident: require NSMP >= 2*(M+L+1)*NOBR - 1"); + } + else + { + if (nsmp < 2*nobr) + error ("slident: require NSMP >= 2*NOBR"); + } + + int ldu; + + if (m == 0) + ldu = 1; + else // m > 0 + ldu = nsmp; + + int ldy = nsmp; + + // arguments out + int n; + int ldr; + + if (metha == 'M' && jobd == 'M') + ldr = max (2*(m+l)*nobr, 3*m*nobr); + else if (metha == 'N' || (metha == 'M' && jobd == 'N')) + ldr = 2*(m+l)*nobr; + else + error ("slib01ad: could not handle 'ldr' case"); + + Matrix r (ldr, 2*(m+l)*nobr); + ColumnVector sv (l*nobr); + + // workspace + int liwork; + + if (metha == 'N') // if METH = 'N' + liwork = (m+l)*nobr; + else if (alg == 'F') // if METH = 'M' and ALG = 'F' + liwork = m+l; + else // if METH = 'M' and ALG = 'C' or 'Q' + liwork = 0; + + // TODO: Handle 'k' for DWORK + + int ldwork; + int ns = nsmp - 2*nobr + 1; + + if (alg == 'C') + { + if (batch == 'F' || batch == 'I') + { + if (conct == 'C') + ldwork = (4*nobr-2)*(m+l); + else // (conct == 'N') + ldwork = 1; + } + else if (metha == 'M') // && (batch == 'L' || batch == 'O') + { + if (conct == 'C' && batch == 'L') + ldwork = max ((4*nobr-2)*(m+l), 5*l*nobr); + else if (jobd == 'M') + ldwork = max ((2*m-1)*nobr, (m+l)*nobr, 5*l*nobr); + else // (jobd == 'N') + ldwork = 5*l*nobr; + } + else // meth == 'N' && (batch == 'L' || batch == 'O') + { + ldwork = 5*(m+l)*nobr + 1; + } + } + else if (alg == 'F') + { + if (batch != 'O' && conct == 'C') + ldwork = (m+l)*2*nobr*(m+l+3); + else if (batch == 'F' || batch == 'I') // && conct == 'N' + ldwork = (m+l)*2*nobr*(m+l+1); + else // (batch == 'L' || '0' && conct == 'N') + ldwork = (m+l)*4*nobr*(m+l+1)+(m+l)*2*nobr; + } + else // (alg == 'Q') + { + // int ns = nsmp - 2*nobr + 1; + + if (ldr >= ns && batch == 'F') + { + ldwork = 4*(m+l)*nobr; + } + else if (ldr >= ns && batch == 'O') + { + if (metha == 'M') + ldwork = max (4*(m+l)*nobr, 5*l*nobr); + else // (meth == 'N') + ldwork = 5*(m+l)*nobr + 1; + } + else if (conct == 'C' && (batch == 'I' || batch == 'L')) + { + ldwork = 4*(nobr+1)*(m+l)*nobr; + } + else // if ALG = 'Q', (BATCH = 'F' or 'O', and LDR < NS), or (BATCH = 'I' or 'L' and CONCT = 'N') + { + ldwork = 6*(m+l)*nobr; + } + } + +/* +IB01AD.f Lines 438-445 +C FURTHER COMMENTS +C +C For ALG = 'Q', BATCH = 'O' and LDR < NS, or BATCH <> 'O', the +C calculations could be rather inefficient if only minimal workspace +C (see argument LDWORK) is provided. It is advisable to provide as +C much workspace as possible. Almost optimal efficiency can be +C obtained for LDWORK = (NS+2)*(2*(M+L)*NOBR), assuming tha... [truncated message content] |
From: <par...@us...> - 2012-05-07 13:38:37
|
Revision: 10375 http://octave.svn.sourceforge.net/octave/?rev=10375&view=rev Author: paramaniac Date: 2012-05-07 13:38:26 +0000 (Mon, 07 May 2012) Log Message: ----------- control-devel: fix miso support, add testcase Modified Paths: -------------- trunk/octave-forge/extra/control-devel/inst/arx.m Added Paths: ----------- trunk/octave-forge/extra/control-devel/devel/pHarx.m Added: trunk/octave-forge/extra/control-devel/devel/pHarx.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/pHarx.m (rev 0) +++ trunk/octave-forge/extra/control-devel/devel/pHarx.m 2012-05-07 13:38:26 UTC (rev 10375) @@ -0,0 +1,67 @@ +%{ +Contributed by: + Jairo Espinosa + K.U.Leuven ESAT-SISTA + K.Mercierlaan 94 + B3001 Heverlee + Jai...@es... + +Description: + Simulation data of a pH neutralization process in a constant volume + stirring tank. + Volume of the tank 1100 liters + Concentration of the acid solution (HAC) 0.0032 Mol/l + Concentration of the base solution (NaOH) 0,05 Mol/l +Sampling: + 10 sec +Number: + 2001 +Inputs: + u1: Acid solution flow in liters + u2: Base solution flow in liters + +Outputs: + y: pH of the solution in the tank + +References: + T.J. Mc Avoy, E.Hsu and S.Lowenthal, Dynamics of pH in controlled + stirred tank reactor, Ind.Eng.Chem.Process Des.Develop.11(1972) + 71-78 + +Properties: + Highly non-linear system. + +Columns: + Column 1: time-steps + Column 2: input u1 + Column 3: input u2 + Column 4: output y + +Category: + Process industry systems + +%} + +clear all, close all, clc + +load pHdata.dat +U=pHdata(:,2:3); +Y=pHdata(:,4); + + +dat = iddata (Y, U) + +% [sys, x0] = ident (dat, 15, 6) % s=15, n=6 +sys = arx (dat, 6, [3,3]) + +% [y, t] = lsim (sys, U, [], x0); +[y, t] = lsim (sys(:, 1:2), U); + + +err = norm (Y - y, 1) / norm (Y, 1) + +figure (1) +plot (t, Y(:,1), 'b', t, y(:,1), 'r') + + + Modified: trunk/octave-forge/extra/control-devel/inst/arx.m =================================================================== --- trunk/octave-forge/extra/control-devel/inst/arx.m 2012-05-07 13:22:01 UTC (rev 10374) +++ trunk/octave-forge/extra/control-devel/inst/arx.m 2012-05-07 13:38:26 UTC (rev 10375) @@ -12,7 +12,8 @@ error ("arx: "); endif - if (! is_real_scalar (na, nb)) +## if (! is_real_scalar (na, nb)) + if (! is_real_vector (na, nb)) error ("arx: "); ## Test for integers ## numel (nb) == size (dat, 3) @@ -33,7 +34,7 @@ ## 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, :) + Phi = Phi(n:end, :); ## Theta = Phi \ Y(n+1:end, :); # naive formula This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2012-05-09 10:36:26
|
Revision: 10387 http://octave.svn.sourceforge.net/octave/?rev=10387&view=rev Author: paramaniac Date: 2012-05-09 10:36:17 +0000 (Wed, 09 May 2012) Log Message: ----------- control-devel: fix arx argument checking for degrees of numerator and denominator Modified Paths: -------------- trunk/octave-forge/extra/control-devel/devel/GlassFurnaceARX.m trunk/octave-forge/extra/control-devel/inst/arx.m Added Paths: ----------- trunk/octave-forge/extra/control-devel/devel/CDplayerARX.m Added: trunk/octave-forge/extra/control-devel/devel/CDplayerARX.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/CDplayerARX.m (rev 0) +++ trunk/octave-forge/extra/control-devel/devel/CDplayerARX.m 2012-05-09 10:36:17 UTC (rev 10387) @@ -0,0 +1,78 @@ +%{ +Contributed by: + Favoreel + KULeuven + Departement Electrotechniek ESAT/SISTA +Kardinaal Mercierlaan 94 +B-3001 Leuven +Belgium + wou...@es... +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 + + Modified: trunk/octave-forge/extra/control-devel/devel/GlassFurnaceARX.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/GlassFurnaceARX.m 2012-05-09 03:49:09 UTC (rev 10386) +++ trunk/octave-forge/extra/control-devel/devel/GlassFurnaceARX.m 2012-05-09 10:36:17 UTC (rev 10387) @@ -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); Modified: trunk/octave-forge/extra/control-devel/inst/arx.m =================================================================== --- trunk/octave-forge/extra/control-devel/inst/arx.m 2012-05-09 03:49:09 UTC (rev 10386) +++ trunk/octave-forge/extra/control-devel/inst/arx.m 2012-05-09 10:36:17 UTC (rev 10387) @@ -9,27 +9,38 @@ endif if (! isa (dat, "iddata")) - error ("arx: "); + error ("arx: first argument must be an iddata dataset"); endif - -## if (! is_real_scalar (na, nb)) - if (! is_real_vector (na, nb)) - error ("arx: "); - ## Test for integers - ## numel (nb) == size (dat, 3) - 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); den = 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 This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2012-05-10 12:09:17
|
Revision: 10399 http://octave.svn.sourceforge.net/octave/?rev=10399&view=rev Author: paramaniac Date: 2012-05-10 12:09:11 +0000 (Thu, 10 May 2012) Log Message: ----------- control-devel: minor changes Modified Paths: -------------- trunk/octave-forge/extra/control-devel/devel/pH.m trunk/octave-forge/extra/control-devel/devel/pHarx.m trunk/octave-forge/extra/control-devel/inst/arx.m Modified: trunk/octave-forge/extra/control-devel/devel/pH.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/pH.m 2012-05-10 11:24:35 UTC (rev 10398) +++ trunk/octave-forge/extra/control-devel/devel/pH.m 2012-05-10 12:09:11 UTC (rev 10399) @@ -57,6 +57,7 @@ [y, t] = lsim (sys, U, [], x0); err = norm (Y - y, 1) / norm (Y, 1) +st = isstable (sys) figure (1) plot (t, Y(:,1), 'b', t, y(:,1), 'r') Modified: trunk/octave-forge/extra/control-devel/devel/pHarx.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/pHarx.m 2012-05-10 11:24:35 UTC (rev 10398) +++ trunk/octave-forge/extra/control-devel/devel/pHarx.m 2012-05-10 12:09:11 UTC (rev 10399) @@ -52,7 +52,7 @@ dat = iddata (Y, U) % [sys, x0] = ident (dat, 15, 6) % s=15, n=6 -sys = arx (dat, 6, [6,6]) % normally na = nb +sys = arx (dat, 6, 6) % normally na = nb % [y, t] = lsim (sys, U, [], x0); [y, t] = lsim (sys(:, 1:2), U); Modified: trunk/octave-forge/extra/control-devel/inst/arx.m =================================================================== --- trunk/octave-forge/extra/control-devel/inst/arx.m 2012-05-10 11:24:35 UTC (rev 10398) +++ trunk/octave-forge/extra/control-devel/inst/arx.m 2012-05-10 12:09:11 UTC (rev 10399) @@ -87,7 +87,7 @@ ## 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 + [U, S, V] = svd (phi{1}, 0); # 0 for "economy size" decomposition S = diag (S); # extract main diagonal r = sum (S > eps*S(1)); V = V(:, 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-12 15:33:00
|
Revision: 10427 http://octave.svn.sourceforge.net/octave/?rev=10427&view=rev Author: paramaniac Date: 2012-05-12 15:32:53 +0000 (Sat, 12 May 2012) Log Message: ----------- control-devel: minor changes Modified Paths: -------------- trunk/octave-forge/extra/control-devel/devel/GlassFurnace.m trunk/octave-forge/extra/control-devel/devel/PowerPlant.m trunk/octave-forge/extra/control-devel/src/slident.cc Modified: trunk/octave-forge/extra/control-devel/devel/GlassFurnace.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/GlassFurnace.m 2012-05-12 15:06:57 UTC (rev 10426) +++ trunk/octave-forge/extra/control-devel/devel/GlassFurnace.m 2012-05-12 15:32:53 UTC (rev 10427) @@ -38,7 +38,7 @@ %} -close all, clc +clear all, close all, clc load glassfurnace.dat T=glassfurnace(:,1); Modified: trunk/octave-forge/extra/control-devel/devel/PowerPlant.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/PowerPlant.m 2012-05-12 15:06:57 UTC (rev 10426) +++ trunk/octave-forge/extra/control-devel/devel/PowerPlant.m 2012-05-12 15:32:53 UTC (rev 10427) @@ -40,7 +40,7 @@ %} -close all, clc +clear all, close all, clc % NB: the code from DaISy is wrong: % powerplant(:,1) is just the sample number Modified: trunk/octave-forge/extra/control-devel/src/slident.cc =================================================================== --- trunk/octave-forge/extra/control-devel/src/slident.cc 2012-05-12 15:06:57 UTC (rev 10426) +++ trunk/octave-forge/extra/control-devel/src/slident.cc 2012-05-12 15:32:53 UTC (rev 10427) @@ -405,8 +405,8 @@ "coefficient matrix"}; - error_msg ("ident", info_a, 2, err_msg); - warning_msg ("ident", iwarn_a, 5, warn_msg); + error_msg ("ident: IB01AD", info_a, 2, err_msg); + warning_msg ("ident: IB01AD", iwarn_a, 5, warn_msg); // resize @@ -586,8 +586,8 @@ "gain matrix is set to zero"}; - error_msg ("ident", info_b, 10, err_msg_b); - warning_msg ("ident", iwarn_b, 5, warn_msg_b); + error_msg ("ident: IB01BD", info_b, 10, err_msg_b); + warning_msg ("ident: IB01BD", iwarn_b, 5, warn_msg_b); // resize a.resize (n, n); @@ -682,8 +682,8 @@ "and/or B and D could be inaccurate"}; - error_msg ("ident", info_c, 2, err_msg_c); - warning_msg ("ident", iwarn_c, 6, warn_msg_c); + error_msg ("ident: IB01CD", info_c, 2, err_msg_c); + warning_msg ("ident: IB01CD", iwarn_c, 6, warn_msg_c); // return values This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |