## octave-cvsupdate

 [Octave-cvsupdate] SF.net SVN: octave:[5987] trunk/octave-forge/main/control/inst/svplot.m From: - 2009-07-06 08:58:16 ```Revision: 5987 http://octave.svn.sourceforge.net/octave/?rev=5987&view=rev Author: slackydeb Date: 2009-07-06 08:58:15 +0000 (Mon, 06 Jul 2009) Log Message: ----------- Little coding style fixes. Modified Paths: -------------- trunk/octave-forge/main/control/inst/svplot.m Modified: trunk/octave-forge/main/control/inst/svplot.m =================================================================== --- trunk/octave-forge/main/control/inst/svplot.m 2009-07-06 08:26:50 UTC (rev 5986) +++ trunk/octave-forge/main/control/inst/svplot.m 2009-07-06 08:58:15 UTC (rev 5987) @@ -51,12 +51,12 @@ function [sigma_min_r, sigma_max_r, w_r] = svplot (sys, w) if (nargin < 1 || nargin > 2) - print_usage(); + print_usage (); endif ## Get State Space Matrices - sys = sysupdate(sys, "ss"); - [A, B, C, D] = sys2ss(sys); + sys = sysupdate (sys, "ss"); + [A, B, C, D] = sys2ss (sys); ## Find interesting Frequency Range w if not specified @@ -64,12 +64,12 @@ j = 1; # Index Counter - for m = 1 : size(B, 2) # For all Inputs - for p = 1 : size(C, 1) # For all Outputs + for m = 1 : size (B, 2) # For all Inputs + for p = 1 : size (C, 1) # For all Outputs ## sys2zp and bode_bounds don't work for MIMO Systems! - sysp = sysprune(sys, p, m); - [zer, pol, k, tsam, inname, outname] = sys2zp(sysp); + sysp = sysprune (sys, p, m); + [zer, pol, k, tsam, inname, outname] = sys2zp (sysp); if (tsam == 0) DIGITAL = 0; @@ -77,7 +77,7 @@ DIGITAL = 1; endif - [wmin, wmax] = bode_bounds(zer, pol, DIGITAL, tsam); + [wmin, wmax] = bode_bounds (zer, pol, DIGITAL, tsam); w_min(j) = wmin; w_max(j) = wmax; @@ -85,37 +85,37 @@ endfor endfor - dec_min = min(w_min); # Begin Plot at 10^dec_min rad/s - dec_max = max(w_max); # End Plot at 10^dec_min rad/s + dec_min = min (w_min); # Begin Plot at 10^dec_min rad/s + dec_max = max (w_max); # End Plot at 10^dec_min rad/s n_steps = 1000; # Number of Steps - - w = logspace(dec_min, dec_max, n_steps); # rad/s + + w = logspace (dec_min, dec_max, n_steps); # rad/s endif ## Repeat for every Frequency w - for k = 1 : size(w, 2) + for k = 1 : size (w, 2) ## Frequency Response Matrix - P = C * inv(i * w(k) * eye(size(A)) - A) * B + D; + P = C * inv (i * w(k) * eye (size(A)) - A) * B + D; ## Singular Value Decomposition - sigma = svd(P); - sigma_min(k) = min(sigma); - sigma_max(k) = max(sigma); + sigma = svd (P); + sigma_min(k) = min (sigma); + sigma_max(k) = max (sigma); endfor if (nargout < 1) # Plot the Information ## Convert to dB for Plotting - sigma_min_db = 20 * log10(sigma_min); - sigma_max_db = 20 * log10(sigma_max); + sigma_min_db = 20 * log10 (sigma_min); + sigma_max_db = 20 * log10 (sigma_max); ## Plot Results - semilogx(w, sigma_min_db, 'b', w, sigma_max_db, 'b') - title('Singular Values') - xlabel('Frequency [rad/s]') - ylabel('Singular Values [dB]') + semilogx (w, sigma_min_db, 'b', w, sigma_max_db, 'b') + title ('Singular Values') + xlabel ('Frequency [rad/s]') + ylabel ('Singular Values [dB]') grid on else # Return Values sigma_min_r = sigma_min; This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. ```
 [Octave-cvsupdate] SF.net SVN: octave:[5990] trunk/octave-forge/main/control/inst/svplot.m From: - 2009-07-08 07:09:36 ```Revision: 5990 http://octave.svn.sourceforge.net/octave/?rev=5990&view=rev Author: paramaniac Date: 2009-07-08 07:09:29 +0000 (Wed, 08 Jul 2009) Log Message: ----------- cosmetic changes in order to improve readability of the code Modified Paths: -------------- trunk/octave-forge/main/control/inst/svplot.m Modified: trunk/octave-forge/main/control/inst/svplot.m =================================================================== --- trunk/octave-forge/main/control/inst/svplot.m 2009-07-07 18:23:28 UTC (rev 5989) +++ trunk/octave-forge/main/control/inst/svplot.m 2009-07-08 07:09:29 UTC (rev 5990) @@ -71,8 +71,8 @@ I = eye (size (A)); - ## Find interesting Frequency Range w if not specified - if (nargin < 2) + ## Find interesting frequency range w if not specified + if (nargin == 1) ## Since no frequency vector w has been specified, the interesting ## decades of the sigma plot have to be found. The already existing @@ -83,38 +83,38 @@ ## Afterwards, the min/max decade can be found by the min()/max() ## command. - j = 1; # Index Counter + j = 1; # Index counter (number of SISO systems) - for m = 1 : size (B, 2) # For all Inputs - for p = 1 : size (C, 1) # For all Outputs + for m = 1 : size (B, 2) # For all inputs + for p = 1 : size (C, 1) # For all outputs - sysp = sysprune (sys, p, m); + sysp = sysprune (sys, p, m); # Create SISO system [zer, pol, k, tsam, inname, outname] = sys2zp (sysp); - if (tsam == 0) + if (tsam == 0) # Continuous system DIGITAL = 0; - else + else # Discrete system DIGITAL = 1; ## FIXME: Discrete systems not yet tested! endif [wmin, wmax] = bode_bounds (zer, pol, DIGITAL, tsam); - w_min(j) = wmin; - w_max(j) = wmax; + w_min(j) = wmin; # Save result for later evaluation + w_max(j) = wmax; # dito j++; endfor endfor - dec_min = min (w_min); # Begin Plot at 10^dec_min rad/s - dec_max = max (w_max); # End Plot at 10^dec_min rad/s - n_steps = 1000; # Number of Frequencies evaluated for Plotting + dec_min = min (w_min); # Begin plot at 10^dec_min [rad/s] + dec_max = max (w_max); # End plot at 10^dec_min [rad/s] + n_steps = 1000; # Number of frequencies evaluated for plotting - w = logspace (dec_min, dec_max, n_steps); # rad/s + w = logspace (dec_min, dec_max, n_steps); # [rad/s] endif - ## Repeat for every Frequency w + ## Repeat for every frequency w for k = 1 : size (w, 2) ## Frequency Response Matrix @@ -126,19 +126,19 @@ sigma_max(k) = max (sigma); endfor - if (nargout == 0) # Plot the Information + if (nargout == 0) # Plot the information - ## Convert to dB for Plotting + ## Convert to dB for plotting sigma_min_db = 20 * log10 (sigma_min); sigma_max_db = 20 * log10 (sigma_max); - ## Plot Results + ## Plot results semilogx (w, sigma_min_db, 'b', w, sigma_max_db, 'b') title ('Singular Values') xlabel ('Frequency [rad/s]') ylabel ('Singular Values [dB]') grid on - else # Return Values + else # Return values sigma_min_r = sigma_min; sigma_max_r = sigma_max; w_r = w; @@ -159,4 +159,4 @@ %! [sigma_min_obs, sigma_max_obs, w_obs] = svplot (ss (A, B, C, D), w); %!assert (sigma_min_obs, sigma_min_exp, 4*eps); # tolerance manually tweaked %!assert (sigma_max_obs, sigma_max_exp, 12*eps); # tolerance manually tweaked -%!assert (w_obs, w_exp, 2*eps); +%!assert (w_obs, w_exp, 2*eps); \ No newline at end of file This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. ```
 [Octave-cvsupdate] SF.net SVN: octave:[5991] trunk/octave-forge/main/control/inst/svplot.m From: - 2009-07-08 08:09:50 ```Revision: 5991 http://octave.svn.sourceforge.net/octave/?rev=5991&view=rev Author: paramaniac Date: 2009-07-08 08:09:47 +0000 (Wed, 08 Jul 2009) Log Message: ----------- removed unclear academic expression Modified Paths: -------------- trunk/octave-forge/main/control/inst/svplot.m Modified: trunk/octave-forge/main/control/inst/svplot.m =================================================================== --- trunk/octave-forge/main/control/inst/svplot.m 2009-07-08 07:09:29 UTC (rev 5990) +++ trunk/octave-forge/main/control/inst/svplot.m 2009-07-08 08:09:47 UTC (rev 5991) @@ -100,8 +100,9 @@ [wmin, wmax] = bode_bounds (zer, pol, DIGITAL, tsam); - w_min(j) = wmin; # Save result for later evaluation - w_max(j) = wmax; # dito + # Save result for later evaluation + w_min(j) = wmin; + w_max(j) = wmax; j++; endfor endfor This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. ```
 [Octave-cvsupdate] SF.net SVN: octave:[5994] trunk/octave-forge/main/control/inst/svplot.m From: - 2009-07-08 19:31:28 ```Revision: 5994 http://octave.svn.sourceforge.net/octave/?rev=5994&view=rev Author: paramaniac Date: 2009-07-08 19:31:26 +0000 (Wed, 08 Jul 2009) Log Message: ----------- added warning for discrete systems, used internal function instead of own code Modified Paths: -------------- trunk/octave-forge/main/control/inst/svplot.m Modified: trunk/octave-forge/main/control/inst/svplot.m =================================================================== --- trunk/octave-forge/main/control/inst/svplot.m 2009-07-08 19:30:15 UTC (rev 5993) +++ trunk/octave-forge/main/control/inst/svplot.m 2009-07-08 19:31:26 UTC (rev 5994) @@ -69,8 +69,15 @@ sys = sysupdate (sys, "ss"); [A, B, C, D] = sys2ss (sys); I = eye (size (A)); + + ## Get system information + [n_c_states, n_d_states, n_in, n_out] = sysdimensions (sys); + + ## Warning for discrete systems, see FIXME below! + if (n_d_states != 0) + warning ("discrete systems have not been tested!"); + endif - ## Find interesting frequency range w if not specified if (nargin == 1) @@ -85,8 +92,8 @@ j = 1; # Index counter (number of SISO systems) - for m = 1 : size (B, 2) # For all inputs - for p = 1 : size (C, 1) # For all outputs + for m = 1 : n_in # For all inputs + for p = 1 : n_out # For all outputs sysp = sysprune (sys, p, m); # Create SISO system [zer, pol, k, tsam, inname, outname] = sys2zp (sysp); @@ -100,7 +107,7 @@ [wmin, wmax] = bode_bounds (zer, pol, DIGITAL, tsam); - # Save result for later evaluation + ## Save result for later evaluation w_min(j) = wmin; w_max(j) = wmax; j++; This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. ```
 [Octave-cvsupdate] SF.net SVN: octave:[5995] trunk/octave-forge/main/control/inst/svplot.m From: - 2009-07-10 05:27:44 ```Revision: 5995 http://octave.svn.sourceforge.net/octave/?rev=5995&view=rev Author: paramaniac Date: 2009-07-10 05:27:40 +0000 (Fri, 10 Jul 2009) Log Message: ----------- svplot: partial rewrite for simplification and discrete systems, bumped version Modified Paths: -------------- trunk/octave-forge/main/control/inst/svplot.m Modified: trunk/octave-forge/main/control/inst/svplot.m =================================================================== --- trunk/octave-forge/main/control/inst/svplot.m 2009-07-08 19:31:26 UTC (rev 5994) +++ trunk/octave-forge/main/control/inst/svplot.m 2009-07-10 05:27:40 UTC (rev 5995) @@ -45,7 +45,7 @@ ## @end deftypefn ## Author: Lukas Reichlin -## Version: 0.1 +## Version: 0.2 function [sigma_min_r, sigma_max_r, w_r] = svplot (sys, w) @@ -65,15 +65,32 @@ endif endif + if (is_siso (sys)) + warning ("sys is a SISO system. You might want to use bode(sys) instead."); + endif + ## Get State Space Matrices sys = sysupdate (sys, "ss"); [A, B, C, D] = sys2ss (sys); I = eye (size (A)); + %## Reduce to upper hessenberg form for efficiency + %[p, A] = hess (A); + %B = p' * B; + %C = C * p; + ## Get system information [n_c_states, n_d_states, n_in, n_out] = sysdimensions (sys); + digital = is_digital(sys, 2); + tsam = sysgettsam(sys); - ## Warning for discrete systems, see FIXME below! + ## Error for mixed systems + if (digital == -1) + error ("system must be either purely continuous or purely discrete"); + endif + + ## FIXME: Discrete systems not tested! + ## Warning for discrete systems if (n_d_states != 0) warning ("discrete systems have not been tested!"); endif @@ -83,57 +100,37 @@ ## Since no frequency vector w has been specified, the interesting ## decades of the sigma plot have to be found. The already existing - ## function bode_bounds does exactly that, unfortunately for SISO - ## systems only. Therefore, a SISO system from every input m to - ## every output p is created. Then for every SISO system the - ## interesting frequency range is calculated by bode_bounds. - ## Afterwards, the min/max decade can be found by the min()/max() - ## command. + ## function bode_bounds does exactly that. - j = 1; # Index counter (number of SISO systems) + zer = tzero (sys); + pol = eig (A); - for m = 1 : n_in # For all inputs - for p = 1 : n_out # For all outputs - - sysp = sysprune (sys, p, m); # Create SISO system - [zer, pol, k, tsam, inname, outname] = sys2zp (sysp); - - if (tsam == 0) # Continuous system - DIGITAL = 0; - else # Discrete system - DIGITAL = 1; - ## FIXME: Discrete systems not yet tested! - endif - - [wmin, wmax] = bode_bounds (zer, pol, DIGITAL, tsam); - - ## Save result for later evaluation - w_min(j) = wmin; - w_max(j) = wmax; - j++; - endfor - endfor - - dec_min = min (w_min); # Begin plot at 10^dec_min [rad/s] - dec_max = max (w_max); # End plot at 10^dec_min [rad/s] + ## Begin plot at 10^dec_min, end plot at 10^dec_min [rad/s] + [dec_min, dec_max] = bode_bounds (zer, pol, digital, tsam); + n_steps = 1000; # Number of frequencies evaluated for plotting w = logspace (dec_min, dec_max, n_steps); # [rad/s] endif + + if (! digital) # continuous system + jw = i * w; + else # discrete system + jw = exp (i*w*tsam); + endif + + ## Repeat for every frequency jw + for k = 1 : size (jw, 2) - - ## Repeat for every frequency w - for k = 1 : size (w, 2) - ## Frequency Response Matrix - P = C * inv (i*w(k)*I - A) * B + D; + P = C * inv (jw(k)*I - A) * B + D; ## Singular Value Decomposition sigma = svd (P); sigma_min(k) = min (sigma); sigma_max(k) = max (sigma); endfor - + if (nargout == 0) # Plot the information ## Convert to dB for plotting This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. ```
 [Octave-cvsupdate] SF.net SVN: octave:[5999] trunk/octave-forge/main/control/inst/svplot.m From: - 2009-07-10 15:27:51 ```Revision: 5999 http://octave.svn.sourceforge.net/octave/?rev=5999&view=rev Author: paramaniac Date: 2009-07-10 15:27:43 +0000 (Fri, 10 Jul 2009) Log Message: ----------- svplot: works now for discrete systems Modified Paths: -------------- trunk/octave-forge/main/control/inst/svplot.m Modified: trunk/octave-forge/main/control/inst/svplot.m =================================================================== --- trunk/octave-forge/main/control/inst/svplot.m 2009-07-10 07:24:24 UTC (rev 5998) +++ trunk/octave-forge/main/control/inst/svplot.m 2009-07-10 15:27:43 UTC (rev 5999) @@ -74,11 +74,6 @@ [A, B, C, D] = sys2ss (sys); I = eye (size (A)); - %## Reduce to upper hessenberg form for efficiency - %[p, A] = hess (A); - %B = p' * B; - %C = C * p; - ## Get system information [n_c_states, n_d_states, n_in, n_out] = sysdimensions (sys); digital = is_digital(sys, 2); @@ -88,12 +83,6 @@ if (digital == -1) error ("system must be either purely continuous or purely discrete"); endif - - ## FIXME: Discrete systems not tested! - ## Warning for discrete systems - if (n_d_states != 0) - warning ("discrete systems have not been tested!"); - endif ## Find interesting frequency range w if not specified if (nargin == 1) This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. ```
 [Octave-cvsupdate] SF.net SVN: octave:[6000] trunk/octave-forge/main/control/inst/svplot.m From: - 2009-07-11 07:30:33 ```Revision: 6000 http://octave.svn.sourceforge.net/octave/?rev=6000&view=rev Author: paramaniac Date: 2009-07-11 07:30:31 +0000 (Sat, 11 Jul 2009) Log Message: ----------- svplot: minor changes in code and documentation Modified Paths: -------------- trunk/octave-forge/main/control/inst/svplot.m Modified: trunk/octave-forge/main/control/inst/svplot.m =================================================================== --- trunk/octave-forge/main/control/inst/svplot.m 2009-07-10 15:27:43 UTC (rev 5999) +++ trunk/octave-forge/main/control/inst/svplot.m 2009-07-11 07:30:31 UTC (rev 6000) @@ -17,7 +17,7 @@ ## @deftypefn{Function File} {[@var{sigma_min}, @var{sigma_max}, @var{w}] =} svplot (@var{sys}) ## @deftypefnx{Function File} {[@var{sigma_min}, @var{sigma_max}, @var{w}] =} svplot (@var{sys}, @var{w}) ## If no output arguments are given, the singular value plot of a MIMO -## system over a range of frequencies is printed on the screen; +## system is printed on the screen; ## otherwise, the singular values of the system data structure are ## computed and returned. ## @@ -41,7 +41,7 @@ ## Vector of frequency values used. ## @end table ## -## @seealso{is_digital} +## @seealso{bode, svd, bode_bounds, is_digital} ## @end deftypefn ## Author: Lukas Reichlin @@ -102,17 +102,17 @@ w = logspace (dec_min, dec_max, n_steps); # [rad/s] endif - if (! digital) # continuous system - jw = i * w; - else # discrete system - jw = exp (i*w*tsam); + if (digital) # discrete system + s = exp (i * w * tsam); + else # continuous system + s = i * w; endif - ## Repeat for every frequency jw - for k = 1 : size (jw, 2) + ## Repeat for every frequency s + for k = 1 : size (s, 2) ## Frequency Response Matrix - P = C * inv (jw(k)*I - A) * B + D; + P = C * inv (s(k)*I - A) * B + D; ## Singular Value Decomposition sigma = svd (P); This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. ```
 [Octave-cvsupdate] SF.net SVN: octave:[6002] trunk/octave-forge/main/control/inst/svplot.m From: - 2009-07-12 02:52:35 ```Revision: 6002 http://octave.svn.sourceforge.net/octave/?rev=6002&view=rev Author: paramaniac Date: 2009-07-12 02:52:26 +0000 (Sun, 12 Jul 2009) Log Message: ----------- svplot: polishing code Modified Paths: -------------- trunk/octave-forge/main/control/inst/svplot.m Modified: trunk/octave-forge/main/control/inst/svplot.m =================================================================== --- trunk/octave-forge/main/control/inst/svplot.m 2009-07-11 19:20:27 UTC (rev 6001) +++ trunk/octave-forge/main/control/inst/svplot.m 2009-07-12 02:52:26 UTC (rev 6002) @@ -24,8 +24,8 @@ ## @strong{Inputs} ## @table @var ## @item sys -## System data structure (must be either purely continuous or discrete; -## see @code{is_digital}). +## System data structure. Must be either purely continuous or discrete; +## see @code{is_digital}. ## @item w ## Optional vector of frequency values. If @var{w} is not specified, it ## will be calculated by @code{bode_bounds}. @@ -56,32 +56,31 @@ endif if(! isstruct (sys)) - error ("first argument sys must be a system data structure"); + error ("svplot: first argument sys must be a system data structure"); endif if (nargin == 2) if (! isvector (w)) - error ("second argument w must be a vector of frequencies"); + error ("svplot: second argument w must be a vector of frequencies"); endif endif if (is_siso (sys)) - warning ("sys is a SISO system. You might want to use bode(sys) instead."); + warning ("svplot: sys is a SISO system. You might want to use bode(sys) instead."); endif - ## Get State Space Matrices + ## Get state space matrices sys = sysupdate (sys, "ss"); [A, B, C, D] = sys2ss (sys); I = eye (size (A)); ## Get system information - [n_c_states, n_d_states, n_in, n_out] = sysdimensions (sys); digital = is_digital(sys, 2); - tsam = sysgettsam(sys); + t_sam = sysgettsam(sys); ## Error for mixed systems if (digital == -1) - error ("system must be either purely continuous or purely discrete"); + error ("svplot: system must be either purely continuous or purely discrete"); endif ## Find interesting frequency range w if not specified @@ -94,17 +93,17 @@ zer = tzero (sys); pol = eig (A); - ## Begin plot at 10^dec_min, end plot at 10^dec_min [rad/s] - [dec_min, dec_max] = bode_bounds (zer, pol, digital, tsam); + ## Begin plot at 10^dec_min, end plot at 10^dec_max [rad/s] + [dec_min, dec_max] = bode_bounds (zer, pol, digital, t_sam); - n_steps = 1000; # Number of frequencies evaluated for plotting + n_freq = 1000; # Number of frequencies evaluated for plotting - w = logspace (dec_min, dec_max, n_steps); # [rad/s] + w = logspace (dec_min, dec_max, n_freq); # [rad/s] endif - if (digital) # discrete system - s = exp (i * w * tsam); - else # continuous system + if (digital) # Discrete system + s = exp (i * w * t_sam); + else # Continuous system s = i * w; endif @@ -112,10 +111,10 @@ for k = 1 : size (s, 2) ## Frequency Response Matrix - P = C * inv (s(k)*I - A) * B + D; + G = C * inv (s(k)*I - A) * B + D; ## Singular Value Decomposition - sigma = svd (P); + sigma = svd (G); sigma_min(k) = min (sigma); sigma_max(k) = max (sigma); endfor @@ -125,11 +124,18 @@ ## Convert to dB for plotting sigma_min_db = 20 * log10 (sigma_min); sigma_max_db = 20 * log10 (sigma_max); + + ## Determine xlabel + if (digital) + xl_str = sprintf ('Frequency [rad/s] Pi / T = %g', pi/t_sam); + else + xl_str = 'Frequency [rad/s]'; + endif ## Plot results semilogx (w, sigma_min_db, 'b', w, sigma_max_db, 'b') title ('Singular Values') - xlabel ('Frequency [rad/s]') + xlabel (xl_str) ylabel ('Singular Values [dB]') grid on else # Return values This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. ```
 [Octave-cvsupdate] SF.net SVN: octave:[6071] trunk/octave-forge/main/control/inst/svplot.m From: - 2009-08-02 20:55:43 ```Revision: 6071 http://octave.svn.sourceforge.net/octave/?rev=6071&view=rev Author: paramaniac Date: 2009-08-02 20:55:25 +0000 (Sun, 02 Aug 2009) Log Message: ----------- svplot: extended features Modified Paths: -------------- trunk/octave-forge/main/control/inst/svplot.m Modified: trunk/octave-forge/main/control/inst/svplot.m =================================================================== --- trunk/octave-forge/main/control/inst/svplot.m 2009-08-02 17:23:45 UTC (rev 6070) +++ trunk/octave-forge/main/control/inst/svplot.m 2009-08-02 20:55:25 UTC (rev 6071) @@ -16,6 +16,7 @@ ## -*- texinfo -*- ## @deftypefn{Function File} {[@var{sigma_min}, @var{sigma_max}, @var{w}] =} svplot (@var{sys}) ## @deftypefnx{Function File} {[@var{sigma_min}, @var{sigma_max}, @var{w}] =} svplot (@var{sys}, @var{w}) +## @deftypefnx{Function File} {[@var{sigma_min}, @var{sigma_max}, @var{w}] =} svplot (@var{sys}, @var{[]}, @var{ptype}) ## If no output arguments are given, the singular value plot of a MIMO ## system is printed on the screen; ## otherwise, the singular values of the system data structure are @@ -29,6 +30,14 @@ ## @item w ## Optional vector of frequency values. If @var{w} is not specified, it ## will be calculated by @code{bode_bounds}. +## @item ptype = 0 +## Singular values of the frequency response H of system sys. Default Value. +## @item ptype = 1 +## Singular values of the frequency response inv(H); i.e. inversed system. +## @item ptype = 2 +## Singular values of the frequency response I + H; i.e. return difference. +## @item ptype = 3 +## Singular values of the frequency response I + inv(H); i.e inversed complementary sensitivity. ## @end table ## ## @strong{Outputs} @@ -45,13 +54,13 @@ ## @end deftypefn ## Author: Lukas Reichlin -## Version: 0.2 +## Version: 0.3 -function [sigma_min_r, sigma_max_r, w_r] = svplot (sys, w) +function [sigma_min_r, sigma_max_r, w_r] = svplot (sys, w, ptype) ## Check whether arguments are OK - if (nargin < 1 || nargin > 2) + if (nargin < 1 || nargin > 3) print_usage (); endif @@ -59,10 +68,10 @@ error ("svplot: first argument sys must be a system data structure"); endif - if (nargin == 2) - if (! isvector (w)) - error ("svplot: second argument w must be a vector of frequencies"); - endif + if (nargin == 1) + w = []; + elseif (! isvector (w) && ! isempty (w)) + error ("svplot: second argument w must be a vector of frequencies"); endif if (is_siso (sys)) @@ -82,9 +91,28 @@ if (digital == -1) error ("svplot: system must be either purely continuous or purely discrete"); endif + + ## Handle plot type + if (nargin == 3) + if (isfloat (ptype)) # Numeric constants like 2 are NOT integers in Octave! + if (ptype != 0 && ptype != 1 && ptype != 2 && ptype != 3) + error ("svplot: third argument ptype must be 0, 1, 2 or 3"); + endif + else + error ("svplot: third argument ptype must be a number"); + endif + + [n_c, n_d, m, p] = sysdimensions (sys); + if (m != p) + error ("svplot: system must be square for ptype 1, 2 or 3"); + endif + J = eye(m); + else + ptype = 0; # Default value + endif ## Find interesting frequency range w if not specified - if (nargin == 1) + if (isempty (w)) ## Since no frequency vector w has been specified, the interesting ## decades of the sigma plot have to be found. The already existing @@ -96,30 +124,52 @@ ## Begin plot at 10^dec_min, end plot at 10^dec_max [rad/s] [dec_min, dec_max] = bode_bounds (zer, pol, digital, t_sam); - n_freq = 1000; # Number of frequencies evaluated for plotting + n_freq = 1000; # Number of frequencies evaluated for plotting - w = logspace (dec_min, dec_max, n_freq); # [rad/s] + w = logspace (dec_min, dec_max, n_freq); # [rad/s] endif - if (digital) # Discrete system + if (digital) # Discrete system s = exp (i * w * t_sam); - else # Continuous system + else # Continuous system s = i * w; endif - - ## Repeat for every frequency s - for k = 1 : size (s, 2) - - ## Frequency Response Matrix - G = C * inv (s(k)*I - A) * B + D; - - ## Singular Value Decomposition - sigma = svd (G); - sigma_min(k) = min (sigma); - sigma_max(k) = max (sigma); - endfor + + switch (ptype) + case 0 + for k = 1 : size (s, 2) # Repeat for every frequency s + H = C * inv (s(k)*I - A) * B + D; # Frequency Response Matrix + sigma = svd (H); # Singular Value Decomposition + sigma_min(k) = min (sigma); + sigma_max(k) = max (sigma); + endfor - if (nargout == 0) # Plot the information + case 1 + for k = 1 : size (s, 2) + H = inv (C * inv (s(k)*I - A) * B + D); + sigma = svd (H); + sigma_min(k) = min (sigma); + sigma_max(k) = max (sigma); + endfor + + case 2 + for k = 1 : size (s, 2) + H = J + C * inv (s(k)*I - A) * B + D; + sigma = svd (H); + sigma_min(k) = min (sigma); + sigma_max(k) = max (sigma); + endfor + + case 3 + for k = 1 : size (s, 2) + H = J + inv (C * inv (s(k)*I - A) * B + D); + sigma = svd (H); + sigma_min(k) = min (sigma); + sigma_max(k) = max (sigma); + endfor + endswitch + + if (nargout == 0) # Plot the information ## Convert to dB for plotting sigma_min_db = 20 * log10 (sigma_min); @@ -138,7 +188,7 @@ xlabel (xl_str) ylabel ('Singular Values [dB]') grid on - else # Return values + else # Return values sigma_min_r = sigma_min; sigma_max_r = sigma_max; w_r = w; This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. ```
 [Octave-cvsupdate] SF.net SVN: octave:[6072] trunk/octave-forge/main/control/inst/svplot.m From: - 2009-08-02 21:17:12 ```Revision: 6072 http://octave.svn.sourceforge.net/octave/?rev=6072&view=rev Author: paramaniac Date: 2009-08-02 21:17:03 +0000 (Sun, 02 Aug 2009) Log Message: ----------- svplot: added forgotten comments Modified Paths: -------------- trunk/octave-forge/main/control/inst/svplot.m Modified: trunk/octave-forge/main/control/inst/svplot.m =================================================================== --- trunk/octave-forge/main/control/inst/svplot.m 2009-08-02 20:55:25 UTC (rev 6071) +++ trunk/octave-forge/main/control/inst/svplot.m 2009-08-02 21:17:03 UTC (rev 6072) @@ -37,7 +37,7 @@ ## @item ptype = 2 ## Singular values of the frequency response I + H; i.e. return difference. ## @item ptype = 3 -## Singular values of the frequency response I + inv(H); i.e inversed complementary sensitivity. +## Singular values of the frequency response I + inv(H); i.e. inversed complementary sensitivity. ## @end table ## ## @strong{Outputs} @@ -136,7 +136,7 @@ endif switch (ptype) - case 0 + case 0 # Default system for k = 1 : size (s, 2) # Repeat for every frequency s H = C * inv (s(k)*I - A) * B + D; # Frequency Response Matrix sigma = svd (H); # Singular Value Decomposition @@ -144,7 +144,7 @@ sigma_max(k) = max (sigma); endfor - case 1 + case 1 # Inversed system for k = 1 : size (s, 2) H = inv (C * inv (s(k)*I - A) * B + D); sigma = svd (H); @@ -152,7 +152,7 @@ sigma_max(k) = max (sigma); endfor - case 2 + case 2 # Return difference for k = 1 : size (s, 2) H = J + C * inv (s(k)*I - A) * B + D; sigma = svd (H); @@ -160,7 +160,7 @@ sigma_max(k) = max (sigma); endfor - case 3 + case 3 # Inversed complementary sensitivity for k = 1 : size (s, 2) H = J + inv (C * inv (s(k)*I - A) * B + D); sigma = svd (H); This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. ```
 [Octave-cvsupdate] SF.net SVN: octave:[6084] trunk/octave-forge/main/control/inst/svplot.m From: - 2009-08-07 21:30:02 ```Revision: 6084 http://octave.svn.sourceforge.net/octave/?rev=6084&view=rev Author: paramaniac Date: 2009-08-07 21:29:53 +0000 (Fri, 07 Aug 2009) Log Message: ----------- svplot: preallocation Modified Paths: -------------- trunk/octave-forge/main/control/inst/svplot.m Modified: trunk/octave-forge/main/control/inst/svplot.m =================================================================== --- trunk/octave-forge/main/control/inst/svplot.m 2009-08-07 11:52:46 UTC (rev 6083) +++ trunk/octave-forge/main/control/inst/svplot.m 2009-08-07 21:29:53 UTC (rev 6084) @@ -135,9 +135,13 @@ s = i * w; endif + l_s = length (s); + sigma_min = zeros (1, l_s); + sigma_max = zeros (1, l_s); + switch (ptype) case 0 # Default system - for k = 1 : size (s, 2) # Repeat for every frequency s + for k = 1 : l_s # Repeat for every frequency s H = C * inv (s(k)*I - A) * B + D; # Frequency Response Matrix sigma = svd (H); # Singular Value Decomposition sigma_min(k) = min (sigma); @@ -145,7 +149,7 @@ endfor case 1 # Inversed system - for k = 1 : size (s, 2) + for k = 1 : l_s H = inv (C * inv (s(k)*I - A) * B + D); sigma = svd (H); sigma_min(k) = min (sigma); @@ -153,7 +157,7 @@ endfor case 2 # Return difference - for k = 1 : size (s, 2) + for k = 1 : l_s H = J + C * inv (s(k)*I - A) * B + D; sigma = svd (H); sigma_min(k) = min (sigma); @@ -161,7 +165,7 @@ endfor case 3 # Inversed complementary sensitivity - for k = 1 : size (s, 2) + for k = 1 : l_s H = J + inv (C * inv (s(k)*I - A) * B + D); sigma = svd (H); sigma_min(k) = min (sigma); This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. ```
 [Octave-cvsupdate] SF.net SVN: octave:[6094] trunk/octave-forge/main/control/inst/svplot.m From: - 2009-08-13 08:54:09 ```Revision: 6094 http://octave.svn.sourceforge.net/octave/?rev=6094&view=rev Author: paramaniac Date: 2009-08-13 08:54:00 +0000 (Thu, 13 Aug 2009) Log Message: ----------- svplot: minor changes in documentation Modified Paths: -------------- trunk/octave-forge/main/control/inst/svplot.m Modified: trunk/octave-forge/main/control/inst/svplot.m =================================================================== --- trunk/octave-forge/main/control/inst/svplot.m 2009-08-13 08:42:20 UTC (rev 6093) +++ trunk/octave-forge/main/control/inst/svplot.m 2009-08-13 08:54:00 UTC (rev 6094) @@ -35,9 +35,11 @@ ## @item ptype = 1 ## Singular values of the frequency response inv(H); i.e. inversed system. ## @item ptype = 2 -## Singular values of the frequency response I + H; i.e. return difference. +## Singular values of the frequency response I + H; i.e. inversed sensitivity +## (or return difference) if H = P * C. ## @item ptype = 3 -## Singular values of the frequency response I + inv(H); i.e. inversed complementary sensitivity. +## Singular values of the frequency response I + inv(H); i.e. inversed complementary +## sensitivity if H = P * C. ## @end table ## ## @strong{Outputs} @@ -156,7 +158,7 @@ sigma_max(k) = max (sigma); endfor - case 2 # Return difference + case 2 # Inversed sensitivity for k = 1 : l_s H = J + C * inv (s(k)*I - A) * B + D; sigma = svd (H); This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. ```
 [Octave-cvsupdate] SF.net SVN: octave:[6106] trunk/octave-forge/main/control/inst/svplot.m From: - 2009-08-17 08:06:04 ```Revision: 6106 http://octave.svn.sourceforge.net/octave/?rev=6106&view=rev Author: slackydeb Date: 2009-08-17 08:05:57 +0000 (Mon, 17 Aug 2009) Log Message: ----------- svplot: clean. Modified Paths: -------------- trunk/octave-forge/main/control/inst/svplot.m Modified: trunk/octave-forge/main/control/inst/svplot.m =================================================================== --- trunk/octave-forge/main/control/inst/svplot.m 2009-08-17 07:42:46 UTC (rev 6105) +++ trunk/octave-forge/main/control/inst/svplot.m 2009-08-17 08:05:57 UTC (rev 6106) @@ -58,7 +58,6 @@ ## Author: Lukas Reichlin ## Version: 0.3 - function [sigma_min_r, sigma_max_r, w_r] = svplot (sys, w, ptype) ## Check whether arguments are OK @@ -84,16 +83,16 @@ sys = sysupdate (sys, "ss"); [A, B, C, D] = sys2ss (sys); I = eye (size (A)); - + ## Get system information digital = is_digital(sys, 2); t_sam = sysgettsam(sys); - + ## Error for mixed systems if (digital == -1) error ("svplot: system must be either purely continuous or purely discrete"); endif - + ## Handle plot type if (nargin == 3) if (isfloat (ptype)) # Numeric constants like 2 are NOT integers in Octave! @@ -103,13 +102,13 @@ else error ("svplot: third argument ptype must be a number"); endif - + [n_c, n_d, m, p] = sysdimensions (sys); if (m != p) error ("svplot: system must be square for ptype 1, 2 or 3"); endif J = eye(m); - else + else ptype = 0; # Default value endif @@ -122,25 +121,25 @@ zer = tzero (sys); pol = eig (A); - + ## Begin plot at 10^dec_min, end plot at 10^dec_max [rad/s] [dec_min, dec_max] = bode_bounds (zer, pol, digital, t_sam); - + n_freq = 1000; # Number of frequencies evaluated for plotting w = logspace (dec_min, dec_max, n_freq); # [rad/s] endif - + if (digital) # Discrete system s = exp (i * w * t_sam); else # Continuous system s = i * w; endif - + l_s = length (s); sigma_min = zeros (1, l_s); sigma_max = zeros (1, l_s); - + switch (ptype) case 0 # Default system for k = 1 : l_s # Repeat for every frequency s @@ -149,7 +148,7 @@ sigma_min(k) = min (sigma); sigma_max(k) = max (sigma); endfor - + case 1 # Inversed system for k = 1 : l_s H = inv (C * inv (s(k)*I - A) * B + D); @@ -157,7 +156,7 @@ sigma_min(k) = min (sigma); sigma_max(k) = max (sigma); endfor - + case 2 # Inversed sensitivity for k = 1 : l_s H = J + C * inv (s(k)*I - A) * B + D; @@ -165,7 +164,7 @@ sigma_min(k) = min (sigma); sigma_max(k) = max (sigma); endfor - + case 3 # Inversed complementary sensitivity for k = 1 : l_s H = J + inv (C * inv (s(k)*I - A) * B + D); @@ -174,20 +173,20 @@ sigma_max(k) = max (sigma); endfor endswitch - + if (nargout == 0) # Plot the information ## Convert to dB for plotting sigma_min_db = 20 * log10 (sigma_min); sigma_max_db = 20 * log10 (sigma_max); - + ## Determine xlabel if (digital) xl_str = sprintf ('Frequency [rad/s] Pi / T = %g', pi/t_sam); else xl_str = 'Frequency [rad/s]'; endif - + ## Plot results semilogx (w, sigma_min_db, 'b', w, sigma_max_db, 'b') title ('Singular Values') This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. ```
 [Octave-cvsupdate] SF.net SVN: octave:[6110] trunk/octave-forge/main/control/inst/svplot.m From: - 2009-08-17 08:45:14 ```Revision: 6110 http://octave.svn.sourceforge.net/octave/?rev=6110&view=rev Author: slackydeb Date: 2009-08-17 08:45:04 +0000 (Mon, 17 Aug 2009) Log Message: ----------- svplot: remove at end of line (dos2unix inst/svplot.m). Modified Paths: -------------- trunk/octave-forge/main/control/inst/svplot.m Modified: trunk/octave-forge/main/control/inst/svplot.m =================================================================== --- trunk/octave-forge/main/control/inst/svplot.m 2009-08-17 08:39:29 UTC (rev 6109) +++ trunk/octave-forge/main/control/inst/svplot.m 2009-08-17 08:45:04 UTC (rev 6110) @@ -1,217 +1,217 @@ -## Copyright (C) 2009 Lukas Reichlin -## -## 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 ;. - -## -*- texinfo -*- -## @deftypefn{Function File} {[@var{sigma_min}, @var{sigma_max}, @var{w}] =} svplot (@var{sys}) -## @deftypefnx{Function File} {[@var{sigma_min}, @var{sigma_max}, @var{w}] =} svplot (@var{sys}, @var{w}) -## @deftypefnx{Function File} {[@var{sigma_min}, @var{sigma_max}, @var{w}] =} svplot (@var{sys}, @var{[]}, @var{ptype}) -## If no output arguments are given, the singular value plot of a MIMO -## system is printed on the screen; -## otherwise, the singular values of the system data structure are -## computed and returned. -## -## @strong{Inputs} -## @table @var -## @item sys -## System data structure. Must be either purely continuous or discrete; -## see @code{is_digital}. -## @item w -## Optional vector of frequency values. If @var{w} is not specified, it -## will be calculated by @code{bode_bounds}. -## @item ptype = 0 -## Singular values of the frequency response H of system sys. Default Value. -## @item ptype = 1 -## Singular values of the frequency response inv(H); i.e. inversed system. -## @item ptype = 2 -## Singular values of the frequency response I + H; i.e. inversed sensitivity -## (or return difference) if H = P * C. -## @item ptype = 3 -## Singular values of the frequency response I + inv(H); i.e. inversed complementary -## sensitivity if H = P * C. -## @end table -## -## @strong{Outputs} -## @table @var -## @item sigma_min -## Vector of minimal singular values. -## @item sigma_max -## Vector of maximal singular values. -## @item w -## Vector of frequency values used. -## @end table -## -## @seealso{bode, svd, bode_bounds, is_digital} -## @end deftypefn - -## Author: Lukas Reichlin -## Version: 0.3 - -function [sigma_min_r, sigma_max_r, w_r] = svplot (sys, w, ptype) - - ## Check whether arguments are OK - if (nargin < 1 || nargin > 3) - print_usage (); - endif - - if(! isstruct (sys)) - error ("svplot: first argument sys must be a system data structure"); - endif - - if (nargin == 1) - w = []; - elseif (! isvector (w) && ! isempty (w)) - error ("svplot: second argument w must be a vector of frequencies"); - endif - - if (is_siso (sys)) - warning ("svplot: sys is a SISO system. You might want to use bode(sys) instead."); - endif - - ## Get state space matrices - sys = sysupdate (sys, "ss"); - [A, B, C, D] = sys2ss (sys); - I = eye (size (A)); - - ## Get system information - digital = is_digital(sys, 2); - t_sam = sysgettsam(sys); - - ## Error for mixed systems - if (digital == -1) - error ("svplot: system must be either purely continuous or purely discrete"); - endif - - ## Handle plot type - if (nargin == 3) - if (isfloat (ptype)) # Numeric constants like 2 are NOT integers in Octave! - if (ptype != 0 && ptype != 1 && ptype != 2 && ptype != 3) - error ("svplot: third argument ptype must be 0, 1, 2 or 3"); - endif - else - error ("svplot: third argument ptype must be a number"); - endif - - [n_c, n_d, m, p] = sysdimensions (sys); - if (m != p) - error ("svplot: system must be square for ptype 1, 2 or 3"); - endif - J = eye(m); - else - ptype = 0; # Default value - endif - - ## Find interesting frequency range w if not specified - if (isempty (w)) - - ## Since no frequency vector w has been specified, the interesting - ## decades of the sigma plot have to be found. The already existing - ## function bode_bounds does exactly that. - - zer = tzero (sys); - pol = eig (A); - - ## Begin plot at 10^dec_min, end plot at 10^dec_max [rad/s] - [dec_min, dec_max] = bode_bounds (zer, pol, digital, t_sam); - - n_freq = 1000; # Number of frequencies evaluated for plotting - - w = logspace (dec_min, dec_max, n_freq); # [rad/s] - endif - - if (digital) # Discrete system - s = exp (i * w * t_sam); - else # Continuous system - s = i * w; - endif - - l_s = length (s); - sigma_min = zeros (1, l_s); - sigma_max = zeros (1, l_s); - - switch (ptype) - case 0 # Default system - for k = 1 : l_s # Repeat for every frequency s - H = C * inv (s(k)*I - A) * B + D; # Frequency Response Matrix - sigma = svd (H); # Singular Value Decomposition - sigma_min(k) = min (sigma); - sigma_max(k) = max (sigma); - endfor - - case 1 # Inversed system - for k = 1 : l_s - H = inv (C * inv (s(k)*I - A) * B + D); - sigma = svd (H); - sigma_min(k) = min (sigma); - sigma_max(k) = max (sigma); - endfor - - case 2 # Inversed sensitivity - for k = 1 : l_s - H = J + C * inv (s(k)*I - A) * B + D; - sigma = svd (H); - sigma_min(k) = min (sigma); - sigma_max(k) = max (sigma); - endfor - - case 3 # Inversed complementary sensitivity - for k = 1 : l_s - H = J + inv (C * inv (s(k)*I - A) * B + D); - sigma = svd (H); - sigma_min(k) = min (sigma); - sigma_max(k) = max (sigma); - endfor - endswitch - - if (nargout == 0) # Plot the information - - ## Convert to dB for plotting - sigma_min_db = 20 * log10 (sigma_min); - sigma_max_db = 20 * log10 (sigma_max); - - ## Determine xlabel - if (digital) - xl_str = sprintf ('Frequency [rad/s] Pi / T = %g', pi/t_sam); - else - xl_str = 'Frequency [rad/s]'; - endif - - ## Plot results - semilogx (w, sigma_min_db, 'b', w, sigma_max_db, 'b') - title ('Singular Values') - xlabel (xl_str) - ylabel ('Singular Values [dB]') - grid on - else # Return values - sigma_min_r = sigma_min; - sigma_max_r = sigma_max; - w_r = w; - endif - -endfunction - - -%!shared sigma_min_exp, sigma_max_exp, w_exp, sigma_min_obs, sigma_max_obs, w_obs -%! A = [1 2; 3 4]; -%! B = [5 6; 7 8]; -%! C = [4 3; 2 1]; -%! D = [8 7; 6 5]; -%! w = [2 3]; -%! sigma_min_exp = [0.698526948925716 0.608629874340667]; -%! sigma_max_exp = [7.91760889977901 8.62745836756994]; -%! w_exp = [2 3]; -%! [sigma_min_obs, sigma_max_obs, w_obs] = svplot (ss (A, B, C, D), w); -%!assert (sigma_min_obs, sigma_min_exp, 4*eps); # tolerance manually tweaked -%!assert (sigma_max_obs, sigma_max_exp, 12*eps); # tolerance manually tweaked -%!assert (w_obs, w_exp, 2*eps); +## Copyright (C) 2009 Lukas Reichlin +## +## 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 ;. + +## -*- texinfo -*- +## @deftypefn{Function File} {[@var{sigma_min}, @var{sigma_max}, @var{w}] =} svplot (@var{sys}) +## @deftypefnx{Function File} {[@var{sigma_min}, @var{sigma_max}, @var{w}] =} svplot (@var{sys}, @var{w}) +## @deftypefnx{Function File} {[@var{sigma_min}, @var{sigma_max}, @var{w}] =} svplot (@var{sys}, @var{[]}, @var{ptype}) +## If no output arguments are given, the singular value plot of a MIMO +## system is printed on the screen; +## otherwise, the singular values of the system data structure are +## computed and returned. +## +## @strong{Inputs} +## @table @var +## @item sys +## System data structure. Must be either purely continuous or discrete; +## see @code{is_digital}. +## @item w +## Optional vector of frequency values. If @var{w} is not specified, it +## will be calculated by @code{bode_bounds}. +## @item ptype = 0 +## Singular values of the frequency response H of system sys. Default Value. +## @item ptype = 1 +## Singular values of the frequency response inv(H); i.e. inversed system. +## @item ptype = 2 +## Singular values of the frequency response I + H; i.e. inversed sensitivity +## (or return difference) if H = P * C. +## @item ptype = 3 +## Singular values of the frequency response I + inv(H); i.e. inversed complementary +## sensitivity if H = P * C. +## @end table +## +## @strong{Outputs} +## @table @var +## @item sigma_min +## Vector of minimal singular values. +## @item sigma_max +## Vector of maximal singular values. +## @item w +## Vector of frequency values used. +## @end table +## +## @seealso{bode, svd, bode_bounds, is_digital} +## @end deftypefn + +## Author: Lukas Reichlin +## Version: 0.3 + +function [sigma_min_r, sigma_max_r, w_r] = svplot (sys, w, ptype) + + ## Check whether arguments are OK + if (nargin < 1 || nargin > 3) + print_usage (); + endif + + if(! isstruct (sys)) + error ("svplot: first argument sys must be a system data structure"); + endif + + if (nargin == 1) + w = []; + elseif (! isvector (w) && ! isempty (w)) + error ("svplot: second argument w must be a vector of frequencies"); + endif + + if (is_siso (sys)) + warning ("svplot: sys is a SISO system. You might want to use bode(sys) instead."); + endif + + ## Get state space matrices + sys = sysupdate (sys, "ss"); + [A, B, C, D] = sys2ss (sys); + I = eye (size (A)); + + ## Get system information + digital = is_digital(sys, 2); + t_sam = sysgettsam(sys); + + ## Error for mixed systems + if (digital == -1) + error ("svplot: system must be either purely continuous or purely discrete"); + endif + + ## Handle plot type + if (nargin == 3) + if (isfloat (ptype)) # Numeric constants like 2 are NOT integers in Octave! + if (ptype != 0 && ptype != 1 && ptype != 2 && ptype != 3) + error ("svplot: third argument ptype must be 0, 1, 2 or 3"); + endif + else + error ("svplot: third argument ptype must be a number"); + endif + + [n_c, n_d, m, p] = sysdimensions (sys); + if (m != p) + error ("svplot: system must be square for ptype 1, 2 or 3"); + endif + J = eye(m); + else + ptype = 0; # Default value + endif + + ## Find interesting frequency range w if not specified + if (isempty (w)) + + ## Since no frequency vector w has been specified, the interesting + ## decades of the sigma plot have to be found. The already existing + ## function bode_bounds does exactly that. + + zer = tzero (sys); + pol = eig (A); + + ## Begin plot at 10^dec_min, end plot at 10^dec_max [rad/s] + [dec_min, dec_max] = bode_bounds (zer, pol, digital, t_sam); + + n_freq = 1000; # Number of frequencies evaluated for plotting + + w = logspace (dec_min, dec_max, n_freq); # [rad/s] + endif + + if (digital) # Discrete system + s = exp (i * w * t_sam); + else # Continuous system + s = i * w; + endif + + l_s = length (s); + sigma_min = zeros (1, l_s); + sigma_max = zeros (1, l_s); + + switch (ptype) + case 0 # Default system + for k = 1 : l_s # Repeat for every frequency s + H = C * inv (s(k)*I - A) * B + D; # Frequency Response Matrix + sigma = svd (H); # Singular Value Decomposition + sigma_min(k) = min (sigma); + sigma_max(k) = max (sigma); + endfor + + case 1 # Inversed system + for k = 1 : l_s + H = inv (C * inv (s(k)*I - A) * B + D); + sigma = svd (H); + sigma_min(k) = min (sigma); + sigma_max(k) = max (sigma); + endfor + + case 2 # Inversed sensitivity + for k = 1 : l_s + H = J + C * inv (s(k)*I - A) * B + D; + sigma = svd (H); + sigma_min(k) = min (sigma); + sigma_max(k) = max (sigma); + endfor + + case 3 # Inversed complementary sensitivity + for k = 1 : l_s + H = J + inv (C * inv (s(k)*I - A) * B + D); + sigma = svd (H); + sigma_min(k) = min (sigma); + sigma_max(k) = max (sigma); + endfor + endswitch + + if (nargout == 0) # Plot the information + + ## Convert to dB for plotting + sigma_min_db = 20 * log10 (sigma_min); + sigma_max_db = 20 * log10 (sigma_max); + + ## Determine xlabel + if (digital) + xl_str = sprintf ('Frequency [rad/s] Pi / T = %g', pi/t_sam); + else + xl_str = 'Frequency [rad/s]'; + endif + + ## Plot results + semilogx (w, sigma_min_db, 'b', w, sigma_max_db, 'b') + title ('Singular Values') + xlabel (xl_str) + ylabel ('Singular Values [dB]') + grid on + else # Return values + sigma_min_r = sigma_min; + sigma_max_r = sigma_max; + w_r = w; + endif + +endfunction + + +%!shared sigma_min_exp, sigma_max_exp, w_exp, sigma_min_obs, sigma_max_obs, w_obs +%! A = [1 2; 3 4]; +%! B = [5 6; 7 8]; +%! C = [4 3; 2 1]; +%! D = [8 7; 6 5]; +%! w = [2 3]; +%! sigma_min_exp = [0.698526948925716 0.608629874340667]; +%! sigma_max_exp = [7.91760889977901 8.62745836756994]; +%! w_exp = [2 3]; +%! [sigma_min_obs, sigma_max_obs, w_obs] = svplot (ss (A, B, C, D), w); +%!assert (sigma_min_obs, sigma_min_exp, 4*eps); # tolerance manually tweaked +%!assert (sigma_max_obs, sigma_max_exp, 12*eps); # tolerance manually tweaked +%!assert (w_obs, w_exp, 2*eps); This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. ```
 [Octave-cvsupdate] SF.net SVN: octave:[6114] trunk/octave-forge/main/control/inst/svplot.m From: - 2009-08-17 09:51:19 ```Revision: 6114 http://octave.svn.sourceforge.net/octave/?rev=6114&view=rev Author: paramaniac Date: 2009-08-17 09:51:10 +0000 (Mon, 17 Aug 2009) Log Message: ----------- svplot: little coding style fixes Modified Paths: -------------- trunk/octave-forge/main/control/inst/svplot.m Modified: trunk/octave-forge/main/control/inst/svplot.m =================================================================== --- trunk/octave-forge/main/control/inst/svplot.m 2009-08-17 08:56:33 UTC (rev 6113) +++ trunk/octave-forge/main/control/inst/svplot.m 2009-08-17 09:51:10 UTC (rev 6114) @@ -85,8 +85,8 @@ I = eye (size (A)); ## Get system information - digital = is_digital(sys, 2); - t_sam = sysgettsam(sys); + digital = is_digital (sys, 2); + t_sam = sysgettsam (sys); ## Error for mixed systems if (digital == -1) @@ -107,7 +107,7 @@ if (m != p) error ("svplot: system must be square for ptype 1, 2 or 3"); endif - J = eye(m); + J = eye (m); else ptype = 0; # Default value endif This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. ```
 [Octave-cvsupdate] SF.net SVN: octave:[6156] trunk/octave-forge/main/control/inst/svplot.m From: - 2009-08-27 11:34:10 ```Revision: 6156 http://octave.svn.sourceforge.net/octave/?rev=6156&view=rev Author: paramaniac Date: 2009-08-27 11:33:52 +0000 (Thu, 27 Aug 2009) Log Message: ----------- svplot: fix broken test Modified Paths: -------------- trunk/octave-forge/main/control/inst/svplot.m Modified: trunk/octave-forge/main/control/inst/svplot.m =================================================================== --- trunk/octave-forge/main/control/inst/svplot.m 2009-08-27 09:28:27 UTC (rev 6155) +++ trunk/octave-forge/main/control/inst/svplot.m 2009-08-27 11:33:52 UTC (rev 6156) @@ -218,6 +218,7 @@ %! sigma_max_exp = [7.91760889977901 8.62745836756994]; %! w_exp = [2 3]; %! [sigma_min_obs, sigma_max_obs, w_obs] = svplot (ss (A, B, C, D), w); -%!assert (sigma_min_obs, sigma_min_exp, 4*eps); # tolerance manually tweaked -%!assert (sigma_max_obs, sigma_max_exp, 12*eps); # tolerance manually tweaked +%!assert (sigma_min_obs, sigma_min_exp, 8*eps); # tolerance manually tweaked +%!assert (sigma_max_obs, sigma_max_exp, 16*eps); # tolerance manually tweaked %!assert (w_obs, w_exp, 2*eps); + This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. ```