Mercurial > forge
changeset 10242:52586f108b57 octave-forge
control-devel: sort variables alphabetically
| author | paramaniac |
|---|---|
| date | Sat, 12 May 2012 16:41:38 +0000 |
| parents | 2eb9871fb3e4 |
| children | 249938380b7a |
| files | extra/control-devel/devel/Destillation_combinations.m extra/control-devel/devel/PowerPlant_combinations.m extra/control-devel/devel/ident_combinations.m extra/control-devel/src/slident.cc |
| diffstat | 4 files changed, 239 insertions(+), 3 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/control-devel/devel/Destillation_combinations.m Sat May 12 16:41:38 2012 +0000 @@ -0,0 +1,90 @@ +%{ +This file describes the data in the destill.dat file. +1. Contributed by: + Peter Van Overschee + K.U.Leuven - ESAT - SISTA + K. Mercierlaan 94 + 3001 Heverlee + Peter.Vanoverschee@esat.kuleuven.ac.be +2. Process/Description: + Data of a simulation (not real !) related to the identification + of an ethane-ethylene destillationcolumn. The series consists of 4 + series: + U_dest, Y_dest: without noise (original series) + U_dest_n10, Y_dest_n10: 10 percent additive white noise + U_dest_n20, Y_dest_n20: 20 percent additive white noise + U_dest_n30, Y_dest_n30: 30 percent additive white noise +3. Sampling time + 15 min. +4. Number of samples: + 90 samples +5. Inputs: + a. ratio between the reboiler duty and the feed flow + b. ratio between the reflux rate and the feed flow + c. ratio between the distillate and the feed flow + d. input ethane composition + e. top pressure +6. Outputs: + a. top ethane composition + b. bottom ethylene composition + c. top-bottom differential pressure. +7. References: + R.P. Guidorzi, M.P. Losito, T. Muratori, The range error test in the + structural identification of linear multivariable systems, + IEEE transactions on automatic control, Vol AC-27, pp 1044-1054, oct. + 1982. +8. Known properties/peculiarities + +9. Some MATLAB-code to retrieve the data + !gunzip destill.dat.Z + load destill.dat + U=destill(:,1:20); + Y=destill(:,21:32); + U_dest=U(:,1:5); + U_dest_n10=U(:,6:10); + U_dest_n20=U(:,11:15); + U_dest_n30=U(:,16:20); + Y_dest=Y(:,1:3); + Y_dest_n10=Y(:,4:6); + Y_dest_n20=Y(:,7:9); + Y_dest_n30=Y(:,10:12); +%} + +clear all, close all, clc + +% DaISy code is wrong, +% first column is sample number +load destill.dat +U=destill(:,2:21); +Y=destill(:,22:33); +U_dest=U(:,1:5); +U_dest_n10=U(:,6:10); +U_dest_n20=U(:,11:15); +U_dest_n30=U(:,16:20); +Y_dest=Y(:,1:3); +Y_dest_n10=Y(:,4:6); +Y_dest_n20=Y(:,7:9); +Y_dest_n30=Y(:,10:12); + + +dat = iddata (Y_dest, U_dest); + +err = zeros (3, 3); + +for meth = 0:2 + for alg = 0:2 + [sys, x0] = ident_combinations (dat, 5, 4, meth, alg); % s=5, n=4 + [y, t] = lsim (sys, U_dest, [], x0); + err(meth+1, alg+1) = norm (Y_dest - y, 1) / norm (Y_dest, 1); + endfor +endfor + +err + +%{ +figure (1) +%plot (t, Y_dest, 'b') +plot (t, Y_dest, 'b', t, y, 'r') +legend ('y measured', 'y simulated', 'location', 'southeast') +%} +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/control-devel/devel/PowerPlant_combinations.m Sat May 12 16:41:38 2012 +0000 @@ -0,0 +1,92 @@ +%{ +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 + +% NB: the code from DaISy is wrong: +% powerplant(:,1) is just the sample number +% therefore increase indices by one +% it took me weeks to find that silly mistake ... +load powerplant.dat +U=powerplant(:,2:6); +Y=powerplant(:,7:9); +Yr=powerplant(:,10:12); + +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) + + +err = zeros (3, 3); + +for meth = 0:2 + for alg = 0:2 + [sys, x0] = ident_combinations (dat, 10, 8, meth, alg); % s=10, n=8 + [y, t] = lsim (sys, U, [], x0); + err(meth+1, alg+1) = norm (Y - y, 1) / norm (Y, 1); + endfor +endfor + +err + +%{ +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) +%}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/control-devel/devel/ident_combinations.m Sat May 12 16:41:38 2012 +0000 @@ -0,0 +1,54 @@ +function [sys, x0] = ident_combinations (dat, s = [], n = [], meth, alg) + + + + %nobr = 15; +% meth = 1; % 2 % geht: meth/alg 1/1, +% alg = 2; % 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 + [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, dat.tsam{1}); + +endfunction
--- a/extra/control-devel/src/slident.cc Sat May 12 15:32:53 2012 +0000 +++ b/extra/control-devel/src/slident.cc Sat May 12 16:41:38 2012 +0000 @@ -137,16 +137,16 @@ switch (imeth) { case 0: + meth_a = 'M'; meth_b = 'M'; - meth_a = 'M'; break; case 1: + meth_a = 'N'; meth_b = 'N'; - meth_a = 'N'; break; case 2: + meth_a = 'N'; // no typo here meth_b = 'C'; - meth_a = 'N'; // no typo here break; default: error ("slib01ad: argument 'meth' invalid");