From: <par...@us...> - 2012-05-15 13:58:23
|
Revision: 10442 http://octave.svn.sourceforge.net/octave/?rev=10442&view=rev Author: paramaniac Date: 2012-05-15 13:58:12 +0000 (Tue, 15 May 2012) Log Message: ----------- control-devel: use svd for multi-experiment data Modified Paths: -------------- trunk/octave-forge/extra/control-devel/devel/CDplayerARX.m trunk/octave-forge/extra/control-devel/inst/arx.m Modified: trunk/octave-forge/extra/control-devel/devel/CDplayerARX.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/CDplayerARX.m 2012-05-15 12:51:58 UTC (rev 10441) +++ trunk/octave-forge/extra/control-devel/devel/CDplayerARX.m 2012-05-15 13:58:12 UTC (rev 10442) @@ -51,7 +51,7 @@ dat = iddata (Y, U) % [sys, x0] = ident (dat, 15, 8) % s=15, n=8 -sys = arx (dat, 4, [4 4]) +sys = arx (dat, 4, 4) %[y, t] = lsim (sys, U, [], x0); %[y, t] = lsim (sys(:,1:2), U); Modified: trunk/octave-forge/extra/control-devel/inst/arx.m =================================================================== --- trunk/octave-forge/extra/control-devel/inst/arx.m 2012-05-15 12:51:58 UTC (rev 10441) +++ trunk/octave-forge/extra/control-devel/inst/arx.m 2012-05-15 13:58:12 UTC (rev 10442) @@ -102,25 +102,7 @@ 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 - S = diag (S); # extract main diagonal - r = sum (S > eps*S(1)); - if (r < length (S)) - warning ("arx: rank-deficient coefficient matrix"); - warning ("sampling time too small"); - warning ("persistence of excitation"); - endif - V = V(:, 1:r); - S = S(1:r); - U = U(:, 1:r); - theta = V * (S .\ (U' * y{1}(n(i)+1:end, i))); # U' is the conjugate transpose + theta = __ls_svd__ (phi{1}, y{1}(n(i)+1:end, i)); else ## multi-experiment dataset ## TODO: find more sophisticated formula than @@ -136,7 +118,31 @@ PhiTY = plus (tmp{:}); ## pseudoinverse Theta = C \ Phi'Y - theta = C \ PhiTY; + theta = __ls_svd__ (C, PhiTY); endif endfunction + + +function x = __ls_svd__ (A, b) + + ## 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 (A, 0); # 0 for "economy size" decomposition + S = diag (S); # extract main diagonal + r = sum (S > eps*S(1)); + if (r < length (S)) + warning ("arx: rank-deficient coefficient matrix"); + warning ("sampling time too small"); + warning ("persistence of excitation"); + endif + V = V(:, 1:r); + S = S(1:r); + U = U(:, 1:r); + x = V * (S .\ (U' * b)); # U' is the conjugate transpose + +endfunction This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |