My first contribution to Octave

[paramaniac]

Dear Octave Community

I wish to make my first contribution to Octave, intended for the  
control package. It's a function called svplot(sys) and it plots the  
singular values of a MIMO system over a range of frequencies, similar  
to Matlab's sigma(sys). You can find the m-file together with a little  
demo here:

http://n.ethz.ch/student/lukasre/download/svplot/

I hope my work is useful to Octave :-)

Regards,
Lukas

BTW: I'm an engineering student, not a professional coder, so please  
let me know about possibilities to improve the function!
_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[paramaniac]

On Fri, 2009-07-03 at 20:28 +0200, Lukas Reichlin wrote:
How about help?

Should we reference Matlab in comments or help files???

Here is a suggested diff:

tomdean

*** svplot.m~ Fri Jul  3 13:05:27 2009
--- svplot.m Fri Jul  3 13:17:58 2009
***************
*** 15,21 ****


 function [sigma_min_r, sigma_max_r, w_r] = svplot (sys, w)
!
 if (nargin < 1 || nargin > 2)
     print_usage();
 endif
--- 15,33 ----


 function [sigma_min_r, sigma_max_r, w_r] = svplot (sys, w)
!   ## Usage: [sigma_min_r, sigma_max_r, w_r] = svplot (sys)
!   ##  -or-  [sigma_min_r, sigma_max_r, w_r] = svplot (sys, w)
!   ## Where:
!   ##       sys is a MIMO system
!   ##       w is an optional frequency vector
!   ##
!   ## If w is not specified, it will be calculated.
!   ##
!   ## Plot the singular values of a MIMO system over a range of
!   ## frequencies, similar to Matlab's sigma(sys)
!   ##
!   ## Lukas Reichlin <lukas.reichlin@...>
!  
 if (nargin < 1 || nargin > 2)
     print_usage();
 endif

You are right, I forgot the most important feature of every software: documentation :-)

For the Package Reference:

[sigma_min, sigma_max, w] = svplot (sys, 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.

Otherwise, the singular values of the system data structure are computed and returned.

Inputs

sys system data structure (must be either purely continuous or discrete; see is_digital)

w Optional vector of frequency values. If w is not specified, it will be calculated by bode_bounds

Outputs

sigma_min Vector of minimal singular values

sigma_max Vector of maximal singular values

w Vector of frequency values used

For the file itself:

 function [sigma_min, sigma_max, w] = svplot (sys, w)
!   ## Usage: svplot (sys)

!   ##  -or-   svplot (sys, w)
!   ## Where:
!   ##       sys is a MIMO system
!   ##       w is an optional frequency vector
!   ##
!   ## If w is not specified, it will be calculated by bode_bounds.
!   ##
!   ## If no output arguments are given:

!   ## The singular value plot of a MIMO system over a range

!   ## of frequencies is printed on the screen.

!   ## Otherwise, the singular values of the system data structure

!   ## are computed and returned

!   ##

!   ## Lukas Reichlin <lukas.reichlin@...>
!  

No, I don't think we should reference Matlab, because it would be unusual in the OCST Manual.

Regards,

Lukas

_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[Luca Favatella]

On 04/07/2009, Lukas Reichlin <lukas.reichlin@...> wrote:
[...]
> [sigma_min, sigma_max, w] = svplot (sys, w)
[...]
> Regards,
> Lukas

Please see this commit
http://octave.svn.sourceforge.net/viewvc/octave?view=rev&revision=5986

I added the file svplot.m refactoring it a bit (doc and coding style),
and added it to the index.
Please tell me if it is ok or if I forgot something.

Can you please post some tests for this function?
(1 or 2 will suffice, even simple ones, without involving graphical output)
You can use this file as a reference (see the last lines)
http://octave.svn.sourceforge.net/viewvc/octave/trunk/octave-forge/main/control/inst/place.m?revision=5857&view=markup


Please consider that I am only an engineering student and I don't
understand very well what your code is doing.


Thanks for all your work,
Luca Favatella
_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[Luca Favatella]

On 06/07/2009, Luca Favatella <slackydeb@...> wrote:

> On 04/07/2009, Lukas Reichlin <lukas.reichlin@...> wrote:
> [...]
>> [sigma_min, sigma_max, w] = svplot (sys, w)
> [...]
>> Regards,
>> Lukas
>
> Please see this commit
> http://octave.svn.sourceforge.net/viewvc/octave?view=rev&revision=5986
>
> I added the file svplot.m refactoring it a bit (doc and coding style),
> and added it to the index.
> Please tell me if it is ok or if I forgot something.

Little coding style fixes committed:
http://octave.svn.sourceforge.net/viewvc/octave?view=rev&revision=5987
_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[paramaniac]

On 06.07.2009, at 10:38, Luca Favatella wrote:

Please see this commit
http://octave.svn.sourceforge.net/viewvc/octave?view=rev&revision=5986

I added the file svplot.m refactoring it a bit (doc and coding style),
and added it to the index.
Please tell me if it is ok or if I forgot something.

Can you please post some tests for this function?
(1 or 2 will suffice, even simple ones, without involving graphical output)
You can use this file as a reference (see the last lines)
http://octave.svn.sourceforge.net/viewvc/octave/trunk/octave-forge/main/control/inst/place.m?revision=5857&view=markup

Thanks for cleaning up things! I have to figure out first how these tests work. Are they executed automatically in some way? Then I need to think of some good examples. Additionally, I shouldn't forget testing with discrete models. I don't know if it works yet. 

Please consider that I am only an engineering student and I don't
understand very well what your code is doing.


Thanks for all your work,
Luca Favatella

Yet another student :-) I will try to improve comments for better understanding of the code.

