This list is closed, nobody may subscribe to it.
2001 |
Jan
|
Feb
|
Mar
|
Apr
|
May
|
Jun
|
Jul
(10) |
Aug
(5) |
Sep
(3) |
Oct
(41) |
Nov
(41) |
Dec
(33) |
---|---|---|---|---|---|---|---|---|---|---|---|---|
2002 |
Jan
(75) |
Feb
(10) |
Mar
(170) |
Apr
(174) |
May
(66) |
Jun
(11) |
Jul
(10) |
Aug
(44) |
Sep
(73) |
Oct
(28) |
Nov
(139) |
Dec
(52) |
2003 |
Jan
(35) |
Feb
(93) |
Mar
(62) |
Apr
(10) |
May
(55) |
Jun
(70) |
Jul
(37) |
Aug
(16) |
Sep
(56) |
Oct
(31) |
Nov
(57) |
Dec
(83) |
2004 |
Jan
(85) |
Feb
(67) |
Mar
(27) |
Apr
(37) |
May
(75) |
Jun
(85) |
Jul
(160) |
Aug
(68) |
Sep
(104) |
Oct
(25) |
Nov
(39) |
Dec
(23) |
2005 |
Jan
(10) |
Feb
(45) |
Mar
(43) |
Apr
(19) |
May
(108) |
Jun
(31) |
Jul
(41) |
Aug
(23) |
Sep
(65) |
Oct
(58) |
Nov
(44) |
Dec
(54) |
2006 |
Jan
(96) |
Feb
(27) |
Mar
(69) |
Apr
(59) |
May
(67) |
Jun
(35) |
Jul
(13) |
Aug
(461) |
Sep
(160) |
Oct
(399) |
Nov
(32) |
Dec
(72) |
2007 |
Jan
(316) |
Feb
(305) |
Mar
(318) |
Apr
(54) |
May
(194) |
Jun
(173) |
Jul
(282) |
Aug
(91) |
Sep
(227) |
Oct
(365) |
Nov
(168) |
Dec
(18) |
2008 |
Jan
(71) |
Feb
(111) |
Mar
(155) |
Apr
(173) |
May
(70) |
Jun
(67) |
Jul
(55) |
Aug
(83) |
Sep
(32) |
Oct
(68) |
Nov
(80) |
Dec
(29) |
2009 |
Jan
(46) |
Feb
(18) |
Mar
(95) |
Apr
(76) |
May
(140) |
Jun
(98) |
Jul
(84) |
Aug
(123) |
Sep
(94) |
Oct
(131) |
Nov
(142) |
Dec
(125) |
2010 |
Jan
(128) |
Feb
(158) |
Mar
(172) |
Apr
(134) |
May
(94) |
Jun
(84) |
Jul
(32) |
Aug
(127) |
Sep
(167) |
Oct
(109) |
Nov
(69) |
Dec
(78) |
2011 |
Jan
(39) |
Feb
(58) |
Mar
(52) |
Apr
(47) |
May
(56) |
Jun
(76) |
Jul
(55) |
Aug
(54) |
Sep
(165) |
Oct
(255) |
Nov
(328) |
Dec
(263) |
2012 |
Jan
(82) |
Feb
(147) |
Mar
(400) |
Apr
(216) |
May
(209) |
Jun
(160) |
Jul
(86) |
Aug
(141) |
Sep
(156) |
Oct
(6) |
Nov
|
Dec
|
2015 |
Jan
|
Feb
|
Mar
|
Apr
|
May
|
Jun
|
Jul
(1) |
Aug
|
Sep
(1) |
Oct
|
Nov
(1) |
Dec
(2) |
2016 |
Jan
|
Feb
(2) |
Mar
(2) |
Apr
(1) |
May
(1) |
Jun
(2) |
Jul
(1) |
Aug
(1) |
Sep
|
Oct
|
Nov
(1) |
Dec
|
2019 |
Jan
|
Feb
|
Mar
(1) |
Apr
|
May
|
Jun
|
Jul
|
Aug
|
Sep
|
Oct
|
Nov
|
Dec
|
2021 |
Jan
|
Feb
|
Mar
|
Apr
(3) |
May
(4) |
Jun
(8) |
Jul
(2) |
Aug
(5) |
Sep
(9) |
Oct
|
Nov
|
Dec
|
From: <car...@us...> - 2012-05-10 11:24:41
|
Revision: 10398 http://octave.svn.sourceforge.net/octave/?rev=10398&view=rev Author: carandraug Date: 2012-05-10 11:24:35 +0000 (Thu, 10 May 2012) Log Message: ----------- inputParser/addSwitch: check for correct number of arguments Modified Paths: -------------- trunk/octave-forge/main/general/inst/@inputParser/subsref.m Modified: trunk/octave-forge/main/general/inst/@inputParser/subsref.m =================================================================== --- trunk/octave-forge/main/general/inst/@inputParser/subsref.m 2012-05-10 11:13:09 UTC (rev 10397) +++ trunk/octave-forge/main/general/inst/@inputParser/subsref.m 2012-05-10 11:24:35 UTC (rev 10398) @@ -217,8 +217,12 @@ endif def = false; case {'addSwitch'} - val = def_val; - def = false; + if ( numel (args) == 0 ) + val = def_val; + def = false; + else + print_usage(func); + endif otherwise error ("invalid index for reference of class %s", class (inPar) ); endswitch This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <car...@us...> - 2012-05-10 11:13:15
|
Revision: 10397 http://octave.svn.sourceforge.net/octave/?rev=10397&view=rev Author: carandraug Date: 2012-05-10 11:13:09 +0000 (Thu, 10 May 2012) Log Message: ----------- inputParser: bug fix on addSwitch (it was being passed to subsref as ParamValue) Modified Paths: -------------- trunk/octave-forge/main/general/inst/@inputParser/addSwitch.m Modified: trunk/octave-forge/main/general/inst/@inputParser/addSwitch.m =================================================================== --- trunk/octave-forge/main/general/inst/@inputParser/addSwitch.m 2012-05-10 07:44:00 UTC (rev 10396) +++ trunk/octave-forge/main/general/inst/@inputParser/addSwitch.m 2012-05-10 11:13:09 UTC (rev 10397) @@ -38,7 +38,7 @@ ## check @inputParser/subsref for the actual code if (nargin == 2) inPar = subsref (inPar, substruct( - '.' , 'addParamValue', + '.' , 'addSwitch', '()', {name} )); else This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2012-05-10 07:44:09
|
Revision: 10396 http://octave.svn.sourceforge.net/octave/?rev=10396&view=rev Author: paramaniac Date: 2012-05-10 07:44:00 +0000 (Thu, 10 May 2012) Log Message: ----------- control-devel: fix test examples from DaISy library, first column is sample number Modified Paths: -------------- trunk/octave-forge/extra/control-devel/devel/Destillation.m trunk/octave-forge/extra/control-devel/devel/PowerPlant.m Modified: trunk/octave-forge/extra/control-devel/devel/Destillation.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/Destillation.m 2012-05-10 07:32:55 UTC (rev 10395) +++ trunk/octave-forge/extra/control-devel/devel/Destillation.m 2012-05-10 07:44:00 UTC (rev 10396) @@ -52,9 +52,11 @@ clear all, close all, clc +% DaISy code is wrong, +% first column is sample number load destill.dat -U=destill(:,1:20); -Y=destill(:,21:32); +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); Modified: trunk/octave-forge/extra/control-devel/devel/PowerPlant.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/PowerPlant.m 2012-05-10 07:32:55 UTC (rev 10395) +++ trunk/octave-forge/extra/control-devel/devel/PowerPlant.m 2012-05-10 07:44:00 UTC (rev 10396) @@ -42,10 +42,14 @@ 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(:,1:5); -Y=powerplant(:,6:8); -Yr=powerplant(:,9:11); +U=powerplant(:,2:6); +Y=powerplant(:,7:9); +Yr=powerplant(:,10:12); inname = {'gas flow', 'turbine valves opening', This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2012-05-10 07:33:05
|
Revision: 10395 http://octave.svn.sourceforge.net/octave/?rev=10395&view=rev Author: paramaniac Date: 2012-05-10 07:32:55 +0000 (Thu, 10 May 2012) Log Message: ----------- control-devel: It looks like I've just found the mistake Modified Paths: -------------- trunk/octave-forge/extra/control-devel/devel/PowerPlantVS.m Modified: trunk/octave-forge/extra/control-devel/devel/PowerPlantVS.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/PowerPlantVS.m 2012-05-10 07:05:49 UTC (rev 10394) +++ trunk/octave-forge/extra/control-devel/devel/PowerPlantVS.m 2012-05-10 07:32:55 UTC (rev 10395) @@ -42,10 +42,14 @@ 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(:,1:5); -Y=powerplant(:,6:8); -Yr=powerplant(:,9:11); +U=powerplant(:,2:6); +Y=powerplant(:,7:9); +Yr=powerplant(:,10:12); inname = {'gas flow', 'turbine valves opening', This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2012-05-10 07:06:00
|
Revision: 10394 http://octave.svn.sourceforge.net/octave/?rev=10394&view=rev Author: paramaniac Date: 2012-05-10 07:05:49 +0000 (Thu, 10 May 2012) Log Message: ----------- control-devel: add debug code for subspace identification routines Added Paths: ----------- trunk/octave-forge/extra/control-devel/devel/CDplayerVS.m trunk/octave-forge/extra/control-devel/devel/Ident_results.zip trunk/octave-forge/extra/control-devel/devel/PowerPlantVS.m trunk/octave-forge/extra/control-devel/devel/identVS.m trunk/octave-forge/extra/control-devel/devel/overview.txt Added: trunk/octave-forge/extra/control-devel/devel/CDplayerVS.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/CDplayerVS.m (rev 0) +++ trunk/octave-forge/extra/control-devel/devel/CDplayerVS.m 2012-05-10 07:05:49 UTC (rev 10394) @@ -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 + +%} + +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] = identVS (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/Ident_results.zip =================================================================== (Binary files differ) Property changes on: trunk/octave-forge/extra/control-devel/devel/Ident_results.zip ___________________________________________________________________ Added: svn:mime-type + application/octet-stream Added: trunk/octave-forge/extra/control-devel/devel/PowerPlantVS.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/PowerPlantVS.m (rev 0) +++ trunk/octave-forge/extra/control-devel/devel/PowerPlantVS.m 2012-05-10 07:05:49 UTC (rev 10394) @@ -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] = identVS (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/identVS.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/identVS.m (rev 0) +++ trunk/octave-forge/extra/control-devel/devel/identVS.m 2012-05-10 07:05:49 UTC (rev 10394) @@ -0,0 +1,66 @@ +function [sys, x0] = identVS (dat, s = [], nn = []) + +n=nn; + + %nobr = 15; + meth = 3-1; %1-1; + alg = 1-1; + jobd = 2-1; + batch = 4-1; + conct = 2-1; + %ctrl = 0; + rcond = 0.0; + tol = -1.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 + +nobr + %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); + +%r +sv + +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 Added: trunk/octave-forge/extra/control-devel/devel/overview.txt =================================================================== --- trunk/octave-forge/extra/control-devel/devel/overview.txt (rev 0) +++ trunk/octave-forge/extra/control-devel/devel/overview.txt 2012-05-10 07:05:49 UTC (rev 10394) @@ -0,0 +1,10 @@ +Working +- CDplayer +- Evaporator +- GlassFurnace +- HeatingSystem + +Broken +- Destillation +- PowerPlant +- pH \ 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: <nir...@us...> - 2012-05-09 21:03:04
|
Revision: 10393 http://octave.svn.sourceforge.net/octave/?rev=10393&view=rev Author: nir-krakauer Date: 2012-05-09 21:02:58 +0000 (Wed, 09 May 2012) Log Message: ----------- modified cspas to extrapolate linearly Modified Paths: -------------- trunk/octave-forge/main/splines/inst/csaps.m Modified: trunk/octave-forge/main/splines/inst/csaps.m =================================================================== --- trunk/octave-forge/main/splines/inst/csaps.m 2012-05-09 19:56:20 UTC (rev 10392) +++ trunk/octave-forge/main/splines/inst/csaps.m 2012-05-09 21:02:58 UTC (rev 10393) @@ -22,6 +22,8 @@ ## ## The chosen cubic spline with natural boundary conditions @var{pp}(@var{x}) minimizes @var{p} Sum_i @var{w}_i*(@var{y}_i - @var{pp}(@var{x}_i))^2 + (1-@var{p}) Int @var{pp}''(@var{x}) d@var{x} ## +## Outside the range of @var{x}, the cubic spline is a straight line +## ## @var{x} and @var{w} should be n by 1 in size; @var{y} should be n by m; @var{xi} should be k by 1; the values in @var{x} should be distinct; the values in @var{w} should be nonzero ## ## @table @asis @@ -61,7 +63,7 @@ [x,i] = sort(x); y = y(i, :); - n = numel(x); + [n m] = size(y); #should also be that n = numel(x); if isempty(w) w = ones(n, 1); @@ -83,17 +85,24 @@ u = (6*(1-p)*QT*diag(1 ./ w)*QT' + p*R) \ (QT*y); a = y - 6*(1-p)*diag(1 ./ w)*QT'*u; -## derivatives at all but the last knot for the piecewise cubic spline - aa = a(1:(end-1), :); - cc = zeros(size(y)); - cc(2:(n-1), :) = 6*p*u; #cc([1 n], :) = 0 [natural spline] - dd = diff(cc) ./ h; - cc = cc(1:(end-1), :); - bb = diff(a) ./ h - (cc/2).*h - (dd/6).*(h.^2); - +## note: add knots to either end of spline pp-form to ensure linear extrapolation + xminus = x(1) - eps(x(1)); + xplus = x(end) + eps(x(end)); + x = [xminus; x; xplus]; + +## derivatives for the piecewise cubic spline + aa = bb = cc = dd = zeros (n+1, m); + aa(2:end, :) = a; + cc(3:n, :) = 6*p*u; #second derivative at endpoints is 0 [natural spline] + dd(2:n, :) = diff(cc(2:(n+1), :)) ./ h; + bb(2:n, :) = diff(a) ./ h - (cc(2:n, :)/2).*h - (dd(2:n, :)/6).*(h.^2); + bb(1, :) = bb(2, :); #linear extension of splines + bb(n + 1, :) = bb(n, :); + aa(1, :) = a(1, :) - eps(x(1))*bb(1, :); + ret = mkpp (x, cat (2, dd'(:)/6, cc'(:)/2, bb'(:), aa'(:)), size(y, 2)); - if ~isempty(xi) + if ~isempty (xi) ret = ppval (ret, xi); endif This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <as...@us...> - 2012-05-09 19:56:26
|
Revision: 10392 http://octave.svn.sourceforge.net/octave/?rev=10392&view=rev Author: asnelt Date: 2012-05-09 19:56:20 +0000 (Wed, 09 May 2012) Log Message: ----------- Update for release 1.1.3. Modified Paths: -------------- trunk/octave-forge/main/statistics/DESCRIPTION trunk/octave-forge/main/statistics/NEWS Modified: trunk/octave-forge/main/statistics/DESCRIPTION =================================================================== --- trunk/octave-forge/main/statistics/DESCRIPTION 2012-05-09 19:54:04 UTC (rev 10391) +++ trunk/octave-forge/main/statistics/DESCRIPTION 2012-05-09 19:56:20 UTC (rev 10392) @@ -1,6 +1,6 @@ Name: Statistics -Version: 1.1.2 -Date: 2012-05-01 +Version: 1.1.3 +Date: 2012-05-09 Author: various authors Maintainer: Arno Onken <as...@as...> Title: Statistics Modified: trunk/octave-forge/main/statistics/NEWS =================================================================== --- trunk/octave-forge/main/statistics/NEWS 2012-05-09 19:54:04 UTC (rev 10391) +++ trunk/octave-forge/main/statistics/NEWS 2012-05-09 19:56:20 UTC (rev 10392) @@ -1,6 +1,10 @@ Summary of important user-visible changes for statistics 1.1.3: ------------------------------------------------------------------- + ** The following functions are new in 1.1.3: + + copularnd mvtrnd + ** The functions mnpdf and mnrnd are now also usable for greater numbers of categories for which the rows do not exactly sum to 1. This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <as...@us...> - 2012-05-09 19:54:11
|
Revision: 10391 http://octave.svn.sourceforge.net/octave/?rev=10391&view=rev Author: asnelt Date: 2012-05-09 19:54:04 +0000 (Wed, 09 May 2012) Log Message: ----------- Initial commit. Modified Paths: -------------- trunk/octave-forge/main/statistics/INDEX Added Paths: ----------- trunk/octave-forge/main/statistics/inst/copularnd.m trunk/octave-forge/main/statistics/inst/mvtrnd.m Modified: trunk/octave-forge/main/statistics/INDEX =================================================================== --- trunk/octave-forge/main/statistics/INDEX 2012-05-09 19:12:01 UTC (rev 10390) +++ trunk/octave-forge/main/statistics/INDEX 2012-05-09 19:54:04 UTC (rev 10391) @@ -5,7 +5,7 @@ binostat chi2stat cl_multinom - copulacdf copulapdf + copulacdf copulapdf copularnd expstat fstat gamlike @@ -17,7 +17,7 @@ lognstat mvnpdf mvnrnd mvncdf mnpdf mnrnd - mvtcdf + mvtcdf mvtrnd nbinstat normalise_distribution normstat Added: trunk/octave-forge/main/statistics/inst/copularnd.m =================================================================== --- trunk/octave-forge/main/statistics/inst/copularnd.m (rev 0) +++ trunk/octave-forge/main/statistics/inst/copularnd.m 2012-05-09 19:54:04 UTC (rev 10391) @@ -0,0 +1,281 @@ +## Copyright (C) 2012 Arno Onken +## +## This program 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. +## +## This program 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 this program. If not, see <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{x} =} copularnd (@var{family}, @var{theta}, @var{n}) +## @deftypefnx {Function File} {} copularnd (@var{family}, @var{theta}, @var{n}, @var{d}) +## @deftypefnx {Function File} {} copularnd ('t', @var{theta}, @var{nu}, @var{n}) +## Generate random samples from a copula family. +## +## @subheading Arguments +## +## @itemize @bullet +## @item +## @var{family} is the copula family name. Currently, @var{family} can be +## @code{'Gaussian'} for the Gaussian family, @code{'t'} for the Student's t +## family, or @code{'Clayton'} for the Clayton family. +## +## @item +## @var{theta} is the parameter of the copula. For the Gaussian and Student's t +## copula, @var{theta} must be a correlation matrix. For bivariate copulas +## @var{theta} can also be a correlation coefficient. For the Clayton family, +## @var{theta} must be a vector with the same number of elements as samples to +## be generated or be scalar. +## +## @item +## @var{nu} is the degrees of freedom for the Student's t family. @var{nu} must +## be a vector with the same number of elements as samples to be generated or +## be scalar. +## +## @item +## @var{n} is the number of rows of the matrix to be generated. @var{n} must be +## a non-negative integer and corresponds to the number of samples to be +## generated. +## +## @item +## @var{d} is the number of columns of the matrix to be generated. @var{d} must +## be a positive integer and corresponds to the dimension of the copula. +## @end itemize +## +## @subheading Return values +## +## @itemize @bullet +## @item +## @var{x} is a matrix of random samples from the copula with @var{n} samples +## of distribution dimension @var{d}. +## @end itemize +## +## @subheading Examples +## +## @example +## @group +## theta = 0.5; +## x = copularnd ("Gaussian", theta); +## @end group +## +## @group +## theta = 0.5; +## nu = 2; +## x = copularnd ("t", theta, nu); +## @end group +## +## @group +## theta = 0.5; +## n = 2; +## x = copularnd ("Clayton", theta, n); +## @end group +## @end example +## +## @subheading References +## +## @enumerate +## @item +## Roger B. Nelsen. @cite{An Introduction to Copulas}. Springer, New York, +## second edition, 2006. +## @end enumerate +## @end deftypefn + +## Author: Arno Onken <as...@as...> +## Description: Random samples from a copula family + +function x = copularnd (family, theta, nu, n) + + # Check arguments + if (nargin < 2) + print_usage (); + endif + + if (! ischar (family)) + error ("copularnd: family must be one of 'Gaussian', 't', and 'Clayton'"); + endif + + lower_family = lower (family); + + # Check family and copula parameters + switch (lower_family) + + case {"gaussian"} + # Gaussian family + if (isscalar (theta)) + # Expand a scalar to a correlation matrix + theta = [1, theta; theta, 1]; + endif + if (! ismatrix (theta) || any (diag (theta) != 1) || any (any (theta != theta')) || min (eig (theta)) <= 0) + error ("copularnd: theta must be a correlation matrix"); + endif + if (nargin > 3) + d = n; + if (! isscalar (d) || d != size (theta, 1)) + error ("copularnd: d must correspond to dimension of theta"); + endif + else + d = size (theta, 1); + endif + if (nargin < 3) + n = 1; + else + n = nu; + if (! isscalar (n) || (n < 0) || round (n) != n) + error ("copularnd: n must be a non-negative integer"); + endif + endif + + case {"t"} + # Student's t family + if (nargin < 3) + print_usage (); + endif + if (isscalar (theta)) + # Expand a scalar to a correlation matrix + theta = [1, theta; theta, 1]; + endif + if (! ismatrix (theta) || any (diag (theta) != 1) || any (any (theta != theta')) || min (eig (theta)) <= 0) + error ("copularnd: theta must be a correlation matrix"); + endif + if (! isscalar (nu) && (! isvector (nu) || length (nu) != n)) + error ("copularnd: nu must be a vector with the same number of rows as x or be scalar"); + endif + nu = nu(:); + if (nargin < 4) + n = 1; + else + if (! isscalar (n) || (n < 0) || round (n) != n) + error ("copularnd: n must be a non-negative integer"); + endif + endif + + case {"clayton"} + # Archimedian one parameter family + if (nargin < 4) + # Default is bivariate + d = 2; + else + d = n; + if (! isscalar (d) || (d < 2) || round (d) != d) + error ("copularnd: d must be an integer greater than 1"); + endif + endif + if (nargin < 3) + # Default is one sample + n = 1; + else + n = nu; + if (! isscalar (n) || (n < 0) || round (n) != n) + error ("copularnd: n must be a non-negative integer"); + endif + endif + if (! isvector (theta) || (! isscalar (theta) && size (theta, 1) != n)) + error ("copularnd: theta must be a column vector with the number of rows equal to n or be scalar"); + endif + if (n > 1 && isscalar (theta)) + theta = repmat (theta, n, 1); + endif + + otherwise + error ("copularnd: unknown copula family '%s'", family); + + endswitch + + if (n == 0) + # Input is empty + x = zeros (0, d); + else + + # Draw random samples according to family + switch (lower_family) + + case {"gaussian"} + # The Gaussian family + x = normcdf (mvnrnd (zeros (1, d), theta, n), 0, 1); + # No parameter bounds check + k = []; + + case {"t"} + # The Student's t family + x = tcdf (mvtrnd (theta, nu, n), nu); + # No parameter bounds check + k = []; + + case {"clayton"} + # The Clayton family + u = rand (n, d); + if (d == 2) + x = zeros (n, 2); + # Conditional distribution method for the bivariate case which also + # works for theta < 0 + x(:, 1) = u(:, 1); + x(:, 2) = (1 + u(:, 1) .^ (-theta) .* (u(:, 2) .^ (-theta ./ (1 + theta)) - 1)) .^ (-1 ./ theta); + else + # Apply the algorithm by Marshall and Olkin: + # Frailty distribution for Clayton copula is gamma + y = randg (1 ./ theta, n, 1); + x = (1 - log (u) ./ repmat (y, 1, d)) .^ (-1 ./ repmat (theta, 1, d)); + endif + k = find (theta == 0); + if (any (k)) + # Produkt copula at columns k + x(k, :) = u(k, :); + endif + # Continue argument check + if (d == 2) + k = find (! (theta >= -1) | ! (theta < inf)); + else + k = find (! (theta >= 0) | ! (theta < inf)); + endif + + endswitch + + # Out of bounds parameters + if (any (k)) + x(k, :) = NaN; + endif + + endif + +endfunction + +%!test +%! theta = 0.5; +%! x = copularnd ("Gaussian", theta); +%! assert (size (x), [1, 2]); +%! assert (all ((x >= 0) & (x <= 1))); + +%!test +%! theta = 0.5; +%! nu = 2; +%! x = copularnd ("t", theta, nu); +%! assert (size (x), [1, 2]); +%! assert (all ((x >= 0) & (x <= 1))); + +%!test +%! theta = 0.5; +%! x = copularnd ("Clayton", theta); +%! assert (size (x), [1, 2]); +%! assert (all ((x >= 0) & (x <= 1))); + +%!test +%! theta = 0.5; +%! n = 2; +%! x = copularnd ("Clayton", theta, n); +%! assert (size (x), [n, 2]); +%! assert (all ((x >= 0) & (x <= 1))); + +%!test +%! theta = [1; 2]; +%! n = 2; +%! d = 3; +%! x = copularnd ("Clayton", theta, n, d); +%! assert (size (x), [n, d]); +%! assert (all ((x >= 0) & (x <= 1))); Added: trunk/octave-forge/main/statistics/inst/mvtrnd.m =================================================================== --- trunk/octave-forge/main/statistics/inst/mvtrnd.m (rev 0) +++ trunk/octave-forge/main/statistics/inst/mvtrnd.m 2012-05-09 19:54:04 UTC (rev 10391) @@ -0,0 +1,137 @@ +## Copyright (C) 2012 Arno Onken <as...@as...> +## +## This program 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. +## +## This program 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 +## this program; if not, see <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{x} =} mvtrnd (@var{sigma}, @var{nu}) +## @deftypefnx {Function File} {@var{x} =} mvtrnd (@var{sigma}, @var{nu}, @var{n}) +## Generate random samples from the multivariate t-distribution. +## +## @subheading Arguments +## +## @itemize @bullet +## @item +## @var{sigma} is the matrix of correlation coefficients. If there are any +## non-unit diagonal elements then @var{sigma} will be normalized. +## +## @item +## @var{nu} is the degrees of freedom for the multivariate t-distribution. +## @var{nu} must be a vector with the same number of elements as samples to be +## generated or be scalar. +## +## @item +## @var{n} is the number of rows of the matrix to be generated. @var{n} must be +## a non-negative integer and corresponds to the number of samples to be +## generated. +## @end itemize +## +## @subheading Return values +## +## @itemize @bullet +## @item +## @var{x} is a matrix of random samples from the multivariate t-distribution +## with @var{n} row samples. +## @end itemize +## +## @subheading Examples +## +## @example +## @group +## sigma = [1, 0.5; 0.5, 1]; +## nu = 3; +## n = 10; +## x = mvtrnd (sigma, nu, n); +## @end group +## +## @group +## sigma = [1, 0.5; 0.5, 1]; +## nu = [2; 3]; +## n = 2; +## x = mvtrnd (sigma, nu, 2); +## @end group +## @end example +## +## @subheading References +## +## @enumerate +## @item +## Wendy L. Martinez and Angel R. Martinez. @cite{Computational Statistics +## Handbook with MATLAB}. Appendix E, pages 547-557, Chapman & Hall/CRC, 2001. +## +## @item +## Samuel Kotz and Saralees Nadarajah. @cite{Multivariate t Distributions and +## Their Applications}. Cambridge University Press, Cambridge, 2004. +## @end enumerate +## @end deftypefn + +## Author: Arno Onken <as...@as...> +## Description: Random samples from the multivariate t-distribution + +function x = mvtrnd (sigma, nu, n) + + # Check arguments + if (nargin < 2) + print_usage (); + endif + + if (! ismatrix (sigma) || any (any (sigma != sigma')) || min (eig (sigma)) <= 0) + error ("mvtrnd: sigma must be a positive definite matrix"); + endif + + if (!isvector (nu) || any (nu <= 0)) + error ("mvtrnd: nu must be a positive scalar or vector"); + endif + nu = nu(:); + + if (nargin > 2) + if (! isscalar (n) || n < 0 | round (n) != n) + error ("mvtrnd: n must be a non-negative integer") + endif + if (isscalar (nu)) + nu = nu * ones (n, 1); + else + if (length (nu) != n) + error ("mvtrnd: n must match the length of nu") + endif + endif + else + n = length (nu); + endif + + # Normalize sigma + if (any (diag (sigma) != 1)) + sigma = sigma ./ sqrt (diag (sigma) * diag (sigma)'); + endif + + # Dimension + d = size (sigma, 1); + # Draw samples + y = mvnrnd (zeros (1, d), sigma, n); + u = repmat (chi2rnd (nu), 1, d); + x = y .* sqrt (repmat (nu, 1, d) ./ u); +endfunction + +%!test +%! sigma = [1, 0.5; 0.5, 1]; +%! nu = 3; +%! n = 10; +%! x = mvtrnd (sigma, nu, n); +%! assert (size (x), [10, 2]); + +%!test +%! sigma = [1, 0.5; 0.5, 1]; +%! nu = [2; 3]; +%! n = 2; +%! x = mvtrnd (sigma, nu, 2); +%! assert (size (x), [2, 2]); This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2012-05-09 19:12:07
|
Revision: 10390 http://octave.svn.sourceforge.net/octave/?rev=10390&view=rev Author: paramaniac Date: 2012-05-09 19:12:01 +0000 (Wed, 09 May 2012) Log Message: ----------- control: reactivate dss minreal test Modified Paths: -------------- trunk/octave-forge/main/control/inst/@lti/minreal.m Modified: trunk/octave-forge/main/control/inst/@lti/minreal.m =================================================================== --- trunk/octave-forge/main/control/inst/@lti/minreal.m 2012-05-09 11:33:45 UTC (rev 10389) +++ trunk/octave-forge/main/control/inst/@lti/minreal.m 2012-05-09 19:12:01 UTC (rev 10390) @@ -92,80 +92,80 @@ ## dss: minreal (SLICOT TG01JD) ## FIXME: Test fails with larger ldwork in sltg01jd.cc -#%!shared Ar, Br, Cr, Dr, Er, Ae, Be, Ce, De, Ee -#%! A = [ -2 -3 0 0 0 0 0 0 0 -#%! 1 0 0 0 0 0 0 0 0 -#%! 0 0 -2 -3 0 0 0 0 0 -#%! 0 0 1 0 0 0 0 0 0 -#%! 0 0 0 0 1 0 0 0 0 -#%! 0 0 0 0 0 1 0 0 0 -#%! 0 0 0 0 0 0 1 0 0 -#%! 0 0 0 0 0 0 0 1 0 -#%! 0 0 0 0 0 0 0 0 1 ]; -#%! -#%! E = [ 1 0 0 0 0 0 0 0 0 -#%! 0 1 0 0 0 0 0 0 0 -#%! 0 0 1 0 0 0 0 0 0 -#%! 0 0 0 1 0 0 0 0 0 -#%! 0 0 0 0 0 0 0 0 0 -#%! 0 0 0 0 1 0 0 0 0 -#%! 0 0 0 0 0 0 0 0 0 -#%! 0 0 0 0 0 0 1 0 0 -#%! 0 0 0 0 0 0 0 1 0 ]; -#%! -#%! B = [ 1 0 -#%! 0 0 -#%! 0 1 -#%! 0 0 -#%! -1 0 -#%! 0 0 -#%! 0 -1 -#%! 0 0 -#%! 0 0 ]; -#%! -#%! C = [ 1 0 1 -3 0 1 0 2 0 -#%! 0 1 1 3 0 1 0 0 1 ]; -#%! -#%! D = zeros (2, 2); -#%! -#%! sys = dss (A, B, C, D, E, "scaled", true); -#%! sysmin = minreal (sys, 0.0); -#%! [Ar, Br, Cr, Dr, Er] = dssdata (sysmin); -#%! -#%! Ae = [ 1.0000 -0.0393 -0.0980 -0.1066 0.0781 -0.2330 0.0777 -#%! 0.0000 1.0312 0.2717 0.2609 -0.1533 0.6758 -0.3553 -#%! 0.0000 0.0000 1.3887 0.6699 -0.4281 1.6389 -0.7615 -#%! 0.0000 0.0000 0.0000 -1.2147 0.2423 -0.9792 0.4788 -#%! 0.0000 0.0000 0.0000 0.0000 -1.0545 0.5035 -0.2788 -#%! 0.0000 0.0000 0.0000 0.0000 0.0000 1.6355 -0.4323 -#%! 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 1.0000 ]; -#%! -#%! Ee = [ 0.4100 0.2590 0.5080 -0.3109 0.0705 0.1429 -0.1477 -#%! -0.7629 -0.3464 0.0992 -0.3007 0.0619 0.2483 -0.0152 -#%! 0.1120 -0.2124 -0.4184 -0.1288 0.0569 -0.4213 -0.6182 -#%! 0.0000 0.1122 -0.0039 0.2771 -0.0758 0.0975 0.3923 -#%! 0.0000 0.0000 0.3708 -0.4290 0.1006 0.1402 -0.2699 -#%! 0.0000 0.0000 0.0000 0.0000 0.9458 -0.2211 0.2378 -#%! 0.0000 0.0000 0.0000 0.5711 0.2648 0.5948 -0.5000 ]; -#%! -#%! Be = [ -0.5597 0.2363 -#%! -0.4843 -0.0498 -#%! -0.4727 -0.1491 -#%! 0.1802 1.1574 -#%! 0.5995 0.1556 -#%! -0.1729 -0.3999 -#%! 0.0000 0.2500 ]; -#%! -#%! Ce = [ 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 4.0000 -#%! 0.0000 0.0000 0.0000 0.0000 0.0000 3.1524 -1.7500 ]; -#%! -#%! De = zeros (2, 2); -#%! -#%!assert (Ar, Ae, 1e-4); -#%!assert (Br, Be, 1e-4); -#%!assert (Cr, Ce, 1e-4); -#%!assert (Dr, De, 1e-4); -#%!assert (Er, Ee, 1e-4); +%!shared Ar, Br, Cr, Dr, Er, Ae, Be, Ce, De, Ee +%! A = [ -2 -3 0 0 0 0 0 0 0 +%! 1 0 0 0 0 0 0 0 0 +%! 0 0 -2 -3 0 0 0 0 0 +%! 0 0 1 0 0 0 0 0 0 +%! 0 0 0 0 1 0 0 0 0 +%! 0 0 0 0 0 1 0 0 0 +%! 0 0 0 0 0 0 1 0 0 +%! 0 0 0 0 0 0 0 1 0 +%! 0 0 0 0 0 0 0 0 1 ]; +%! +%! E = [ 1 0 0 0 0 0 0 0 0 +%! 0 1 0 0 0 0 0 0 0 +%! 0 0 1 0 0 0 0 0 0 +%! 0 0 0 1 0 0 0 0 0 +%! 0 0 0 0 0 0 0 0 0 +%! 0 0 0 0 1 0 0 0 0 +%! 0 0 0 0 0 0 0 0 0 +%! 0 0 0 0 0 0 1 0 0 +%! 0 0 0 0 0 0 0 1 0 ]; +%! +%! B = [ 1 0 +%! 0 0 +%! 0 1 +%! 0 0 +%! -1 0 +%! 0 0 +%! 0 -1 +%! 0 0 +%! 0 0 ]; +%! +%! C = [ 1 0 1 -3 0 1 0 2 0 +%! 0 1 1 3 0 1 0 0 1 ]; +%! +%! D = zeros (2, 2); +%! +%! sys = dss (A, B, C, D, E, "scaled", true); +%! sysmin = minreal (sys, 0.0); +%! [Ar, Br, Cr, Dr, Er] = dssdata (sysmin); +%! +%! Ae = [ 1.0000 -0.0393 -0.0980 -0.1066 0.0781 -0.2330 0.0777 +%! 0.0000 1.0312 0.2717 0.2609 -0.1533 0.6758 -0.3553 +%! 0.0000 0.0000 1.3887 0.6699 -0.4281 1.6389 -0.7615 +%! 0.0000 0.0000 0.0000 -1.2147 0.2423 -0.9792 0.4788 +%! 0.0000 0.0000 0.0000 0.0000 -1.0545 0.5035 -0.2788 +%! 0.0000 0.0000 0.0000 0.0000 0.0000 1.6355 -0.4323 +%! 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 1.0000 ]; +%! +%! Ee = [ 0.4100 0.2590 0.5080 -0.3109 0.0705 0.1429 -0.1477 +%! -0.7629 -0.3464 0.0992 -0.3007 0.0619 0.2483 -0.0152 +%! 0.1120 -0.2124 -0.4184 -0.1288 0.0569 -0.4213 -0.6182 +%! 0.0000 0.1122 -0.0039 0.2771 -0.0758 0.0975 0.3923 +%! 0.0000 0.0000 0.3708 -0.4290 0.1006 0.1402 -0.2699 +%! 0.0000 0.0000 0.0000 0.0000 0.9458 -0.2211 0.2378 +%! 0.0000 0.0000 0.0000 0.5711 0.2648 0.5948 -0.5000 ]; +%! +%! Be = [ -0.5597 0.2363 +%! -0.4843 -0.0498 +%! -0.4727 -0.1491 +%! 0.1802 1.1574 +%! 0.5995 0.1556 +%! -0.1729 -0.3999 +%! 0.0000 0.2500 ]; +%! +%! Ce = [ 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 4.0000 +%! 0.0000 0.0000 0.0000 0.0000 0.0000 3.1524 -1.7500 ]; +%! +%! De = zeros (2, 2); +%! +%!assert (Ar, Ae, 1e-4); +%!assert (Br, Be, 1e-4); +%!assert (Cr, Ce, 1e-4); +%!assert (Dr, De, 1e-4); +%!assert (Er, Ee, 1e-4); ## tf: minreal This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2012-05-09 11:33:52
|
Revision: 10389 http://octave.svn.sourceforge.net/octave/?rev=10389&view=rev Author: paramaniac Date: 2012-05-09 11:33:45 +0000 (Wed, 09 May 2012) Log Message: ----------- control: update news file Modified Paths: -------------- trunk/octave-forge/main/control/NEWS Modified: trunk/octave-forge/main/control/NEWS =================================================================== --- trunk/octave-forge/main/control/NEWS 2012-05-09 11:20:59 UTC (rev 10388) +++ trunk/octave-forge/main/control/NEWS 2012-05-09 11:33:45 UTC (rev 10389) @@ -4,11 +4,25 @@ control-2.3.51 Release Date: 2012-xx-yy Release Manager: Lukas Reichlin =============================================================================== -** filt, filtdata +** filt, filtdata, tf -- Added function "filt" to specify disrete-time transfer functions in DSP - format. + format, i.e. z^-1. -- Added function "filtdata" to return any type of discrete-time LTI model in DSP format. + -- tf models have a new property "inv". To display a discrete-time TF sys + in z^-1, set sys.inv=true. In order to switch to z, set sys.inv=false. + "filt" sets property "inv" to true (z^-1) by default, while "tf" uses + false (z) as default value. + +** ctranspose + Conjugate transpose or pertransposition of LTI objects. + Used by Octave for "sys’". For a transfer-function matrix G, G’ denotes the + conjugate of G given by G.’(-s) for a continuous-time system or G.’(1/z) for + a discrete-time system. The frequency response of the pertransposition of G + is the Hermitian (conjugate) transpose of G(jw), + i.e. freqresp (G’, w) = freqresp (G, w)’. + WARNING: Do NOT use this for dual problems, use the transpose "sys.’" + (note the dot) instead. =============================================================================== This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2012-05-09 11:21:10
|
Revision: 10388 http://octave.svn.sourceforge.net/octave/?rev=10388&view=rev Author: paramaniac Date: 2012-05-09 11:20:59 +0000 (Wed, 09 May 2012) Log Message: ----------- control: make DSP format dominant (due to implementation, the sum of two TF in z^-1 format would be in z again if it were recessive) Modified Paths: -------------- trunk/octave-forge/main/control/inst/@tf/__sys_group__.m Modified: trunk/octave-forge/main/control/inst/@tf/__sys_group__.m =================================================================== --- trunk/octave-forge/main/control/inst/@tf/__sys_group__.m 2012-05-09 10:36:17 UTC (rev 10387) +++ trunk/octave-forge/main/control/inst/@tf/__sys_group__.m 2012-05-09 11:20:59 UTC (rev 10388) @@ -61,7 +61,7 @@ retsys.tfvar = sys1.tfvar; endif - if (sys1.inv && sys2.inv) + if (sys1.inv || sys2.inv) retsys.inv = true; endif 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: <mtm...@us...> - 2012-05-09 03:49:16
|
Revision: 10386 http://octave.svn.sourceforge.net/octave/?rev=10386&view=rev Author: mtmiller Date: 2012-05-09 03:49:09 +0000 (Wed, 09 May 2012) Log Message: ----------- signal: make window length argument consistent across all functions Modified Paths: -------------- trunk/octave-forge/main/signal/inst/chebwin.m trunk/octave-forge/main/signal/inst/flattopwin.m trunk/octave-forge/main/signal/inst/gausswin.m trunk/octave-forge/main/signal/inst/kaiser.m trunk/octave-forge/main/signal/inst/rectwin.m trunk/octave-forge/main/signal/inst/triang.m trunk/octave-forge/main/signal/inst/tukeywin.m Modified: trunk/octave-forge/main/signal/inst/chebwin.m =================================================================== --- trunk/octave-forge/main/signal/inst/chebwin.m 2012-05-09 01:25:06 UTC (rev 10385) +++ trunk/octave-forge/main/signal/inst/chebwin.m 2012-05-09 03:49:09 UTC (rev 10386) @@ -13,9 +13,9 @@ ## You should have received a copy of the GNU General Public License along with ## this program; if not, see <http://www.gnu.org/licenses/>. -## Usage: chebwin (n, at) +## Usage: chebwin (L, at) ## -## Returns the filter coefficients of the n-point Dolph-Chebyshev window +## Returns the filter coefficients of the L-point Dolph-Chebyshev window ## with at dB of attenuation in the stop-band of the corresponding ## Fourier transform. ## @@ -31,13 +31,13 @@ ## ## The window is described in frequency domain by the expression: ## -## Cheb(n-1, beta * cos(pi * k/n)) +## Cheb(L-1, beta * cos(pi * k/L)) ## W(k) = ------------------------------- -## Cheb(n-1, beta) +## Cheb(L-1, beta) ## ## with ## -## beta = cosh(1/(n-1) * acosh(10^(at/20)) +## beta = cosh(1/(L-1) * acosh(10^(at/20)) ## ## and Cheb(m,x) denoting the m-th order Chebyshev polynomial calculated ## at the point x. @@ -48,38 +48,38 @@ ## ## See also: kaiser -function w = chebwin (n, at) +function w = chebwin (L, at) if (nargin != 2) print_usage; - elseif !(isscalar (n) && (n == round(n)) && (n > 0)) - error ("chebwin: n has to be a positive integer"); + elseif !(isscalar (L) && (L == round(L)) && (L > 0)) + error ("chebwin: L has to be a positive integer"); elseif !(isscalar (at) && (at == real (at))) error ("chebwin: at has to be a real scalar"); endif - if (n == 1) + if (L == 1) w = 1; else # beta calculation gamma = 10^(-at/20); - beta = cosh(1/(n-1) * acosh(1/gamma)); + beta = cosh(1/(L-1) * acosh(1/gamma)); # freq. scale - k = (0:n-1); - x = beta*cos(pi*k/n); + k = (0:L-1); + x = beta*cos(pi*k/L); # Chebyshev window (freq. domain) - p = cheb(n-1, x); + p = cheb(L-1, x); # inverse Fourier transform - if (rem(n,2)) + if (rem(L,2)) w = real(fft(p)); - M = (n+1)/2; + M = (L+1)/2; w = w(1:M)/w(1); w = [w(M:-1:2) w]'; else # half-sample delay (even order) - p = p.*exp(j*pi/n * (0:n-1)); + p = p.*exp(j*pi/L * (0:L-1)); w = real(fft(p)); - M = n/2+1; + M = L/2+1; w = w/w(2); w = [w(M:-1:2) w(2:M)]'; endif Modified: trunk/octave-forge/main/signal/inst/flattopwin.m =================================================================== --- trunk/octave-forge/main/signal/inst/flattopwin.m 2012-05-09 01:25:06 UTC (rev 10385) +++ trunk/octave-forge/main/signal/inst/flattopwin.m 2012-05-09 03:49:09 UTC (rev 10386) @@ -1,15 +1,15 @@ ## Author: Paul Kienzle <pki...@us...> (2004) ## This program is granted to the public domain. -## flattopwin(n, [periodic|symmetric]) +## flattopwin(L, [periodic|symmetric]) ## ## Return the window f(w): ## ## f(w) = 1 - 1.93 cos(2 pi w) + 1.29 cos(4 pi w) ## - 0.388 cos(6 pi w) + 0.0322cos(8 pi w) ## -## where w = i/(n-1) for i=0:n-1 for a symmetric window, or -## w = i/n for i=0:n-1 for a periodic window. The default +## where w = i/(L-1) for i=0:L-1 for a symmetric window, or +## w = i/L for i=0:L-1 for a periodic window. The default ## is symmetric. The returned window is normalized to a peak ## of 1 at w = 0.5. ## @@ -23,23 +23,23 @@ ## [1] Gade, S; Herlufsen, H; (1987) "Use of weighting functions in DFT/FFT ## analysis (Part I)", Bruel & Kjaer Technical Review No.3. -function w = flattopwin (n, sym) +function w = flattopwin (L, sym) if nargin == 0 || nargin > 2 print_usage; endif - divisor = n-1; + divisor = L-1; if nargin > 1 match = strmatch(sym,['periodic';'symmetric']); if isempty(match), error("window type must be periodic or symmetric"); elseif match == 1 - divisor = n; + divisor = L; else - divisor = n-1; + divisor = L-1; endif endif - x = 2*pi*[0:n-1]'/divisor; + x = 2*pi*[0:L-1]'/divisor; w = (1-1.93*cos(x)+1.29*cos(2*x)-0.388*cos(3*x)+0.0322*cos(4*x))/4.6402; endfunction Modified: trunk/octave-forge/main/signal/inst/gausswin.m =================================================================== --- trunk/octave-forge/main/signal/inst/gausswin.m 2012-05-09 01:25:06 UTC (rev 10385) +++ trunk/octave-forge/main/signal/inst/gausswin.m 2012-05-09 03:49:09 UTC (rev 10386) @@ -13,21 +13,21 @@ ## You should have received a copy of the GNU General Public License along with ## this program; if not, see <http://www.gnu.org/licenses/>. -## usage: w = gausswin(n, a) +## usage: w = gausswin(L, a) ## -## Generate an n-point gaussian window of the given width. Use larger a -## for a narrow window. Use larger n for a smoother curve. +## Generate an L-point gaussian window of the given width. Use larger a +## for a narrow window. Use larger L for a smoother curve. ## ## w = exp ( -(a*x)^2/2 ) ## -## for x = linspace(-(n-1)/n, (n-1)/n, n) +## for x = linspace(-(L-1)/L, (L-1)/L, L) -function x = gausswin(n, w) +function x = gausswin(L, w) if nargin < 1 || nargin > 2 print_usage; end if nargin == 1, w = 2.5; endif - x = exp ( -0.5 * ( w/n * [ -(n-1) : 2 : n-1 ]' ) .^ 2 ); + x = exp ( -0.5 * ( w/L * [ -(L-1) : 2 : L-1 ]' ) .^ 2 ); endfunction Modified: trunk/octave-forge/main/signal/inst/kaiser.m =================================================================== --- trunk/octave-forge/main/signal/inst/kaiser.m 2012-05-09 01:25:06 UTC (rev 10385) +++ trunk/octave-forge/main/signal/inst/kaiser.m 2012-05-09 03:49:09 UTC (rev 10386) @@ -14,36 +14,36 @@ ## You should have received a copy of the GNU General Public License along with ## this program; if not, see <http://www.gnu.org/licenses/>. -## usage: kaiser (n, beta) +## usage: kaiser (L, beta) ## -## Returns the filter coefficients of the n-point Kaiser window with +## Returns the filter coefficients of the L-point Kaiser window with ## parameter beta. ## ## For the definition of the Kaiser window, see A. V. Oppenheim & ## R. W. Schafer, "Discrete-Time Signal Processing". ## -## The continuous version of width n centered about x=0 is: +## The continuous version of width L centered about x=0 is: ## -## besseli(0, beta * sqrt(1-(2*x/n).^2)) -## k(x) = -------------------------------------, n/2 <= x <= n/2 +## besseli(0, beta * sqrt(1-(2*x/L).^2)) +## k(x) = -------------------------------------, L/2 <= x <= L/2 ## besseli(0, beta) ## ## See also: kaiserord -function w = kaiser (n, beta = 0.5) +function w = kaiser (L, beta = 0.5) if (nargin < 1) print_usage; - elseif !(isscalar (n) && (n == round (n)) && (n > 0)) - error ("kaiser: n has to be a positive integer"); + elseif !(isscalar (L) && (L == round (L)) && (L > 0)) + error ("kaiser: L has to be a positive integer"); elseif !(isscalar (beta) && (beta == real (beta))) error ("kaiser: beta has to be a real scalar"); endif - if (n == 1) + if (L == 1) w = 1; else - m = n - 1; + m = L - 1; k = (0 : m)'; k = 2 * beta / m * sqrt (k .* (m - k)); w = besseli (0, k) / besseli (0, beta); Modified: trunk/octave-forge/main/signal/inst/rectwin.m =================================================================== --- trunk/octave-forge/main/signal/inst/rectwin.m 2012-05-09 01:25:06 UTC (rev 10385) +++ trunk/octave-forge/main/signal/inst/rectwin.m 2012-05-09 03:49:09 UTC (rev 10386) @@ -14,12 +14,12 @@ ## this program; if not, see <http://www.gnu.org/licenses/>. ## -*- texinfo -*- -## @deftypefn {Function File} {[@var{w}] =} rectwin(@var{n}) -## Return the filter coefficients of a rectangle window of length N. +## @deftypefn {Function File} {[@var{w}] =} rectwin(@var{L}) +## Return the filter coefficients of a rectangle window of length L. ## @seealso{hamming, hanning} ## @end deftypefn -function w = rectwin(n) +function w = rectwin(L) if (nargin < 1); print_usage; end - w = ones(round(n),1); + w = ones(round(L),1); endfunction Modified: trunk/octave-forge/main/signal/inst/triang.m =================================================================== --- trunk/octave-forge/main/signal/inst/triang.m 2012-05-09 01:25:06 UTC (rev 10385) +++ trunk/octave-forge/main/signal/inst/triang.m 2012-05-09 03:49:09 UTC (rev 10386) @@ -13,20 +13,20 @@ ## You should have received a copy of the GNU General Public License along with ## this program; if not, see <http://www.gnu.org/licenses/>. -## usage: w = triang (n) +## usage: w = triang (L) ## -## Returns the filter coefficients of a triangular window of length n. +## Returns the filter coefficients of a triangular window of length L. ## Unlike the bartlett window, triang does not go to zero at the edges -## of the window. For odd n, triang(n) is equal to bartlett(n+2) except +## of the window. For odd L, triang(L) is equal to bartlett(L+2) except ## for the zeros at the edges of the window. -function w = triang(n) +function w = triang(L) if (nargin != 1) print_usage; - elseif (!isscalar(n) || n != fix (n) || n < 1) - error("triang(n): n has to be an integer > 0"); + elseif (!isscalar(L) || L != fix (L) || L < 1) + error("triang: L has to be an integer > 0"); endif - w = 1 - abs ([-(n-1):2:(n-1)]' / (n+rem(n,2))); + w = 1 - abs ([-(L-1):2:(L-1)]' / (L+rem(L,2))); endfunction %!error triang Modified: trunk/octave-forge/main/signal/inst/tukeywin.m =================================================================== --- trunk/octave-forge/main/signal/inst/tukeywin.m 2012-05-09 01:25:06 UTC (rev 10385) +++ trunk/octave-forge/main/signal/inst/tukeywin.m 2012-05-09 03:49:09 UTC (rev 10386) @@ -14,9 +14,9 @@ ## this program; if not, see <http://www.gnu.org/licenses/>. ## -*- texinfo -*- -## @deftypefn {Function File} {@var{w} =} tukeywin (@var{m}, @var{r}) +## @deftypefn {Function File} {@var{w} =} tukeywin (@var{L}, @var{r}) ## Return the filter coefficients of a Tukey window (also known as the -## cosine-tapered window) of length @var{m}. @var{r} defines the ratio +## cosine-tapered window) of length @var{L}. @var{r} defines the ratio ## between the constant section and and the cosine section. It has to be ## between 0 and 1. The function returns a Hanning window for @var{r} ## egals 0 and a full box for @var{r} egals 1. By default @var{r} is set @@ -28,7 +28,7 @@ ## Page 67, Equation 38. ## @end deftypefn -function w = tukeywin (m, r = 1/2) +function w = tukeywin (L, r = 1/2) if (nargin < 1 || nargin > 2) print_usage; @@ -45,23 +45,23 @@ switch r case 0, ## full box - w = ones (m, 1); + w = ones (L, 1); case 1, ## Hanning window - w = hanning (m); + w = hanning (L); otherwise ## cosine-tapered window - t = linspace(0,1,m)(1:end/2)'; + t = linspace(0,1,L)(1:end/2)'; w = (1 + cos(pi*(2*t/r-1)))/2; - w(floor(r*(m-1)/2)+2:end) = 1; - w = [w; ones(mod(m,2)); flipud(w)]; + w(floor(r*(L-1)/2)+2:end) = 1; + w = [w; ones(mod(L,2)); flipud(w)]; endswitch endfunction %!demo -%! m = 100; +%! L = 100; %! r = 1/3; -%! w = tukeywin (m, r); -%! title(sprintf("%d-point Tukey window, R = %d/%d", m, [p, q] = rat(r), q)); +%! w = tukeywin (L, r); +%! title(sprintf("%d-point Tukey window, R = %d/%d", L, [p, q] = rat(r), q)); %! plot(w); This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <mtm...@us...> - 2012-05-09 01:25:14
|
Revision: 10385 http://octave.svn.sourceforge.net/octave/?rev=10385&view=rev Author: mtmiller Date: 2012-05-09 01:25:06 +0000 (Wed, 09 May 2012) Log Message: ----------- blackmanharris: fix window symmetry Submitted by Mike Gross <mi...@ap...>, resubmitted by Daniel Sebald <dan...@ie...>. Closes artifact #3524658 Modified Paths: -------------- trunk/octave-forge/main/signal/inst/blackmanharris.m Modified: trunk/octave-forge/main/signal/inst/blackmanharris.m =================================================================== --- trunk/octave-forge/main/signal/inst/blackmanharris.m 2012-05-08 19:55:06 UTC (rev 10384) +++ trunk/octave-forge/main/signal/inst/blackmanharris.m 2012-05-09 01:25:06 UTC (rev 10385) @@ -31,6 +31,6 @@ a1 = 0.48829; a2 = 0.14128; a3 = 0.01168; - n = -ceil(N/2):N/2; - w = a0 + a1.*cos(2.*pi.*n./N) + a2.*cos(4.*pi.*n./N) + a3.*cos(6.*pi.*n./N); + n = 0:N; + w = a0 - a1.*cos(2.*pi.*n./N) + a2.*cos(4.*pi.*n./N) - a3.*cos(6.*pi.*n./N); endfunction This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <as...@us...> - 2012-05-08 19:55:12
|
Revision: 10384 http://octave.svn.sourceforge.net/octave/?rev=10384&view=rev Author: asnelt Date: 2012-05-08 19:55:06 +0000 (Tue, 08 May 2012) Log Message: ----------- Inserted missing comma. Modified Paths: -------------- trunk/octave-forge/main/statistics/inst/mnrnd.m trunk/octave-forge/main/statistics/inst/raylrnd.m Modified: trunk/octave-forge/main/statistics/inst/mnrnd.m =================================================================== --- trunk/octave-forge/main/statistics/inst/mnrnd.m 2012-05-08 14:56:54 UTC (rev 10383) +++ trunk/octave-forge/main/statistics/inst/mnrnd.m 2012-05-08 19:55:06 UTC (rev 10384) @@ -169,7 +169,7 @@ %! n = 10 * ones (3, 1); %! p = [0.2, 0.5, 0.3]; %! x = mnrnd (n, p); -%! assert (size (x), [length(n) length(p)]); +%! assert (size (x), [length(n), length(p)]); %! assert (all (x >= 0)); %! assert (all (round (x) == x)); %! assert (all (sum (x, 2) == n)); Modified: trunk/octave-forge/main/statistics/inst/raylrnd.m =================================================================== --- trunk/octave-forge/main/statistics/inst/raylrnd.m 2012-05-08 14:56:54 UTC (rev 10383) +++ trunk/octave-forge/main/statistics/inst/raylrnd.m 2012-05-08 19:55:06 UTC (rev 10384) @@ -153,5 +153,5 @@ %! r = 2; %! c = 3; %! x = raylrnd (sigma, r, c); -%! assert (size (x), [r c]); +%! assert (size (x), [r, c]); %! assert (all (x >= 0)); This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <nir...@us...> - 2012-05-08 14:57:05
|
Revision: 10383 http://octave.svn.sourceforge.net/octave/?rev=10383&view=rev Author: nir-krakauer Date: 2012-05-08 14:56:54 +0000 (Tue, 08 May 2012) Log Message: ----------- added cspas to INDEX Modified Paths: -------------- trunk/octave-forge/main/splines/INDEX Modified: trunk/octave-forge/main/splines/INDEX =================================================================== --- trunk/octave-forge/main/splines/INDEX 2012-05-08 13:19:19 UTC (rev 10382) +++ trunk/octave-forge/main/splines/INDEX 2012-05-08 14:56:54 UTC (rev 10383) @@ -2,6 +2,7 @@ Spline functions csapi csape + csaps fnplt fnder fnval This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2012-05-08 13:19:30
|
Revision: 10382 http://octave.svn.sourceforge.net/octave/?rev=10382&view=rev Author: paramaniac Date: 2012-05-08 13:19:19 +0000 (Tue, 08 May 2012) Log Message: ----------- control-devel: test example for mimo arx, pretty nonsensical Added Paths: ----------- trunk/octave-forge/extra/control-devel/devel/GlassFurnaceARX.m Added: trunk/octave-forge/extra/control-devel/devel/GlassFurnaceARX.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/GlassFurnaceARX.m (rev 0) +++ trunk/octave-forge/extra/control-devel/devel/GlassFurnaceARX.m 2012-05-08 13:19:19 UTC (rev 10382) @@ -0,0 +1,69 @@ +%{ +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 +sys = arx (dat, 5, [5 5 5]) + +%[y, t] = lsim (sys, U, [], x0); +[y, t] = lsim (sys(:, 1:3), U); + + +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') + + This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2012-05-08 13:13:49
|
Revision: 10381 http://octave.svn.sourceforge.net/octave/?rev=10381&view=rev Author: paramaniac Date: 2012-05-08 13:13:38 +0000 (Tue, 08 May 2012) Log Message: ----------- control-devel: mimo support for arx Modified Paths: -------------- trunk/octave-forge/extra/control-devel/inst/arx.m Modified: trunk/octave-forge/extra/control-devel/inst/arx.m =================================================================== --- trunk/octave-forge/extra/control-devel/inst/arx.m 2012-05-08 12:06:50 UTC (rev 10380) +++ trunk/octave-forge/extra/control-devel/inst/arx.m 2012-05-08 13:13:38 UTC (rev 10381) @@ -22,46 +22,96 @@ ## TODO: handle MIMO and multi-experiment data [~, p, m, ex] = size (dat); - Y = dat.y{1}; - U = dat.u{1}; + Y = dat.y; + U = dat.u; Ts = dat.tsam{1}; max_nb = max (nb); n = max (na, max_nb); + + num = cell (p, m+1); + den = cell (p, m+1); + + for i = 1 : p # for every output + Phi = cell (1, ex); + 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, :); + 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); + B = reshape (B, 1, []); + B = cellfun (@(x) [0; x], B, "uniformoutput", false); + + num(i, :) = [B, {1}]; + den(i, :) = repmat ({A}, 1, m+1); + endfor + + sys = filt (num, den, Ts); + +endfunction + + +function theta = __theta__ (phi, y, i, n) + + ## Theta = Phi \ Y(n+1:end, :); # naive formula + + ex = numel (phi); # number of experiments + + if (ex == 1) + ## single-experiment dataset + + ## solve linear least squares problem by pseudoinverse + ## the pseudoinverse is computed by singular value decomposition + ## M = U S V* ---> M+ = V S+ U* + ## Th = Ph \ Y = Ph+ Y + ## Th = V S+ U* Y, S+ = 1 ./ diag (S) + + [U, S, V] = svd (phi{1}, 0); # 0 for "economy size" decomposition, U overwrites input U + S = diag (S); # extract main diagonal + r = sum (S > eps*S(1)); + V = V(:, 1:r); + S = S(1:r); + U = U(:, 1:r); + theta = V * (S .\ (U' * y{1}(n+1:end, i))); # U' is the conjugate transpose + else + ## multi-experiment dataset + ## TODO: find more sophisticated formula than + ## Theta = (Phi1' Phi + Phi2' Phi2 + ...) \ (Phi1' Y1 + Phi2' Y2 + ...) + + ## covariance matrix C = (Phi1' Phi + Phi2' Phi2 + ...) + tmp = cellfun (@(Phi) Phi.' * Phi, phi, "uniformoutput", false); + C = plus (tmp{:}); + + ## PhiTY = (Phi1' Y1 + Phi2' Y2 + ...) + tmp = cellfun (@(Phi, Y) Phi.' * Y(n+1:end, i), phi, y, "uniformoutput", false); + PhiTY = plus (tmp{:}); + + ## pseudoinverse Theta = C \ Phi'Y + theta = C \ PhiTY; + 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, :); - - ## Theta = Phi \ Y(n+1:end, :); # naive formula - - ## solve linear least squares problem by pseudoinverse - ## the pseudoinverse is computed by singular value decomposition - ## M = U S V* ---> M+ = V S+ U* - ## Th = Ph \ Y = Ph+ Y - ## Th = V S+ U* Y, S+ = 1 ./ diag (S) - - [U, S, V] = svd (Phi, 0); # 0 for "economy size" decomposition, U overwrites input U - S = diag (S); # extract main diagonal - r = sum (S > eps*S(1)); - V = V(:, 1:r); - S = S(1:r); - U = U(:, 1:r); - Theta = V * (S .\ (U' * Y(n+1:end, :))); # U' is the conjugate transpose - - A = [1; Theta(1:na)]; # a0 = 1, a1 = Theta(1), an = Theta(n) - ## B = [0; Theta(na+1:end)]; # b0 = 0 (leading zero required by filt) - ThetaB = Theta(na+1:end); - B = mat2cell (ThetaB, nb); - B = reshape (B, 1, []); - B = cellfun (@(x) [0; x], B, "uniformoutput", false); - - ## sys = filt ({B, 1}, {A, A}, Ts); - num = [B, {1}]; - den = repmat ({A}, 1, m+1); - sys = filt (num, den, Ts); -endfunction \ No newline at end of file +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: <mtm...@us...> - 2012-05-08 12:07:02
|
Revision: 10380 http://octave.svn.sourceforge.net/octave/?rev=10380&view=rev Author: mtmiller Date: 2012-05-08 12:06:50 +0000 (Tue, 08 May 2012) Log Message: ----------- comm: simplify Makefile, remove dependency generation and related variables and targets Modified Paths: -------------- trunk/octave-forge/main/comm/src/Makefile Modified: trunk/octave-forge/main/comm/src/Makefile =================================================================== --- trunk/octave-forge/main/comm/src/Makefile 2012-05-08 10:47:23 UTC (rev 10379) +++ trunk/octave-forge/main/comm/src/Makefile 2012-05-08 12:06:50 UTC (rev 10380) @@ -7,64 +7,36 @@ op-gm-gm.cc op-gm-m.cc op-gm-s.cc op-m-gm.cc op-s-gm.cc \ ov-galois.cc GALOISOBJECTS = $(patsubst %.cc,%.o,$(GALOISSOURCES)) +GALOISHEADERS = galois.h galois-def.h galois-ops.h galoisfield.h ov-galois.h -GALOISDEPENDS = $(patsubst %.cc,%.d,$(GALOISSOURCES)) OTHERSOURCES = primpoly.cc isprimitive.cc __errcore__.cc cyclpoly.cc \ cyclgen.cc syndtable.cc __gfweight__.cc genqamdemod.cc OTHERTARGETS = $(patsubst %.cc,%.oct,$(OTHERSOURCES)) -OTHEROBJECTS = $(patsubst %.cc,%.o,$(OTHERSOURCES)) -OTHERDEPENDS = $(patsubst %.cc,%.d,$(OTHERSOURCES)) -TARGETS = $(GALOISTARGET) $(OTHERTARGETS) -SOURCES = $(GALOISSOURCES) $(OTHERSOURCES) -OBJECTS = $(GALOISOBJECTS) $(OTHEROBJECTS) -ifeq ($(MAKECMDGOALS),all) -DEPENDS = $(GALOISDEPENDS) $(OTHERDEPENDS) -endif -ifeq ($(MAKECMDGOALS),) -DEPENDS = $(GALOISDEPENDS) $(OTHERDEPENDS) -endif - -DELETES = $(OBJECTS) $(GALOISDEPENDS) $(OTHERDEPENDS) *~ $(TARGETS) core octave-core - DEFINES = -DGALOIS_DISP_PRIVATES MOFLAGS = .PHONY: all dist clean realclean count .SUFFIXES: -all : $(DEPENDS) $(OTHERTARGETS) $(GALOISTARGET) +all : $(OTHERTARGETS) $(GALOISTARGET) install : @$(INSTALL) -d $(DESTDIR)$(MPATH)/comm -$(GALOISTARGET) : $(DEPENDS) $(GALOISOBJECTS) - @echo "Linking $@"; \ +$(GALOISTARGET) : $(GALOISOBJECTS) $(MKOCTFILE) $(MOFLAGS) $(GALOISOBJECTS) -o $@ $(HDF5_LIBS) -ifneq (,$(DEPENDS)) - sinclude $(DEPENDS) -endif +$(GALOISOBJECTS): $(GALOISHEADERS) -%.oct : %.d %.o - @echo "Linking $@"; \ - $(MKOCTFILE) $(MOFLAGS) $(@:.oct=.o) -o $@ $(HDF5_LIBS) - -%.d: %.cc - @echo "Depending $<"; \ - $(MKOCTFILE) $(MOFLAGS) $(DEFINES) -M $< - %.o:%.cc - @echo "Compiling $@"; \ $(MKOCTFILE) $(MOFLAGS) $(DEFINES) -c $< clean: - @echo "Cleaning..."; \ - $(RM) -f $(DELETES) + $(RM) -f *.oct *.o core octave-core *~ -realclean: - @echo "Cleaning..."; \ - $(RM) -f $(DELETES) Makeconf config.log config.status +realclean: clean + $(RM) -f Makeconf config.log config.status dist: This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <cde...@us...> - 2012-05-08 10:47:34
|
Revision: 10379 http://octave.svn.sourceforge.net/octave/?rev=10379&view=rev Author: cdemills Date: 2012-05-08 10:47:23 +0000 (Tue, 08 May 2012) Log Message: ----------- Removes single and double quotes Modified Paths: -------------- trunk/octave-forge/extra/dataframe/inst/@dataframe/dataframe.m trunk/octave-forge/extra/dataframe/inst/@dataframe/private/df_matassign.m Modified: trunk/octave-forge/extra/dataframe/inst/@dataframe/dataframe.m =================================================================== --- trunk/octave-forge/extra/dataframe/inst/@dataframe/dataframe.m 2012-05-08 09:18:11 UTC (rev 10378) +++ trunk/octave-forge/extra/dataframe/inst/@dataframe/dataframe.m 2012-05-08 10:47:23 UTC (rev 10379) @@ -303,7 +303,7 @@ if (unquot) try %# remove quotes and leading space(s) - x(indj, indm) = regexp (dummy{indk}, '[^'' ].*[^'']', 'match'){1}; + x(indj, indm) = regexp (dummy{indk}, '[^''" ].*[^''"]', 'match'){1}; catch %# if the previous test fails, try a simpler one in = regexp (dummy{indk}, '[^'' ]+', 'match'); Modified: trunk/octave-forge/extra/dataframe/inst/@dataframe/private/df_matassign.m =================================================================== --- trunk/octave-forge/extra/dataframe/inst/@dataframe/private/df_matassign.m 2012-05-08 09:18:11 UTC (rev 10378) +++ trunk/octave-forge/extra/dataframe/inst/@dataframe/private/df_matassign.m 2012-05-08 10:47:23 UTC (rev 10379) @@ -67,17 +67,17 @@ df._ridx(indr, :, :) = []; %# to remove a line, iterate on each column df._data = cellfun (@(x) feval(@subsasgn, x, S, []), \ - df._data, "UniformOutPut", false); + df._data, "UniformOutPut", false); if (isa (indr, 'char')) - df._cnt(1) = 0; - else - df._cnt(1) = df._cnt(1) - length (indr); - endif + df._cnt(1) = 0; + else + df._cnt(1) = df._cnt(1) - length (indr); + endif endif df = df_thirddim (df); return; endif - + indc_was_set = ~isempty (indc); if (~indc_was_set) %# initial dataframe was empty ncol = size (RHS, 2); indc = 1:ncol; @@ -116,10 +116,10 @@ else inds = []; endif - + rname = cell(0, 0); rname_width = max (1, size (df._name{2}, 2)); ridx = []; cname = rname; ctype = rname; - + if (iscell (RHS)) if ((length (indc) == df._cnt(2) && size (RHS, 2) >= df._cnt(2)) \ || 0 == df._cnt(2) || isempty (S.subs{1}) || isempty (S.subs{2})) @@ -133,10 +133,10 @@ dummy = strcmp (dummy, 'char'); if (all (dummy)) if (length (df._over{2}) >= max (indc) \ - && ~all (df._over{2}(indc)) && ~isempty (S.subs{2})) + && ~all (df._over{2}(indc)) && ~isempty (S.subs{2})) warning("Trying to overwrite colum names"); endif - + cname = RHS(1, :).'; RHS = RHS(2:end, :); if (~indr_was_set) nrow = nrow - 1; indr = 1:nrow; @@ -157,19 +157,19 @@ %# at this stage, verify that the first line doesn't contain %# chars only; use them for column types dummy = cellfun ('class', \ - RHS(1, ~cellfun ('isempty', RHS(1, :))), \ - 'UniformOutput', false); + RHS(1, ~cellfun ('isempty', RHS(1, :))), \ + 'UniformOutput', false); dummy = strcmp (dummy, 'char'); if (all (dummy)) if (length (df._over{2}) >= max (indc) \ - && ~all (df._over{2}(indc))) + && ~all (df._over{2}(indc))) warning ("Trying to overwrite colum names"); endif if (sum (~cellfun ('isempty', RHS(1, indc))) == ncol) ctype = RHS(1, :); endif - + RHS = RHS(2:end, :); if (~indr_was_set) nrow = nrow - 1; indr = 1:nrow; @@ -181,7 +181,7 @@ %# row index and/or row name if (size (RHS, 1) > 1) dummy = all (cellfun ('isnumeric', \ - RHS(~cellfun ('isempty', RHS(:, 1)), 1))); + RHS(~cellfun ('isempty', RHS(:, 1)), 1))); else dummy = isnumeric(RHS{1, 1}); endif @@ -207,18 +207,18 @@ ridx = []; endif endif - + if (size (RHS, 2) > df._cnt(2)) %# verify the the first row doesn't contain chars only, use them %# for row names dummy = cellfun ('class', \ - RHS(~cellfun ('isempty', RHS(:, 1)), 1), \ - 'UniformOutput', false); + RHS(~cellfun ('isempty', RHS(:, 1)), 1), \ + 'UniformOutput', false); dummy = strcmp (dummy, 'char') \ && (~isempty (cname) && size (cname{1}, 2) < 1); if (all (dummy)) if (length (df._over{1}) >= max (indr) \ - && ~all (df._over{1}(indr))) + && ~all (df._over{1}(indr))) warning("Trying to overwrite row names"); else rname = RHS(:, 1); @@ -241,14 +241,14 @@ endif endif endif - + %# perform row resizing if columns are already filled if (~isempty (indr) && isnumeric(indr)) if (max (indr) > df._cnt(1) && size (df._data, 2) == df._cnt(2)) df = df_pad (df, 1, max (indr)-df._cnt(1), rname_width); endif endif - + if (iscell(RHS)) %# we must pad on a column-by-column basis %# verify that each cell contains a non-empty vector, and that sizes %# are compatible @@ -267,11 +267,15 @@ %# if 1 < size (RHS, 1) && any (dummy > 1), %# error("cells may only contain scalar"); %# endif - + if (size (RHS, 2) > indc) - keyboard + if (size (cname, 1) > indc) + ncol = size (RHS, 2); indc = 1:ncol; + else + keyboard + endif endif - + %# try to detect and remove bottom garbage eff_len = zeros (nrow, 1); if (size (RHS, 1) > 1) @@ -295,7 +299,7 @@ endwhile clear eff_len; endif - + %# the real assignement if (1 == size (RHS, 1)) %# each cell contains one vector fillfunc = @(x) RHS{x}; @@ -304,7 +308,7 @@ fillfunc = @(x) cell2mat (RHS(:, x)); endif - indj = 1; + indj = 1; for indi = (1:ncol) if (indc(indi) > df._cnt(2)) %# perform dynamic resizing one-by-one, to get type right This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2012-05-08 09:18:23
|
Revision: 10378 http://octave.svn.sourceforge.net/octave/?rev=10378&view=rev Author: paramaniac Date: 2012-05-08 09:18:11 +0000 (Tue, 08 May 2012) Log Message: ----------- control-devel: clean up example Modified Paths: -------------- trunk/octave-forge/extra/control-devel/devel/HeatingSystem.m trunk/octave-forge/extra/control-devel/devel/ident.m Modified: trunk/octave-forge/extra/control-devel/devel/HeatingSystem.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/HeatingSystem.m 2012-05-08 08:22:45 UTC (rev 10377) +++ trunk/octave-forge/extra/control-devel/devel/HeatingSystem.m 2012-05-08 09:18:11 UTC (rev 10378) @@ -75,29 +75,22 @@ 'outunit', '°C') % s=15, n=7 -sys = arx (dat, 7, 7) % normally na = nb +[sys1, x0] = ident (dat, 15, 7) +sys2 = arx (dat, 7, 7) % normally na = nb -[y, t] = lsim (sys(:, 1), U); +[y1, t1] = lsim (sys1, U, [], x0); +[y2, t] = lsim (sys2(:, 1), U); +err1 = norm (Y - y1, 1) / norm (Y, 1) +err2 = norm (Y - y2, 1) / norm (Y, 1) -err = norm (Y - y, 1) / norm (Y, 1) - figure (1) -plot (t, Y(:,1), 'b', t, y(:,1), 'r') +plot (t, Y, t, y1, t, y2) title ('DaISy: Heating System [99-001]') -legend ('measured temperature', 'simulated temperature', 'location', 'southeast') +xlim ([t(1), t(end)]) +xlabel ('Time [s]') +ylabel ('Temperature [Degree Celsius]') +legend ('measured', 'simulated subspace', 'simulated ARX', 'location', 'southeast') -[sys2, x0] = ident (dat, 15, 7) -[y2, t2] = lsim (sys2, U, [], x0); - -err2 = norm (Y - y2, 1) / norm (Y, 1) - -figure (2) -plot (t, Y, 'b', t, y, 'r', t, y2, 'g') -title ('DaISy: Heating System [99-001]') -legend ('measured temperature', 'simulated temperature', 'slicot', 'location', 'southeast') - - - Modified: trunk/octave-forge/extra/control-devel/devel/ident.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/ident.m 2012-05-08 08:22:45 UTC (rev 10377) +++ trunk/octave-forge/extra/control-devel/devel/ident.m 2012-05-08 09:18:11 UTC (rev 10378) @@ -49,6 +49,6 @@ %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, -1); + sys = ss (a, b, c, d, dat.tsam{1}); endfunction This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2012-05-08 08:22:52
|
Revision: 10377 http://octave.svn.sourceforge.net/octave/?rev=10377&view=rev Author: paramaniac Date: 2012-05-08 08:22:45 +0000 (Tue, 08 May 2012) Log Message: ----------- control-devel: add draft of another DaISy example Modified Paths: -------------- trunk/octave-forge/extra/control-devel/devel/pHarx.m Added Paths: ----------- trunk/octave-forge/extra/control-devel/devel/HeatingSystem.m trunk/octave-forge/extra/control-devel/devel/heating_system.dat Added: trunk/octave-forge/extra/control-devel/devel/HeatingSystem.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/HeatingSystem.m (rev 0) +++ trunk/octave-forge/extra/control-devel/devel/HeatingSystem.m 2012-05-08 08:22:45 UTC (rev 10377) @@ -0,0 +1,103 @@ +%{ +1. Contributed by: + + Roy Smith + Dept. of Electrical & Computer Engineering + University of California, + Santa Barbara, CA 93106 + U.S.A. + ro...@ec... + +2. Process/Description: + + The experiment is a simple SISO heating system. + The input drives a 300 Watt Halogen lamp, suspended + several inches above a thin steel plate. The output + is a thermocouple measurement taken from the back of + the plate. + +3. Sampling interval: + + 2.0 seconds + +4. Number of samples + + 801 + +5. Inputs: + + u: input drive voltage + ... +6. Outputs: + + y: temperature (deg. C) + ... +7. References: + + The use of this experiment and data for robust + control model validation is described in: + + "Sampled Data Model Validation: an Algorithm and + Experimental Application," Geir Dullerud & Roy Smith, + International Journal of Robust and Nonlinear Control, + Vol. 6, No. 9/10, pp. 1065-1078, 1996. + +8. Known properties/peculiarities + + The data (and nominal model) is the above paper have the + output expressed in 10's deg. C. This has been rescaled + to the original units of deg. C. in the DaISy data set. + There is also a -1 volt offset in u in the data shown plotted + in the original paper. This has been removed in the + DaISy dataset. + + The data shows evidence of discrepancies. One of the + issues studied in the above paper is the size of these + discrepancies - measured in this case in terms of the norm + of the smallest perturbation required to account for the + difference between the nominal model and the data. + + The steady state input (prior to the start of the experiment) + is u = 6.0 Volts. + +%} + +clear all, close all, clc + +load heating_system.dat +U=heating_system(:,2); +Y=heating_system(:,3); + + +dat = iddata (Y, U, 2.0, 'inname', 'input drive voltage', \ + 'inunit', 'Volt', \ + 'outname', 'temperature', \ + 'outunit', '°C') + +% s=15, n=7 +sys = arx (dat, 7, 7) % normally na = nb + +[y, t] = lsim (sys(:, 1), U); + + +err = norm (Y - y, 1) / norm (Y, 1) + +figure (1) +plot (t, Y(:,1), 'b', t, y(:,1), 'r') +title ('DaISy: Heating System [99-001]') +legend ('measured temperature', 'simulated temperature', 'location', 'southeast') + + +[sys2, x0] = ident (dat, 15, 7) + +[y2, t2] = lsim (sys2, U, [], x0); + +err2 = norm (Y - y2, 1) / norm (Y, 1) + +figure (2) +plot (t, Y, 'b', t, y, 'r', t, y2, 'g') +title ('DaISy: Heating System [99-001]') +legend ('measured temperature', 'simulated temperature', 'slicot', 'location', 'southeast') + + + Added: trunk/octave-forge/extra/control-devel/devel/heating_system.dat =================================================================== --- trunk/octave-forge/extra/control-devel/devel/heating_system.dat (rev 0) +++ trunk/octave-forge/extra/control-devel/devel/heating_system.dat 2012-05-08 08:22:45 UTC (rev 10377) @@ -0,0 +1,801 @@ + 0.00 9.000 129.320 + 2.00 9.000 129.320 + 4.00 9.000 129.564 + 6.00 9.000 131.516 + 8.00 9.000 134.200 + 10.00 9.000 136.640 + 12.00 9.000 139.324 + 14.00 9.000 141.764 + 16.00 9.000 144.204 + 18.00 9.000 146.156 + 20.00 9.000 148.596 + 22.00 9.000 150.548 + 24.00 9.000 152.256 + 26.00 9.000 154.208 + 28.00 9.000 155.428 + 30.00 9.000 156.648 + 32.00 9.000 157.624 + 34.00 9.000 158.600 + 36.00 9.000 159.576 + 38.00 9.000 160.796 + 40.00 9.000 161.772 + 42.00 9.000 162.260 + 44.00 9.000 162.748 + 46.00 9.000 163.724 + 48.00 9.000 164.212 + 50.00 9.000 165.188 + 52.00 9.000 166.164 + 54.00 9.000 167.140 + 56.00 9.000 167.872 + 58.00 9.000 168.604 + 60.00 9.000 169.092 + 62.00 9.000 169.336 + 64.00 9.000 169.824 + 66.00 9.000 170.312 + 68.00 9.000 170.556 + 70.00 9.000 171.044 + 72.00 9.000 171.532 + 74.00 9.000 172.020 + 76.00 9.000 172.508 + 78.00 9.000 172.752 + 80.00 9.000 172.996 + 82.00 9.000 173.240 + 84.00 9.000 173.484 + 86.00 9.000 173.972 + 88.00 9.000 173.728 + 90.00 9.000 173.728 + 92.00 9.000 174.216 + 94.00 9.000 174.948 + 96.00 9.000 175.680 + 98.00 9.000 175.924 + 100.00 9.000 175.924 + 102.00 9.000 176.168 + 104.00 9.000 176.900 + 106.00 9.000 177.144 + 108.00 9.000 177.632 + 110.00 9.000 178.120 + 112.00 9.000 177.876 + 114.00 9.000 178.364 + 116.00 9.000 178.364 + 118.00 9.000 178.364 + 120.00 9.000 178.852 + 122.00 9.000 179.096 + 124.00 9.000 179.584 + 126.00 9.000 179.828 + 128.00 9.000 180.072 + 130.00 9.000 180.316 + 132.00 9.000 180.316 + 134.00 9.000 180.804 + 136.00 9.000 180.804 + 138.00 9.000 181.292 + 140.00 9.000 181.536 + 142.00 9.000 182.024 + 144.00 9.000 182.024 + 146.00 9.000 182.024 + 148.00 9.000 182.024 + 150.00 9.000 182.268 + 152.00 9.000 182.512 + 154.00 9.000 182.756 + 156.00 9.000 183.000 + 158.00 9.000 183.000 + 160.00 9.000 183.000 + 162.00 9.000 183.488 + 164.00 9.000 183.732 + 166.00 9.000 183.976 + 168.00 9.000 184.220 + 170.00 9.000 184.220 + 172.00 9.000 184.464 + 174.00 9.000 184.708 + 176.00 9.000 184.952 + 178.00 9.000 184.464 + 180.00 9.000 184.220 + 182.00 9.000 184.464 + 184.00 9.000 184.708 + 186.00 9.000 185.196 + 188.00 9.000 185.440 + 190.00 9.000 185.684 + 192.00 9.000 185.440 + 194.00 9.000 185.684 + 196.00 9.000 185.928 + 198.00 9.000 186.172 + 200.00 9.000 186.172 + 202.00 9.000 186.416 + 204.00 9.000 186.416 + 206.00 9.000 186.416 + 208.00 9.000 186.416 + 210.00 9.000 186.172 + 212.00 9.000 186.172 + 214.00 9.000 186.416 + 216.00 9.000 186.904 + 218.00 9.000 187.392 + 220.00 9.000 187.636 + 222.00 9.000 187.636 + 224.00 9.000 187.880 + 226.00 9.000 187.880 + 228.00 9.000 187.880 + 230.00 9.000 187.880 + 232.00 9.000 187.880 + 234.00 9.000 187.636 + 236.00 9.000 187.880 + 238.00 9.000 187.880 + 240.00 9.000 187.880 + 242.00 9.000 187.880 + 244.00 9.000 188.124 + 246.00 9.000 188.368 + 248.00 9.000 188.368 + 250.00 9.000 188.368 + 252.00 9.000 188.612 + 254.00 9.000 187.880 + 256.00 9.000 187.880 + 258.00 9.000 187.880 + 260.00 9.000 188.124 + 262.00 9.000 188.124 + 264.00 9.000 188.368 + 266.00 9.000 189.100 + 268.00 9.000 189.344 + 270.00 9.000 189.344 + 272.00 9.000 189.344 + 274.00 9.000 189.588 + 276.00 9.000 189.100 + 278.00 9.000 189.344 + 280.00 9.000 189.832 + 282.00 9.000 190.320 + 284.00 9.000 190.076 + 286.00 9.000 189.832 + 288.00 9.000 189.832 + 290.00 9.000 190.320 + 292.00 9.000 190.564 + 294.00 9.000 190.564 + 296.00 9.000 190.320 + 298.00 9.000 189.588 + 300.00 9.000 189.588 + 302.00 9.000 190.076 + 304.00 9.000 190.076 + 306.00 9.000 190.076 + 308.00 9.000 190.076 + 310.00 9.000 190.320 + 312.00 9.000 190.076 + 314.00 9.000 190.320 + 316.00 9.000 190.320 + 318.00 9.000 190.564 + 320.00 9.000 190.564 + 322.00 9.000 190.564 + 324.00 9.000 190.320 + 326.00 9.000 190.320 + 328.00 9.000 190.320 + 330.00 9.000 190.320 + 332.00 9.000 190.564 + 334.00 9.000 191.052 + 336.00 9.000 191.540 + 338.00 9.000 191.296 + 340.00 9.000 191.296 + 342.00 9.000 191.296 + 344.00 9.000 191.296 + 346.00 9.000 191.296 + 348.00 9.000 191.784 + 350.00 9.000 191.296 + 352.00 9.000 191.296 + 354.00 9.000 191.540 + 356.00 9.000 192.028 + 358.00 9.000 191.540 + 360.00 9.000 191.784 + 362.00 9.000 191.784 + 364.00 9.000 191.540 + 366.00 9.000 191.052 + 368.00 9.000 191.052 + 370.00 9.000 191.052 + 372.00 9.000 190.808 + 374.00 9.000 190.808 + 376.00 9.000 191.052 + 378.00 9.000 191.052 + 380.00 9.000 191.052 + 382.00 9.000 191.296 + 384.00 9.000 191.540 + 386.00 9.000 191.540 + 388.00 9.000 191.540 + 390.00 9.000 191.784 + 392.00 9.000 192.028 + 394.00 9.000 192.028 + 396.00 9.000 192.028 + 398.00 9.000 192.028 + 400.00 4.000 192.028 + 402.00 4.000 191.540 + 404.00 4.000 189.588 + 406.00 4.000 186.172 + 408.00 4.000 181.780 + 410.00 4.000 177.632 + 412.00 4.000 173.972 + 414.00 4.000 170.068 + 416.00 4.000 166.164 + 418.00 4.000 162.992 + 420.00 4.000 160.064 + 422.00 4.000 157.624 + 424.00 4.000 155.184 + 426.00 4.000 153.232 + 428.00 4.000 151.036 + 430.00 4.000 149.328 + 432.00 4.000 147.132 + 434.00 4.000 145.180 + 436.00 4.000 143.228 + 438.00 4.000 142.008 + 440.00 4.000 140.544 + 442.00 4.000 139.324 + 444.00 4.000 138.348 + 446.00 4.000 137.128 + 448.00 4.000 136.152 + 450.00 4.000 135.420 + 452.00 4.000 134.444 + 454.00 4.000 133.468 + 456.00 4.000 132.736 + 458.00 4.000 131.516 + 460.00 4.000 130.296 + 462.00 4.000 129.320 + 464.00 4.000 128.588 + 466.00 4.000 127.856 + 468.00 4.000 127.124 + 470.00 4.000 126.392 + 472.00 4.000 125.416 + 474.00 4.000 124.684 + 476.00 4.000 123.952 + 478.00 4.000 123.464 + 480.00 4.000 122.732 + 482.00 4.000 122.000 + 484.00 4.000 121.512 + 486.00 4.000 120.780 + 488.00 4.000 120.536 + 490.00 4.000 119.804 + 492.00 4.000 119.316 + 494.00 4.000 118.584 + 496.00 4.000 117.852 + 498.00 4.000 117.364 + 500.00 4.000 116.876 + 502.00 4.000 116.388 + 504.00 4.000 116.144 + 506.00 4.000 115.656 + 508.00 4.000 115.412 + 510.00 4.000 114.680 + 512.00 4.000 114.436 + 514.00 4.000 113.948 + 516.00 4.000 113.704 + 518.00 4.000 112.972 + 520.00 4.000 112.728 + 522.00 4.000 112.240 + 524.00 4.000 111.996 + 526.00 4.000 111.752 + 528.00 4.000 111.752 + 530.00 4.000 111.752 + 532.00 4.000 111.508 + 534.00 4.000 111.264 + 536.00 4.000 111.020 + 538.00 4.000 111.264 + 540.00 4.000 110.776 + 542.00 4.000 110.288 + 544.00 4.000 110.044 + 546.00 4.000 109.800 + 548.00 4.000 109.312 + 550.00 4.000 109.556 + 552.00 4.000 109.556 + 554.00 4.000 109.556 + 556.00 4.000 109.312 + 558.00 4.000 108.824 + 560.00 4.000 108.092 + 562.00 4.000 107.360 + 564.00 4.000 106.872 + 566.00 4.000 106.628 + 568.00 4.000 106.140 + 570.00 4.000 105.896 + 572.00 4.000 105.652 + 574.00 4.000 105.408 + 576.00 4.000 104.920 + 578.00 4.000 104.676 + 580.00 4.000 104.188 + 582.00 4.000 103.700 + 584.00 4.000 103.456 + 586.00 4.000 103.212 + 588.00 4.000 102.968 + 590.00 4.000 102.480 + 592.00 4.000 102.236 + 594.00 4.000 101.748 + 596.00 4.000 101.504 + 598.00 4.000 101.260 + 600.00 4.000 101.260 + 602.00 4.000 101.016 + 604.00 4.000 100.772 + 606.00 4.000 100.528 + 608.00 4.000 100.528 + 610.00 4.000 100.284 + 612.00 4.000 99.796 + 614.00 4.000 99.308 + 616.00 4.000 99.308 + 618.00 4.000 99.308 + 620.00 4.000 99.064 + 622.00 4.000 98.820 + 624.00 4.000 98.576 + 626.00 4.000 98.576 + 628.00 4.000 98.088 + 630.00 4.000 97.844 + 632.00 4.000 97.600 + 634.00 4.000 97.112 + 636.00 4.000 96.868 + 638.00 4.000 96.624 + 640.00 4.000 96.868 + 642.00 4.000 96.624 + 644.00 4.000 96.624 + 646.00 4.000 96.380 + 648.00 4.000 95.892 + 650.00 4.000 95.892 + 652.00 4.000 95.892 + 654.00 4.000 95.404 + 656.00 4.000 95.648 + 658.00 4.000 95.404 + 660.00 4.000 95.160 + 662.00 4.000 94.916 + 664.00 4.000 94.672 + 666.00 4.000 94.672 + 668.00 4.000 94.428 + 670.00 4.000 94.428 + 672.00 4.000 94.184 + 674.00 4.000 93.940 + 676.00 4.000 93.696 + 678.00 4.000 93.696 + 680.00 4.000 93.696 + 682.00 4.000 93.208 + 684.00 4.000 92.964 + 686.00 4.000 92.720 + 688.00 4.000 92.720 + 690.00 4.000 92.476 + 692.00 4.000 92.476 + 694.00 4.000 92.476 + 696.00 4.000 92.476 + 698.00 4.000 92.476 + 700.00 4.000 92.476 + 702.00 4.000 91.988 + 704.00 4.000 91.988 + 706.00 4.000 91.744 + 708.00 4.000 91.744 + 710.00 4.000 91.500 + 712.00 4.000 91.012 + 714.00 4.000 91.012 + 716.00 4.000 90.768 + 718.00 4.000 90.768 + 720.00 4.000 90.768 + 722.00 4.000 90.768 + 724.00 4.000 91.012 + 726.00 4.000 90.768 + 728.00 4.000 90.524 + 730.00 4.000 90.280 + 732.00 4.000 90.280 + 734.00 4.000 90.036 + 736.00 4.000 89.792 + 738.00 4.000 89.548 + 740.00 4.000 89.548 + 742.00 4.000 89.548 + 744.00 4.000 89.304 + 746.00 4.000 89.304 + 748.00 4.000 89.060 + 750.00 4.000 89.060 + 752.00 4.000 88.816 + 754.00 4.000 88.572 + 756.00 4.000 88.572 + 758.00 4.000 88.816 + 760.00 4.000 88.572 + 762.00 4.000 88.328 + 764.00 4.000 88.084 + 766.00 4.000 87.840 + 768.00 4.000 87.352 + 770.00 4.000 87.596 + 772.00 4.000 87.352 + 774.00 4.000 87.352 + 776.00 4.000 87.352 + 778.00 4.000 87.108 + 780.00 4.000 87.108 + 782.00 4.000 87.108 + 784.00 4.000 87.108 + 786.00 4.000 86.864 + 788.00 4.000 87.108 + 790.00 4.000 86.864 + 792.00 4.000 86.620 + 794.00 4.000 86.620 + 796.00 4.000 86.620 + 798.00 4.000 86.376 + 800.00 8.000 86.132 + 802.00 8.000 85.888 + 804.00 8.000 86.864 + 806.00 8.000 89.548 + 808.00 8.000 93.208 + 810.00 8.000 97.356 + 812.00 8.000 101.504 + 814.00 8.000 105.164 + 816.00 8.000 108.824 + 818.00 8.000 111.996 + 820.00 8.000 114.924 + 822.00 8.000 117.364 + 824.00 8.000 119.316 + 826.00 8.000 121.512 + 828.00 8.000 123.464 + 830.00 8.000 124.928 + 832.00 8.000 126.636 + 834.00 8.000 128.100 + 836.00 8.000 129.564 + 838.00 8.000 130.784 + 840.00 8.000 132.004 + 842.00 8.000 133.224 + 844.00 8.000 134.200 + 846.00 8.000 135.176 + 848.00 8.000 136.152 + 850.00 8.000 136.884 + 852.00 8.000 138.104 + 854.00 8.000 138.836 + 856.00 8.000 139.324 + 858.00 8.000 140.300 + 860.00 8.000 141.276 + 862.00 8.000 141.764 + 864.00 8.000 142.008 + 866.00 8.000 142.740 + 868.00 8.000 143.960 + 870.00 8.000 144.448 + 872.00 8.000 145.180 + 874.00 8.000 145.912 + 876.00 8.000 146.644 + 878.00 8.000 147.376 + 880.00 8.000 147.620 + 882.00 8.000 148.352 + 884.00 8.000 148.840 + 886.00 8.000 149.084 + 888.00 8.000 149.572 + 890.00 8.000 150.060 + 892.00 8.000 150.304 + 894.00 8.000 150.548 + 896.00 8.000 151.036 + 898.00 8.000 151.524 + 900.00 8.000 152.500 + 902.00 8.000 152.988 + 904.00 8.000 153.232 + 906.00 8.000 153.720 + 908.00 8.000 154.208 + 910.00 8.000 154.696 + 912.00 8.000 155.184 + 914.00 8.000 155.672 + 916.00 8.000 155.672 + 918.00 8.000 155.672 + 920.00 8.000 155.916 + 922.00 8.000 156.160 + 924.00 8.000 156.404 + 926.00 8.000 156.404 + 928.00 8.000 156.892 + 930.00 8.000 157.380 + 932.00 8.000 157.624 + 934.00 8.000 158.112 + 936.00 8.000 158.600 + 938.00 8.000 158.600 + 940.00 8.000 159.088 + 942.00 8.000 159.576 + 944.00 8.000 159.332 + 946.00 8.000 159.332 + 948.00 8.000 159.576 + 950.00 8.000 159.820 + 952.00 8.000 159.820 + 954.00 8.000 159.820 + 956.00 8.000 160.064 + 958.00 8.000 160.308 + 960.00 8.000 160.552 + 962.00 8.000 160.796 + 964.00 8.000 160.796 + 966.00 8.000 160.796 + 968.00 8.000 161.040 + 970.00 8.000 161.528 + 972.00 8.000 161.772 + 974.00 8.000 161.772 + 976.00 8.000 162.016 + 978.00 8.000 162.504 + 980.00 8.000 162.504 + 982.00 8.000 162.748 + 984.00 8.000 162.748 + 986.00 8.000 162.748 + 988.00 8.000 162.992 + 990.00 8.000 163.236 + 992.00 8.000 163.236 + 994.00 8.000 162.992 + 996.00 8.000 163.236 + 998.00 8.000 163.480 + 1000.00 3.000 163.724 + 1002.00 3.000 163.968 + 1004.00 3.000 162.992 + 1006.00 3.000 160.308 + 1008.00 3.000 155.916 + 1010.00 3.000 151.280 + 1012.00 3.000 146.644 + 1014.00 3.000 142.740 + 1016.00 3.000 139.080 + 1018.00 3.000 136.152 + 1020.00 3.000 133.224 + 1022.00 3.000 130.784 + 1024.00 3.000 128.588 + 1026.00 3.000 126.636 + 1028.00 3.000 125.172 + 1030.00 3.000 123.220 + 1032.00 3.000 121.756 + 1034.00 3.000 120.292 + 1036.00 3.000 118.828 + 1038.00 3.000 117.364 + 1040.00 3.000 115.900 + 1042.00 3.000 114.680 + 1044.00 3.000 113.460 + 1046.00 3.000 112.240 + 1048.00 3.000 111.264 + 1050.00 3.000 110.288 + 1052.00 3.000 109.068 + 1054.00 3.000 107.604 + 1056.00 3.000 106.628 + 1058.00 3.000 105.652 + 1060.00 3.000 104.676 + 1062.00 3.000 103.944 + 1064.00 3.000 102.968 + 1066.00 3.000 102.236 + 1068.00 3.000 101.504 + 1070.00 3.000 100.284 + 1072.00 3.000 99.796 + 1074.00 3.000 99.308 + 1076.00 3.000 98.332 + 1078.00 3.000 97.844 + 1080.00 3.000 97.356 + 1082.00 3.000 96.868 + 1084.00 3.000 96.380 + 1086.00 3.000 95.892 + 1088.00 3.000 95.160 + 1090.00 3.000 94.672 + 1092.00 3.000 94.184 + 1094.00 3.000 93.696 + 1096.00 3.000 93.208 + 1098.00 3.000 92.964 + 1100.00 3.000 92.476 + 1102.00 3.000 91.988 + 1104.00 3.000 91.256 + 1106.00 3.000 90.768 + 1108.00 3.000 90.280 + 1110.00 3.000 89.792 + 1112.00 3.000 89.548 + 1114.00 3.000 89.060 + 1116.00 3.000 88.084 + 1118.00 3.000 87.596 + 1120.00 3.000 87.352 + 1122.00 3.000 87.108 + 1124.00 3.000 86.864 + 1126.00 3.000 86.864 + 1128.00 3.000 86.376 + 1130.00 3.000 85.888 + 1132.00 3.000 85.888 + 1134.00 3.000 85.644 + 1136.00 3.000 84.912 + 1138.00 3.000 84.180 + 1140.00 3.000 84.180 + 1142.00 3.000 83.936 + 1144.00 3.000 83.692 + 1146.00 3.000 83.448 + 1148.00 3.000 83.204 + 1150.00 3.000 82.716 + 1152.00 3.000 82.228 + 1154.00 3.000 81.984 + 1156.00 3.000 81.984 + 1158.00 3.000 81.740 + 1160.00 3.000 81.740 + 1162.00 3.000 81.496 + 1164.00 3.000 81.008 + 1166.00 3.000 80.764 + 1168.00 3.000 80.276 + 1170.00 3.000 80.032 + 1172.00 3.000 79.788 + 1174.00 3.000 79.544 + 1176.00 3.000 79.300 + 1178.00 3.000 79.056 + 1180.00 3.000 78.568 + 1182.00 3.000 78.568 + 1184.00 3.000 78.324 + 1186.00 3.000 78.324 + 1188.00 3.000 78.080 + 1190.00 3.000 77.836 + 1192.00 3.000 77.836 + 1194.00 3.000 77.592 + 1196.00 3.000 77.104 + 1198.00 3.000 76.860 + 1200.00 7.000 76.372 + 1202.00 7.000 76.372 + 1204.00 7.000 77.104 + 1206.00 7.000 79.300 + 1208.00 7.000 82.472 + 1210.00 7.000 85.888 + 1212.00 7.000 89.060 + 1214.00 7.000 92.232 + 1216.00 7.000 95.160 + 1218.00 7.000 97.600 + 1220.00 7.000 99.552 + 1222.00 7.000 101.260 + 1224.00 7.000 103.212 + 1226.00 7.000 105.164 + 1228.00 7.000 106.628 + 1230.00 7.000 107.604 + 1232.00 7.000 108.580 + 1234.00 7.000 109.800 + 1236.00 7.000 110.776 + 1238.00 7.000 111.996 + 1240.00 7.000 112.972 + 1242.00 7.000 113.948 + 1244.00 7.000 114.680 + 1246.00 7.000 115.656 + 1248.00 7.000 116.144 + 1250.00 7.000 117.120 + 1252.00 7.000 117.608 + 1254.00 7.000 118.340 + 1256.00 7.000 118.828 + 1258.00 7.000 119.804 + 1260.00 7.000 120.292 + 1262.00 7.000 121.024 + 1264.00 7.000 121.512 + 1266.00 7.000 122.000 + 1268.00 7.000 122.244 + 1270.00 7.000 122.976 + 1272.00 7.000 123.464 + 1274.00 7.000 123.952 + 1276.00 7.000 124.440 + 1278.00 7.000 124.928 + 1280.00 7.000 125.660 + 1282.00 7.000 125.904 + 1284.00 7.000 126.148 + 1286.00 7.000 126.636 + 1288.00 7.000 126.636 + 1290.00 7.000 127.124 + 1292.00 7.000 127.368 + 1294.00 7.000 127.368 + 1296.00 7.000 127.368 + 1298.00 7.000 127.368 + 1300.00 7.000 127.856 + 1302.00 7.000 128.344 + 1304.00 7.000 128.588 + 1306.00 7.000 129.076 + 1308.00 7.000 129.320 + 1310.00 7.000 129.564 + 1312.00 7.000 130.052 + 1314.00 7.000 130.296 + 1316.00 7.000 130.296 + 1318.00 7.000 130.540 + 1320.00 7.000 130.784 + 1322.00 7.000 131.272 + 1324.00 7.000 131.516 + 1326.00 7.000 131.760 + 1328.00 7.000 131.760 + 1330.00 7.000 132.004 + 1332.00 7.000 132.004 + 1334.00 7.000 132.492 + 1336.00 7.000 132.492 + 1338.00 7.000 132.492 + 1340.00 7.000 132.248 + 1342.00 7.000 132.492 + 1344.00 7.000 132.980 + 1346.00 7.000 132.980 + 1348.00 7.000 132.980 + 1350.00 7.000 133.468 + 1352.00 7.000 133.956 + 1354.00 7.000 133.956 + 1356.00 7.000 133.956 + 1358.00 7.000 133.956 + 1360.00 7.000 134.200 + 1362.00 7.000 134.444 + 1364.00 7.000 134.932 + 1366.00 7.000 135.176 + 1368.00 7.000 135.176 + 1370.00 7.000 135.176 + 1372.00 7.000 135.176 + 1374.00 7.000 135.420 + 1376.00 7.000 135.664 + 1378.00 7.000 135.908 + 1380.00 7.000 136.396 + 1382.00 7.000 136.396 + 1384.00 7.000 136.396 + 1386.00 7.000 136.640 + 1388.00 7.000 136.396 + 1390.00 7.000 136.884 + 1392.00 7.000 136.640 + 1394.00 7.000 136.884 + 1396.00 7.000 136.640 + 1398.00 7.000 136.640 + 1400.00 6.000 136.884 + 1402.00 6.000 136.884 + 1404.00 6.000 136.884 + 1406.00 6.000 136.396 + 1408.00 6.000 135.664 + 1410.00 6.000 134.932 + 1412.00 6.000 133.956 + 1414.00 6.000 132.980 + 1416.00 6.000 132.248 + 1418.00 6.000 131.760 + 1420.00 6.000 131.028 + 1422.00 6.000 130.784 + 1424.00 6.000 130.296 + 1426.00 6.000 130.052 + 1428.00 6.000 129.320 + 1430.00 6.000 128.588 + 1432.00 6.000 128.100 + 1434.00 6.000 127.612 + 1436.00 6.000 127.368 + 1438.00 6.000 127.368 + 1440.00 6.000 127.368 + 1442.00 6.000 127.368 + 1444.00 6.000 127.124 + 1446.00 6.000 127.124 + 1448.00 6.000 126.636 + 1450.00 6.000 126.636 + 1452.00 6.000 126.392 + 1454.00 6.000 126.148 + 1456.00 6.000 126.148 + 1458.00 6.000 126.148 + 1460.00 6.000 125.904 + 1462.00 6.000 125.904 + 1464.00 6.000 125.416 + 1466.00 6.000 125.172 + 1468.00 6.000 124.928 + 1470.00 6.000 124.928 + 1472.00 6.000 124.928 + 1474.00 6.000 125.172 + 1476.00 6.000 125.660 + 1478.00 6.000 126.148 + 1480.00 6.000 126.636 + 1482.00 6.000 126.636 + 1484.00 6.000 126.636 + 1486.00 6.000 126.392 + 1488.00 6.000 126.636 + 1490.00 6.000 126.880 + 1492.00 6.000 126.636 + 1494.00 6.000 125.904 + 1496.00 6.000 125.416 + 1498.00 6.000 125.416 + 1500.00 6.000 126.148 + 1502.00 6.000 126.880 + 1504.00 6.000 126.880 + 1506.00 6.000 126.880 + 1508.00 6.000 126.880 + 1510.00 6.000 126.636 + 1512.00 6.000 126.392 + 1514.00 6.000 126.148 + 1516.00 6.000 125.904 + 1518.00 6.000 125.660 + 1520.00 6.000 125.416 + 1522.00 6.000 125.172 + 1524.00 6.000 124.928 + 1526.00 6.000 124.440 + 1528.00 6.000 123.952 + 1530.00 6.000 123.708 + 1532.00 6.000 123.220 + 1534.00 6.000 122.976 + 1536.00 6.000 122.244 + 1538.00 6.000 122.000 + 1540.00 6.000 122.488 + 1542.00 6.000 122.732 + 1544.00 6.000 122.732 + 1546.00 6.000 122.732 + 1548.00 6.000 122.732 + 1550.00 6.000 122.244 + 1552.00 6.000 122.244 + 1554.00 6.000 122.000 + 1556.00 6.000 121.756 + 1558.00 6.000 121.756 + 1560.00 6.000 121.756 + 1562.00 6.000 121.756 + 1564.00 6.000 121.756 + 1566.00 6.000 121.268 + 1568.00 6.000 121.268 + 1570.00 6.000 121.268 + 1572.00 6.000 121.024 + 1574.00 6.000 120.780 + 1576.00 6.000 120.780 + 1578.00 6.000 120.536 + 1580.00 6.000 120.780 + 1582.00 6.000 120.780 + 1584.00 6.000 121.024 + 1586.00 6.000 121.024 + 1588.00 6.000 121.268 + 1590.00 6.000 121.268 + 1592.00 6.000 120.780 + 1594.00 6.000 120.780 + 1596.00 6.000 120.536 + 1598.00 6.000 120.536 + 1600.00 6.000 120.536 Modified: trunk/octave-forge/extra/control-devel/devel/pHarx.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/pHarx.m 2012-05-08 07:56:26 UTC (rev 10376) +++ trunk/octave-forge/extra/control-devel/devel/pHarx.m 2012-05-08 08:22:45 UTC (rev 10377) @@ -52,7 +52,7 @@ dat = iddata (Y, U) % [sys, x0] = ident (dat, 15, 6) % s=15, n=6 -sys = arx (dat, 6, [3,3]) +sys = arx (dat, 6, [6,6]) % normally na = nb % [y, t] = lsim (sys, U, [], x0); [y, t] = lsim (sys(:, 1:2), U); This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2012-05-08 07:56:33
|
Revision: 10376 http://octave.svn.sourceforge.net/octave/?rev=10376&view=rev Author: paramaniac Date: 2012-05-08 07:56:26 +0000 (Tue, 08 May 2012) Log Message: ----------- control: list tf property 'inv', require discrete-time systems Modified Paths: -------------- trunk/octave-forge/main/control/inst/@tf/__property_names__.m trunk/octave-forge/main/control/inst/@tf/__set__.m Modified: trunk/octave-forge/main/control/inst/@tf/__property_names__.m =================================================================== --- trunk/octave-forge/main/control/inst/@tf/__property_names__.m 2012-05-07 13:38:26 UTC (rev 10375) +++ trunk/octave-forge/main/control/inst/@tf/__property_names__.m 2012-05-08 07:56:26 UTC (rev 10376) @@ -30,12 +30,14 @@ ## cell vector of tf-specific properties props = {"num"; "den"; - "tfvar"}; + "tfvar"; + "inv"}; ## cell vector of tf-specific assignable values vals = {"p-by-m cell array of row vectors (m = number of inputs)"; "p-by-m cell array of row vectors (p = number of outputs)"; - "string (usually s or z)"}; + "string (usually s or z)"; + "logical (true for negative powers of TF variable)"}; if (nargin == 1) [ltiprops, ltivals] = __property_names__ (sys.lti); Modified: trunk/octave-forge/main/control/inst/@tf/__set__.m =================================================================== --- trunk/octave-forge/main/control/inst/@tf/__set__.m 2012-05-07 13:38:26 UTC (rev 10375) +++ trunk/octave-forge/main/control/inst/@tf/__set__.m 2012-05-08 07:56:26 UTC (rev 10376) @@ -43,10 +43,12 @@ endif case "inv" - if (isscalar (val)) + if (! isdt (sys)) + error ("tf: set: property 'inv' requires discrete-time system"); + elseif (! isscalar (val)) + error ("tf: set: property 'inv' must be a scalar logical"); + else sys.inv = logical (val); - else - error ("tf: set: property 'inv' must be a scalar logical"); endif otherwise This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
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-07 13:22:07
|
Revision: 10374 http://octave.svn.sourceforge.net/octave/?rev=10374&view=rev Author: paramaniac Date: 2012-05-07 13:22:01 +0000 (Mon, 07 May 2012) Log Message: ----------- control-devel: support miso arx identification Modified Paths: -------------- trunk/octave-forge/extra/control-devel/inst/arx.m Modified: trunk/octave-forge/extra/control-devel/inst/arx.m =================================================================== --- trunk/octave-forge/extra/control-devel/inst/arx.m 2012-05-07 12:39:54 UTC (rev 10373) +++ trunk/octave-forge/extra/control-devel/inst/arx.m 2012-05-07 13:22:01 UTC (rev 10374) @@ -15,20 +15,24 @@ if (! is_real_scalar (na, nb)) error ("arx: "); ## Test for integers + ## numel (nb) == size (dat, 3) endif ## TODO: handle MIMO and multi-experiment data + [~, p, m, ex] = size (dat); Y = dat.y{1}; U = dat.u{1}; Ts = dat.tsam{1}; - n = max (na, nb); + max_nb = max (nb); + n = max (na, max_nb); ## avoid warning: toeplitz: column wins anti-diagonal conflict PhiY = toeplitz (Y(1:end-1, :), [Y(1, :); zeros(na-1, 1)]); - PhiU = toeplitz (U(1:end-1, :), [U(1, :); zeros(nb-1, 1)]); - Phi = [-PhiY, PhiU]; + ## 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, :) ## Theta = Phi \ Y(n+1:end, :); # naive formula @@ -48,8 +52,15 @@ Theta = V * (S .\ (U' * Y(n+1:end, :))); # U' is the conjugate transpose A = [1; Theta(1:na)]; # a0 = 1, a1 = Theta(1), an = Theta(n) - B = [0; Theta(na+1:end)]; # b0 = 0 (leading zero required by filt) + ## B = [0; Theta(na+1:end)]; # b0 = 0 (leading zero required by filt) + ThetaB = Theta(na+1:end); + B = mat2cell (ThetaB, nb); + B = reshape (B, 1, []); + B = cellfun (@(x) [0; x], B, "uniformoutput", false); - sys = filt ({B, 1}, {A, A}, Ts); + ## sys = filt ({B, 1}, {A, A}, Ts); + num = [B, {1}]; + den = repmat ({A}, 1, m+1); + sys = filt (num, den, Ts); 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. |