Mercurial > forge
changeset 9907:ea2023de7d8f octave-forge
control-devel: debugging stuff for ib01ad
| author | paramaniac |
|---|---|
| date | Thu, 29 Mar 2012 15:07:17 +0000 |
| parents | d9c65ce686f4 |
| children | 56ae3eaf8fde |
| files | extra/control-devel/devel/PowerPlantTest.m extra/control-devel/devel/identtest.m extra/control-devel/inst/@iddata/set.m extra/control-devel/src/devel_slicot_functions.cc extra/control-devel/src/slitest.cc |
| diffstat | 5 files changed, 554 insertions(+), 2 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/control-devel/devel/PowerPlantTest.m Thu Mar 29 15:07:17 2012 +0000 @@ -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 + Peter.Vanoverschee@esat.kuleuven.ac.be +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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/control-devel/devel/identtest.m Thu Mar 29 15:07:17 2012 +0000 @@ -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
--- a/extra/control-devel/inst/@iddata/set.m Thu Mar 29 14:33:05 2012 +0000 +++ b/extra/control-devel/inst/@iddata/set.m Thu Mar 29 15:07:17 2012 +0000 @@ -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"}
--- a/extra/control-devel/src/devel_slicot_functions.cc Thu Mar 29 14:33:05 2012 +0000 +++ b/extra/control-devel/src/devel_slicot_functions.cc Thu Mar 29 15:07:17 2012 +0000 @@ -1,2 +1,3 @@ // #include "slib01ad.cc" // preprocess the input-output data #include "slident.cc" // system identification +#include "slitest.cc" // debug system identification
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/control-devel/src/slitest.cc Thu Mar 29 15:07:17 2012 +0000 @@ -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 <lukas.reichlin@gmail.com> +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; +}