Regards,

Lukas

_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[Luca Favatella]

[sorry for the double email, I forgot to keep cc]


On 06/07/2009, Lukas Reichlin <lukas.reichlin@...> wrote:
> On 06.07.2009, at 10:38, Luca Favatella wrote:
[...]
>> Can you please post some tests for this function?
>> (1 or 2 will suffice, even simple ones, without involving graphical
>> output)
>> You can use this file as a reference (see the last lines)
>> http://octave.svn.sourceforge.net/viewvc/octave/trunk/octave-forge/main/control/inst/place.m?revision=5857&view=markup
>
> Thanks for cleaning up things! I have to figure out first how these
> tests work. Are they executed automatically in some way? Then I need

All tests of a function are executed when you type "test function_name".


> to think of some good examples. Additionally, I shouldn't forget
> testing with discrete models. I don't know if it works yet.
[...]
> Yet another student :-) I will try to improve comments for better
> understanding of the code.

Good :).


> Regards,
> Lukas
_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[paramaniac]

I improved the code of svplot and added a test procedure. Although I  
got a little problem with the assert command. I don't know how to use  
it for several return values.

Regards
Lukas




## Copyright (C) 2009 Lukas Reichlin <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 <http://www.gnu.org/licenses/>.

## -*- 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})
## If no output arguments are given, the singular value plot of a MIMO
## system over a range of frequencies 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}.
## @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{is_digital}
## @end deftypefn

## Author: Lukas Reichlin <lukas.reichlin@...>
## Version: 0.1alpha


function [sigma_min_r, sigma_max_r, w_r] = svplot (sys, w)

   ## Check whether arguments are OK
   if (nargin < 1 || nargin > 2)
     print_usage ();
   endif

   if(! isstruct (sys))
     error ("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");
     endif
   endif

   ## Get State Space Matrices
   sys = sysupdate (sys, "ss");
   [A, B, C, D] = sys2ss (sys);
   I = eye (size (A));


   ## Find interesting Frequency Range w if not specified
   if (nargin < 2)
     ## 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.

     j = 1; # Index Counter

     for m = 1 : size (B, 2) # For all Inputs
       for p = 1 : size (C, 1) # For all Outputs

         sysp = sysprune (sys, p, m);
         [zer, pol, k, tsam, inname, outname] = sys2zp (sysp);

         if (tsam == 0)
           DIGITAL = 0;
         else
           DIGITAL = 1;
           ## Discrete systems not yet tested!
         endif

         [wmin, wmax] = bode_bounds (zer, pol, DIGITAL, tsam);

         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
     n_steps = 1000; # Number of Frequencies evaluated for Plotting

     w = logspace (dec_min, dec_max, n_steps); # rad/s
   endif


   ## Repeat for every Frequency w
   for k = 1 : size (w, 2)

     ## Frequency Response Matrix
     P = C * inv (i*w(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
     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]')
     grid on
   else # Return Values
     sigma_min_r = sigma_min;
     sigma_max_r = sigma_max;
     w_r = w;
   endif

endfunction


%!shared A, B, C, D, w, sigma_min_exp, sigma_max_exp, w_exp
%! 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];
%! assert (svplot (ss (A, B, C, D), w), {sigma_min_exp, sigma_max_exp,  
w_exp}, 2*eps);


_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[paramaniac]

The mailing list accidentally made some line breaks. Here is the  
intact version of my code:

http://n.ethz.ch/student/lukasre/download/svplot/

Sorry for the noise.

Regards,
Lukas
_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[Luca Favatella]

On 07/07/2009, Lukas Reichlin <lukas.reichlin@...> wrote:
[...]
> http://n.ethz.ch/student/lukasre/download/svplot/

Committed, with little changes. Please see
http://octave.svn.sourceforge.net/viewvc/octave/trunk/octave-forge/main/control/inst/svplot.m?revision=5989

About the test, I refactored it to deal better with multiple return variables.
I had to manually tweak tolerances on my pc. Please take a look at it.

If you can, feel free to propose other tests (even on other functions,
even simple ones). They are greatly appreciated.


[...]
> Regards,
> Lukas

Thanks for all your work,
Luca Favatella
_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[paramaniac]

On 07.07.2009, at 20:32, Luca Favatella wrote:

> On 07/07/2009, Lukas Reichlin <lukas.reichlin@...> wrote:
> [...]
>> http://n.ethz.ch/student/lukasre/download/svplot/
>
> Committed, with little changes. Please see
> http://octave.svn.sourceforge.net/viewvc/octave/trunk/octave-forge/main/control/inst/svplot.m?revision=5989
>
> About the test, I refactored it to deal better with multiple return  
> variables.
> I had to manually tweak tolerances on my pc. Please take a look at it.
>
> If you can, feel free to propose other tests (even on other functions,
> even simple ones). They are greatly appreciated.
>
>
> [...]
>> Regards,
>> Lukas
>
> Thanks for all your work,
> Luca Favatella

Thanks for your work and patience :-)

I made some cosmetic changes to improve (hopefully) readability of the  
code.

http://n.ethz.ch/student/lukasre/download/svplot/

Do you know how discrete systems work? Until now, we tackled only  
continuous systems in our control systems lectures at university.  
Therefore I'm not sure if my code works for discrete systems too.

BTW: I've just created my Source Forge account. My user name is  
paramaniac (I'm addicted to paragliding ;-). I realized that MacOSX  
has all the necessary tools (svn, ssh, ...) already on board. That's  
convenience! Do you think it's safe when I update my code by myself?  
This could save you some work, although I don't want to mess things up!

Regards,
Lukas
_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[Luca Favatella]

On 07/07/2009, Lukas Reichlin <lukas.reichlin@...> wrote:
[...]
> I made some cosmetic changes to improve (hopefully) readability of the
> code.
>
> http://n.ethz.ch/student/lukasre/download/svplot/

Good.
I didn't review/commit it.
Please read below.


> Do you know how discrete systems work? Until now, we tackled only

I studied them, if I remember well :P.


> continuous systems in our control systems lectures at university.
> Therefore I'm not sure if my code works for discrete systems too.

I think the main question is:
is there someone (at least one, e.g. you) out there using this
function for discrete systems?
If even you don't need discrete systems support in the function (at
least at the moment), I think the right thing to do is disabling
discrete systems support (at least at the moment).


> BTW: I've just created my Source Forge account. My user name is
> paramaniac (I'm addicted to paragliding ;-). I realized that MacOSX
> has all the necessary tools (svn, ssh, ...) already on board. That's
> convenience! Do you think it's safe when I update my code by myself?

Yours is a great great great idea!

I'm not the author of the control package.
I only contributed a few fixes, tests and a pair of new functions, as
you are doing :).
As a mantainer, my main goal is (if not able to actively improve the
package) at least reviewing and merging community contributions.
If who contributes also commits, from my POV is the best :).


So, I think it is good you email octave-dev (clean thread) asking for
commit access on octave-forge svn, saying that you are interested in
committing and improving your contributions in the control package.
I hope they'll say to add you as a mantainer/co-mantainer of control
package , or similar.


> This could save you some work, although I don't want to mess things up!

Don't worry, you are improving :).
You can commit and I could review your commits when you ask me.


> Regards,
> Lukas

Cheers,
Luca Favatella
_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave