Mercurial > hg > octave-max
changeset 3236:98e15955107e
[project @ 1999-03-05 07:17:10 by jwe]
line wrap: on
line diff
--- a/scripts/control/DEMOcontrol.m +++ b/scripts/control/DEMOcontrol.m @@ -21,7 +21,6 @@ # Demo programs: bddemo.m, frdemo.m, analdemo.m, moddmeo.m, rldemo.m # # Written by David Clem August 15, 1994 -# $Revision: 2.0.0.2 $ disp(' O C T A V E C O N T R O L S Y S T E M S T O O L B O X')
--- a/scripts/control/abcddim.m +++ b/scripts/control/abcddim.m @@ -35,54 +35,6 @@ # Written by A. S. Hodel (scotte@eng.auburn.edu) August 1993. # a s hodel: modified to accept pure-gain systems aug 1996 -# $Revision: 1.17 $ -# $Log: abcddim.m,v $ -# Revision 1.17 1998-12-10 03:06:31 jwe -# *** empty log message *** -# -# Revision 2.0.0.2 1998/12/08 23:29:21 hodel -# Octave-Marsyas Interface updated for signals-as-lists -# -# Revision 2.0.0.1 1998/12/08 21:40:44 hodel -# Dummy version to match ftp.eng.auburn.edu version number -# -# Revision 2.0.0.0 1998/12/08 21:36:51 hodel -# Branch for beta release patches -# -# Revision 2.0 1998/12/08 21:34:56 hodel -# Initial beta release of signals-as-lists rewrite; -# sysdimensions now takes opt as an argument -# -# Revision 2.0.0.1 1998/12/08 20:54:18 hodel -# sysdimensions takes opt parameter now -# -# Revision 2.0.0.0 1998/12/08 20:30:08 hodel -# beta release revision -# -# Revision 2.0 1998/12/08 20:27:55 hodel -# Initial list rewrite of OCST -# -# Revision 1.2 1998/10/05 17:12:56 hodelas -# various bug changes -# -# Revision 1.1.1.1 1998/05/19 20:24:05 jwe -# -# Revision 1.4 1997/12/01 16:44:22 scotte -# *** empty log message *** -# -# Revision 1.3 1997/02/12 15:38:14 hodel -# *** empty log message *** -# -# -# fixed typo -# -# Revision 1.1 1997/02/12 11:34:53 hodel -# Initial revision -# -# Revision 1.7 1997/02/07 15:29:56 scotte -# fixed is_square check to allow for empty a matrix -# (this allows for pure gain blocks) -# if (nargin != 4) error ("abcddim: four arguments required");
--- a/scripts/control/abcddims.m +++ b/scripts/control/abcddims.m @@ -24,7 +24,6 @@ # get set to 0. my and ny are the row and column dimensions of the result. # Written by A. S. Hodel (scotte@eng.auburn.edu) Feb 1997 -# $Revision: 2.0.0.2 $ y = x; if(isempty(y))
--- a/scripts/control/analdemo.m +++ b/scripts/control/analdemo.m @@ -20,7 +20,6 @@ # Octave Controls toolbox demo: State Space analysis demo # Written by David Clem August 15, 1994 # Updated by John Ingram December 1996 -# $Revision: 2.0.0.2 $ while (1) clc
--- a/scripts/control/are.m +++ b/scripts/control/are.m @@ -35,7 +35,6 @@ # See also: balance # Written by A. S. Hodel (scotte@eng.auburn.edu) August 1993. -# $Revision: 1.20 $ if (nargin == 3 || nargin == 4) if (nargin == 4)
--- a/scripts/control/bddemo.m +++ b/scripts/control/bddemo.m @@ -21,17 +21,6 @@ # Written by David Clem August 15, 1994 # Modified by A S Hodel Summer-Fall 1996 -# Revision 1.6 1998/05/05 17:02:45 scotte -# updated May 5 by Kai Meuller -# -# Revision 1.2 1998/05/05 13:21:51 mueller -# buildssic command description -# buildssic example added -# -# Revision 1.1 1998/05/05 13:20:26 mueller -# Initial revision -# - str_sav = implicit_str_to_num_ok; sav_page = page_screen_output; implicit_str_to_num_ok = 1;
--- a/scripts/control/bode.m +++ b/scripts/control/bode.m @@ -57,64 +57,6 @@ # Modified by David Clem November 13, 1994 # again by A. S. Hodel July 1995 (smart plot range, etc.) # Modified by Kai P. Mueller September 28, 1997 (multiplot mode) -# $Revision: 2.0.0.2 $ -# $Log: bode.m,v $ -# Revision 2.0.0.2 1998/12/08 23:29:23 hodel -# Octave-Marsyas Interface updated for signals-as-lists -# -# Revision 2.0.0.1 1998/12/08 21:40:45 hodel -# Dummy version to match ftp.eng.auburn.edu version number -# -# Revision 2.0.0.0 1998/12/08 21:36:52 hodel -# Branch for beta release patches -# -# Revision 2.0 1998/12/08 21:34:57 hodel -# Initial beta release of signals-as-lists rewrite; -# sysdimensions now takes opt as an argument -# -# Revision 2.0.0.1 1998/12/08 20:54:19 hodel -# sysdimensions takes opt parameter now -# -# Revision 2.0.0.0 1998/12/08 20:30:09 hodel -# beta release revision -# -# Revision 2.0 1998/12/08 20:27:56 hodel -# Initial list rewrite of OCST -# -# Revision 1.7 1998/10/21 12:46:59 hodelas -# moved grid command so that grid appears in plots -# -# Revision 1.6 1998/09/04 20:57:18 hodelas -# fixed bodquist bug (use reshape instead of transpose); removed extraneous -# output from bode. -# -# Bodquist is now much faster -# -# Revision 1.4 1998/08/24 15:50:03 hodelas -# updated documentation -# -# Revision 1.3 1998/08/13 20:11:27 hodelas -# Added calls to axis2dlim for flat-plots -# -# Revision 1.2 1998/07/24 18:16:51 hodelas -# rewrote bodquist as a function call. nyquist interactive plot now optional -# -# Revision 1.1.1.1 1998/05/19 20:24:05 jwe -# -# Revision 1.7 1998/02/09 13:04:11 scotte -# fixed oneplot/gset nokey to function only if gnuplot_has_multiplot -# -# Revision 1.6 1997/12/01 16:51:50 scotte -# updated by Mueller 27 Nov 97 -# -# Revision 1.2 1997/11/24 18:53:01 mueller -# gset autoscale prevents the following error message: -# line 0: x range must be greater than 0 for log scale! -# gset nokey and call to oneplot() added -# -# Revision 1.1 1997/11/11 17:31:27 mueller -# Initial revision -# # check number of input arguments given if (nargin < 1 | nargin > 4) @@ -158,7 +100,6 @@ tistr = "(jw)"; endif xlabel(xlstr); - ylabel("Gain in dB"); if(is_siso(sys)) if (gnuplot_has_multiplot) subplot(2,1,1); @@ -172,7 +113,13 @@ disp(outlist(outname," ")); endif wv = [min(w), max(w)]; - md = 20*log10(mag); + if(max(mag) > 0) + ylabel("Gain in dB"); + md = 20*log10(mag); + else + ylabel("Gain |Y/U|") + md = mag; + endif axvec = axis2dlim([vec(w),vec(md)]); axvec(1:2) = wv;
--- a/scripts/control/bode_bounds.m +++ b/scripts/control/bode_bounds.m @@ -24,8 +24,6 @@ # # used internally in freqresp -# $Revision: 2.0.0.2 $ - # make sure zer,pol are row vectors if(!isempty(pol)) pol = reshape(pol,1,length(pol)); endif if(!isempty(zer)) zer = reshape(zer,1,length(zer)); endif
--- a/scripts/control/buildssic.m +++ b/scripts/control/buildssic.m @@ -108,7 +108,6 @@ # # Written by Kai Mueller April 1998 -# $Revision: 2.0.0.2 $ if((nargin < 5) || (nargin > 12)) usage("[sys] = buildssic(Clst,Ulst,Olst,Ilst,s1,s2,s3,s4,s5,s6,s7,s8)");
--- a/scripts/control/c2d.m +++ b/scripts/control/c2d.m @@ -49,7 +49,6 @@ # Written by R.B. Tenison (btenison@eng.auburn.edu) # October 1993 # Updated by John Ingram for system data structure August 1996 -# $Revision: 1.15 $ save_val = implicit_str_to_num_ok; # save for later implicit_str_to_num_ok = 1;
--- a/scripts/control/com2str.m +++ b/scripts/control/com2str.m @@ -25,7 +25,6 @@ # 0 (default): -1, 0, 1, 1i, 1 + 0.5i # 1 (for use with zpout): -1, 0, + 1, + 1i, + 1 + 0.5i # -# $Revision: 2.0.0.2 $ if (nargin < 1 | nargin > 2) usage("com2str(zz{,flg})");
--- a/scripts/control/controldemo.m +++ b/scripts/control/controldemo.m @@ -21,7 +21,6 @@ # Demo programs: bddemo.m, frdemo.m, analdemo.m, moddmeo.m, rldemo.m # # Written by David Clem August 15, 1994 -# $Revision: 2.0.0.2 $ disp(' O C T A V E C O N T R O L S Y S T E M S T O O L B O X')
--- a/scripts/control/ctrb.m +++ b/scripts/control/ctrb.m @@ -35,7 +35,6 @@ # Written by Kai P. Mueller November 4, 1997 # based on is_controllable.m of Scottedward Hodel # modified by - # $Revision: 2.0.0.2 $ if (nargin == 2) a = sys;
--- a/scripts/control/d2c.m +++ b/scripts/control/d2c.m @@ -53,7 +53,6 @@ # Written by R. Bruce Tenison August 23, 1994 # Updated by John Ingram for system data structure August 1996 # SYS_INTERNAL accesses members of system data structure -# $Revision: 2.0.0.2 $ save_val = implicit_str_to_num_ok; # save for later implicit_str_to_num_ok = 1;
--- a/scripts/control/damp.m +++ b/scripts/control/damp.m @@ -29,7 +29,6 @@ # Written by Kai P. Mueller September 29, 1997. # Update -# $Revision: 2.0.0.2 $ # assume a continuous system DIGITAL = 0;
--- a/scripts/control/dare.m +++ b/scripts/control/dare.m @@ -47,7 +47,6 @@ ## Author: A. S. Hodel <scotte@eng.auburn.edu> ## Created: August 1993 ## Adapted-By: jwe -## $Revision: 1.17 $ function x = dare (a, b, c, r, opt)
--- a/scripts/control/dcgain.m +++ b/scripts/control/dcgain.m @@ -26,7 +26,6 @@ # # Written by Kai P Mueller (mueller@ifr.ing.tu-bs.de) October 1, 1997 -# $Revision: 2.0.0.2 $ if((nargin < 1) || (nargin > 2) || (nargout > 1)) usage("[gm, ok] = dcgain(sys[, tol])");
--- a/scripts/control/dezero.m +++ b/scripts/control/dezero.m @@ -26,7 +26,6 @@ ## Adapted from deblank by A. S. Hodel (a.s.hodel@eng.auburn.edu) ## (the name dezero is a reference to the Fermilab D0 experiment, ## where my sister did her PhD research) -## $Revision: 2.0.0.2 $ function t = dezero (s)
--- a/scripts/control/dgkfdemo.m +++ b/scripts/control/dgkfdemo.m @@ -19,7 +19,6 @@ function dgkfdemo() # Octave Controls toolbox demo: H2/Hinfinity options demos # Written by A. S. Hodel June 1995 -# $Revision: 2.0.0.2 $ save_val = page_screen_output; page_screen_output = 1;
--- a/scripts/control/dgram.m +++ b/scripts/control/dgram.m @@ -25,7 +25,6 @@ # a m a' - m + b*b' = 0 # Written by A. S. Hodel July 1995 - # $Revision: 1.13 $ # let dlyap do the error checking... m = dlyap(a,b*b');
--- a/scripts/control/dhinfdemo.m +++ b/scripts/control/dhinfdemo.m @@ -1,3 +1,21 @@ +# Copyright (C) 1996,1998 Kai Mueller +# +# This file is part of Octave. +# +# Octave 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 2, or (at your option) any +# later version. +# +# Octave 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 Octave; see the file COPYING. If not, write to the Free +# Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. + # ------------------------------------------------------------ # dhinfdemo Design of a discrete H_infinity controller. # This is not a true discrete design. The design @@ -40,8 +58,6 @@ # K. Mueller, <mueller@ifr.ing.tu-bs.de> # Technical University of Braunschweig, IfR -# $Revision: 2.0.0.2 $ $Date: 1998/12/08 23:29:27 $ -# echo off disp(" ");
--- a/scripts/control/dlqe.m +++ b/scripts/control/dlqe.m @@ -45,8 +45,6 @@ # Modified for discrete time by R. Bruce Tenison (btenison@eng.auburn.edu) # October, 1993 -# $Revision: 1.16 $ - if (nargin != 5 && nargin != 6) error ("dlqe: invalid number of arguments"); endif
--- a/scripts/control/dlqg.m +++ b/scripts/control/dlqg.m @@ -47,7 +47,6 @@ # See also: lqg, dlqe, dlqr # Written by A. S. Hodel August 1995 -# $Revision: 2.0.0.2 $ warning("dlqg: obsolete. use lqg instead (system data structure format)");
--- a/scripts/control/dlqr.m +++ b/scripts/control/dlqr.m @@ -46,7 +46,6 @@ # Written by A. S. Hodel (scotte@eng.auburn.edu) August 1993. # Converted to discrete time by R. B. Tenison # (btenison@eng.auburn.edu) October 1993 -# $Revision: 1.17 $ if (nargin != 4 && nargin != 5) error ("dlqr: invalid number of arguments");
--- a/scripts/control/dlyap.m +++ b/scripts/control/dlyap.m @@ -34,7 +34,6 @@ # (1977). # Written by A. S. Hodel (scotte@eng.auburn.edu) August 1993. -# $Revision: 1.15 $ if ((n = is_square (a)) == 0) warning ("dlyap: a must be square");
new file mode 100644 --- /dev/null +++ b/scripts/control/dre.m @@ -0,0 +1,126 @@ +# Copyright (C) 1998 A. Scottedward Hodel +# +# This file is part of Octave. +# +# Octave 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 2, or (at your option) any +# later version. +# +# Octave 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 Octave; see the file COPYING. If not, write to the Free +# Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. + +function [tvals,Plist] = dre(sys,Q,R,Qf,t0,tf,Ptol,maxits) +# [tvals,Plist] = dre(sys,Q,R,Qf,t0,tf{,Ptol,maxits}); +# Solve the differential Riccati equation +# -d P/dt = A'P + P A - P B inv(R) B' P + Q +# P(tf) = Qf +# for the LTI system sys. Solution of standard LTI +# state feedback optimization +# min \int_{t_0}^{t_f} x' Q x + u' R u dt + x(t_f)' Qf x(t_f) +# optimal input is +# u = - inv(R) B' P(t) x +# inputs: +# sys: continuous time system data structure +# Q: state integral penalty +# R: input integral penalty +# Qf: state terminal penalty +# t0,tf: limits on the integral +# Ptol: tolerance (used to select time samples; see below); default = 0.1 +# max number of refinement iterations (default=10) +# outputs: +# tvals: time values at which P(t) is computed +# Plist: list values of P(t); nth(Plist,ii) is P(tvals(ii)). +# +# tvals is selected so that || nth(Plist,ii) - nth(Plist,ii-1) || < Ptol +# for ii=2:length(tvals) +# +# Reference: + +if(nargin < 6 | nargin > 8 | nargout != 2) + usage("[tvals,Plist] = dre(sys,Q,R,Qf,t0,tf{,Ptol})"); +elseif(!is_struct(sys)) + error("sys must be a system data structure") +elseif(is_digital(sys)) + error("sys must be a continuous time system") +elseif(!is_matrix(Q) | !is_matrix(R) | !is_matrix(Qf)) + error("Q, R, and Qf must be matrices."); +elseif(!is_scalar(t0) | !is_scalar(tf)) + error("t0 and tf must be scalars") +elseif(t0 >= tf) error("t0=%e >= tf=%e",t0,tf); +elseif(nargin == 6) Ptol = 0.1; +elseif(!is_scalar(Ptol)) error("Ptol must be a scalar"); +elseif(Ptol <= 0) error("Ptol must be positive"); +endif + +if(nargin < 8) maxits = 10; +elseif(!is_scalar(maxits)) error("maxits must be a scalar"); +elseif(maxits <= 0) error("maxits must be positive"); +endif +maxits = ceil(maxits); + +[aa,bb] = sys2ss(sys); +nn = sysdimensions(sys,"cst"); +mm = sysdimensions(sys,"in"); +pp = sysdimensions(sys,"out"); + +if(size(Q) != [nn nn]) + error("Q(%dx%d); sys has %d states",rows(Q),columns(Q),nn); +elseif(size(Qf) != [nn nn]) + error("Qf(%dx%d); sys has %d states",rows(Qf),columns(Qf),nn); +elseif(size(R) != [mm mm]) + error("R(%dx%d); sys has %d inputs",rows(R),columns(R),mm); +endif + +# construct Hamiltonian matrix +H = [aa , -(bb/R)*bb' ; -Q, -aa']; + +# select time step to avoid numerical overflow +fast_eig = max(abs(eig(H))); +tc = log(10)/fast_eig; +nst = ceil((tf-t0)/tc); +tvals = -linspace(-tf,-t0,nst); +Plist = list(Qf); +In = eye(nn); +n1 = nn+1; +n2 = nn+nn; +done = 0; +while(!done) + done = 1; # assume this pass will do the job + # sort time values in reverse order + tvals = -sort(-tvals); + tvlen = length(tvals); + maxerr = 0; + # compute new values of P(t); recompute old values just in case + for ii=2:tvlen + uv_i_minus_1 = [ In ; nth(Plist,ii-1) ]; + delta_t = tvals(ii-1) - tvals(ii); + uv = expm(-H*delta_t)*uv_i_minus_1; + Qi = uv(n1:n2,1:nn)/uv(1:nn,1:nn); + Plist(ii) = (Qi+Qi')/2; + # check error + Perr = norm(nth(Plist,ii) - nth(Plist,ii-1))/norm(nth(Plist,ii)); + maxerr = max(maxerr,Perr); + if(Perr > Ptol) + new_t = mean(tvals([ii,ii-1])); + tvals = [tvals new_t]; + done = 0; + endif + endfor + + # check number of iterations + maxits = maxits - 1; + done = done+(maxits==0); +endwhile +if(maxerr > Ptol) + warning("dre: \n\texiting with%4d points, max rel chg. =%e, Ptol=%e\n", ... + tvlen,maxerr,Ptol); + tvals = tvals(1:length(Plist)); +endif +endfunction
--- a/scripts/control/fir2sys.m +++ b/scripts/control/fir2sys.m @@ -32,7 +32,6 @@ # Name changed to TF2SYS July 1995 # updated for new system data structure format July 1996 # adapted from tf2sys july 1996 - # $Revision: 2.0.0.2 $ save_val = implicit_str_to_num_ok; implicit_str_to_num_ok = 1;
--- a/scripts/control/frdemo.m +++ b/scripts/control/frdemo.m @@ -20,7 +20,6 @@ # Octave Controls toolbox demo: Frequency Response demo # Written by David Clem August 15, 1994 -# $Revision: 2.0.0.2 $ # a s hodel: updated to match new order of ss2zp outputs # J Ingram: updated for system data structure format August 1996
--- a/scripts/control/freqchkw.m +++ b/scripts/control/freqchkw.m @@ -21,7 +21,6 @@ # used by freqresp to check that input frequency vector is legal # A S Hodel July 1996 - # $Revision: 2.0.0.2 $ if(isempty(w)) USEW = 0;
--- a/scripts/control/freqresp.m +++ b/scripts/control/freqresp.m @@ -28,7 +28,6 @@ # ff and w are both returned as row vectors # Written by: R. Bruce Tenison July 11, 1994 - # $Revision: 2.0.0.2 $ # SYS_INTERNAL accesses members of system data structure save_val = empty_list_elements_ok;
--- a/scripts/control/gram.m +++ b/scripts/control/gram.m @@ -25,7 +25,6 @@ # a m + a' + b*b' = 0 # Written by A. S. Hodel - # $Revision: 2.0.0.2 $ # let lyap do the error checking... m = lyap(a,b*b');
--- a/scripts/control/h2norm.m +++ b/scripts/control/h2norm.m @@ -30,7 +30,6 @@ # A. S. Hodel Aug 1995 # updated for system data structure by John Ingram November 1996 - # $Revision: 2.0.0.2 $ if((nargin != 1)) usage("h2gain = h2norm(sys)");
--- a/scripts/control/h2syn.m +++ b/scripts/control/h2syn.m @@ -42,7 +42,6 @@ # Pf: ARE solution matrix for filter subproblem # Updated for System structure December 1996 by John Ingram - # $Revision: 2.0.0.2 $ if ((nargin < 3) | (nargin > 4)) usage("[K,gain, Kc, Kf, Pc, Pf] = h2syn(Asys,nu,ny[,tol])");
--- a/scripts/control/hinf_ctr.m +++ b/scripts/control/hinf_ctr.m @@ -31,11 +31,9 @@ # Do not attempt to use this at home; no argument checking performed. # A. S. Hodel August 1995 - # $Revision: 2.0.0.2 $ # Revised by Kai P Mueller April 1998 to solve the general H_infinity # problem using unitary transformations Q (on w and z) # and non-singular transformations R (on u and y). - # $Revision: 2.0.0.2 $ nw = dgs.nw; nu = dgs.nu;
--- a/scripts/control/hinfdemo.m +++ b/scripts/control/hinfdemo.m @@ -1,3 +1,21 @@ +# Copyright (C) 1996,1998 Kai Mueller +# +# This file is part of Octave. +# +# Octave 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 2, or (at your option) any +# later version. +# +# Octave 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 Octave; see the file COPYING. If not, write to the Free +# Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. + # hinfdemo H_infinity design demos for continuous SISO and MIMO # systems and a discrete system. # The SISO system is difficult to control because @@ -99,7 +117,6 @@ # Kai P. Mueller 30-APR-1998 <mueller@ifr.ing.tu-bs.de -# $Revision: 2.0.0.2 $ yn = []; while (length(yn) < 1)
--- a/scripts/control/hinfnorm.m +++ b/scripts/control/hinfnorm.m @@ -39,7 +39,6 @@ # Iglesias and Glover, "State-Space approach to discrete-time Hinf control," # Int. J. Control, vol 54, #5, 1991 # Zhou, Doyle, Glover, "Robust and Optimal Control," Prentice-Hall, 1996 - # $Revision: 2.0.0.2 $ if((nargin == 0) || (nargin > 4)) usage("[g gmin gmax] = hinfnorm(sys[,tol,gmin,gmax,ptol])");
--- a/scripts/control/hinfsyn.m +++ b/scripts/control/hinfsyn.m @@ -1,4 +1,4 @@ -# Copyright (C) 1996 A. Scottedward Hodel +# Copyright (C) 1996,1998 A. Scottedward Hodel # # This file is part of Octave. # @@ -56,15 +56,10 @@ # A. S. Hodel August 1995 # Updated for Packed system structures December 1996 by John Ingram - # $Revision: 2.0.0.2 $ # # Revised by Kai P Mueller April 1998 to solve the general H_infinity # problem using unitary transformations Q (on w and z) # and non-singular transformations R (on u and y). - # $Revision: 2.0.0.2 $ - - old_page_val = page_screen_output; - page_screen_output = 0; if( (nargin < 1) | (nargin > 8) ) usage("[K,g,GW,Xinf,Yinf] = hinfsyn(Asys,nu,ny,gmin,gmax,gtol,ptol,tol)"); @@ -151,17 +146,38 @@ printf("----------------------------------------"); printf("--------------------------------------\n"); + # set up error messages + errmesg = list(" o o o o o ", ... + " # - - - - ", ... + " o # - - - ", ... + " o o # - - ", ... + " o o o # - ", ... + " o o o o # ", ... + " - - - - - "); + errdesx = list("", ... + "X im eig.", ... + "Hx not Ham.", ... + "X inf.eig", ... + "X not symm.", ... + "X not pos", ... + "R singular"); + + errdesy = list(" ", ... + "Y im eig.", ... + "Hy not Ham.", ... + "Y inf.eig", ... + "Y not symm.", ... + "Y not pos", ... + "Rtilde singular"); + + # now do the search while (!iteration_finished) switch (search_state) - case (0) - g = ghi; - case (1) - g = glo; - case (2) - g = 0.5 * (ghi + glo); - otherwise - error(" *** This should never happen!"); + case (0) g = ghi; + case (1) g = glo; + case (2) g = 0.5 * (ghi + glo); + otherwise error(" *** This should never happen!"); endswitch printf("%10.4f ", g); @@ -175,128 +191,39 @@ Rtilde(1:nz,1:nz) = -g*g*eye(nz); Rtilde = Rtilde + ddot1 * ddot1'; - # build hamiltonian Ha for X_inf - xx = ([BB; -C1'*d1dot]/R) * [d1dot'*C1 BB']; - Ha = [A 0*A; -C1'*C1 -A'] - xx; - x_ha_err = 0; - # copied from from are(...)... - [d, Ha] = balance(Ha); - [u, s] = schur(Ha, "A"); - rev = real(eig(s)); - if (any(abs(rev) <= ptol)) - # eigenvalues near the imaginary axis - x_ha_err = 1; - elseif (sum(rev > 0) != sum(rev < 0)) - # unequal number of positive and negative eigenvalues - x_ha_err = 2; - else - # compute positive Riccati equation solution - u = d * u; - Xinf = u(n+1:2*n,1:n) / u(1:n,1:n); - if (!all(all(finite(Xinf)))) - x_ha_err = 3; - elseif (norm(Xinf-Xinf') >= 10*ptol) - # solution not symmetric - x_ha_err = 4; - else - # positive semidefinite? - rev = eig(Xinf); - if (any(rev <= -ptol)) - x_ha_err = 5; - endif - endif - endif - - # build hamiltonian Ha for Y_inf - xx = ([CC'; -B1*ddot1']/Rtilde) * [ddot1*B1' CC]; - Ha = [A' 0*A; -B1*B1' -A] - xx; - y_ha_err = 0; - # copied from from are(...)... - [d, Ha] = balance(Ha); - [u, s] = schur(Ha, "A"); - rev = real(eig(s)); - if (any(abs(rev) <= ptol)) - # eigenvalues near the imaginary axis - y_ha_err = 1; - elseif (sum(rev > 0) != sum(rev < 0)) - # unequal number of positive and negative eigenvalues - y_ha_err = 2; - else - # compute positive Riccati equation solution - u = d * u; - Yinf = u(n+1:2*n,1:n) / u(1:n,1:n); - if (!all(all(finite(Yinf)))) - y_ha_err = 3; - elseif (norm(Yinf-Yinf') >= 10*ptol) - # solution not symmetric - x_ha_err = 4; - else - # positive semidefinite? - rev = eig(Yinf); - if (any(rev <= -ptol)) - y_ha_err = 5; - endif - endif - endif + [Xinf,x_ha_err] = hinfsyn_ric(A,BB,C1,d1dot,R,ptol); + [Yinf,y_ha_err] = hinfsyn_ric(A',CC',B1',ddot1',Rtilde,ptol); # assume failure for this gamma passed = 0; + rerr=""; if (!x_ha_err && !y_ha_err) # test spectral radius condition rho = max(abs(eig(Xinf * Yinf))); if (rho < g*g) # spectral radius condition passed passed = 1; + else + rerr = sprintf("rho=%f",rho); endif endif - switch (x_ha_err) - case (0) - # iax nev ene sym pos - printf(" o o o o o "); - xerr = " "; - case (1) - printf(" # - - - - "); - xerr = "X im.eig."; - case (2) - printf(" o # - - - "); - xerr = "Hx not Ham."; - case (3) - printf(" o o # - - "); - xerr = "X inf.eig."; - case (4) - printf(" o o o # - "); - xerr = "X not symm."; - case (5) - printf(" o o o o # "); - xerr = "X not pos."; - otherwise - error(" *** Xinf fail: this should never happen!"); - endswitch - switch (y_ha_err) - case (0) - # iax nev ene sym pos rho - if (passed) - printf(" o o o o o o y all tests passed.\n"); - elseif (x_ha_err) - printf(" o o o o o - n %s\n", xerr); - else - printf(" o o o o o # n rho=%f.\n", rho); - endif - case (1) - printf(" # - - - - - n %s/Y im. eig.", xerr); - case (2) - printf(" o # - - - - n %s/Hy not Ham.", xerr); - case (3) - printf(" o o # - - - n %s/Y inf.eig.", xerr); - case (4) - printf(" o o o # - - n %s/Y not symm.", xerr); - case (5) - printf(" o o o o # - n %s/Y not pos.", xerr); - otherwise - error(" *** Yinf fail: this should never happen!"); - endswitch + if(x_ha_err >= 0 & x_ha_err <= 6) + printf("%s",nth(errmesg,x_ha_err+1)); + xerr = nth(errdesx,x_ha_err+1); + else + error(" *** Xinf fail: this should never happen!"); + endif + if(y_ha_err >= 0 & y_ha_err <= 6) + printf("%s",nth(errmesg,y_ha_err+1)); + yerr = nth(errdesy,y_ha_err+1); + else + error(" *** Yinf fail: this should never happen!"); + endif + + if(passed) printf(" y all tests passed.\n"); + else printf(" n %s/%s%s\n",xerr,yerr,rerr); endif if (passed && (de/g < gtol)) search_state = 3; @@ -312,25 +239,19 @@ search_state = 1; endif case (1) - if (!passed) - search_state = 2; + if (!passed) search_state = 2; else # lower bound must not pass but passed fprintf(" *** the lower bound of gamma (%f) passed.\n", g); iteration_finished = 3; endif case (2) - if (!passed) - glo = g; - else - ghi = g; - endif - de = ghi - glo; - case (3) - # done - iteration_finished = 1; - otherwise - error(" *** This should never happen!"); + # Normal case; must check that singular R, Rtilde wasn't the problem. + if ((!passed) & (x_ha_err != 6) & (y_ha_err != 6) ) glo = g; + else ghi = g; endif + de = ghi - glo; + case (3) iteration_finished = 1; # done + otherwise error(" *** This should never happen!"); endswitch endwhile @@ -366,5 +287,4 @@ endif - page_screen_output = old_page_val; endfunction
--- a/scripts/control/hinfsyn_chk.m +++ b/scripts/control/hinfsyn_chk.m @@ -35,7 +35,6 @@ # Do not attempt to use this at home; no argument checking performed. # A. S. Hodel August 1995 - # $Revision: 2.0.0.2 $ Pc = Pf = [];
new file mode 100644 --- /dev/null +++ b/scripts/control/hinfsyn_ric.m @@ -0,0 +1,72 @@ +# Copyright (C) 1996,1998 A. Scottedward Hodel +# +# This file is part of Octave. +# +# Octave 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 2, or (at your option) any +# later version. +# +# Octave 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 Octave; see the file COPYING. If not, write to the Free +# Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. + +function [Xinf,x_ha_err] = hinfsyn_ric(A,BB,C1,d1dot,R,ptol) +# +# forms +# xx = ([BB; -C1'*d1dot]/R) * [d1dot'*C1 BB']; +# Ha = [A 0*A; -C1'*C1 -A'] - xx; +# and solves associated Riccati equation +# returns error code +# x_ha_err: +# 0: successful +# 1: Xinf has imaginary eigenvalues +# 2: Hx not Hamiltonian +# 3: Xinf has inf. eigenvalues (numerical overflow) +# 4: Xinf not symmetric +# 5: Xinf not positive definite +# 6: R is singular + +x_ha_err = 0; # assume success +Xinf = []; # default return value +n = is_square(A); +nw = is_square(R); +if(rank(R) != nw) x_ha_err = 6; +else # build hamiltonian Ha for X_inf + xx = ([BB; -C1'*d1dot]/R) * [d1dot'*C1 BB']; + Ha = [A 0*A; -C1'*C1 -A'] - xx; + x_ha_err = 0; + [d, Ha] = balance(Ha); + [u, s] = schur(Ha, "A"); + rev = real(eig(s)); + + if (any(abs(rev) <= ptol)) # eigenvalues near the imaginary axis + x_ha_err = 1; + elseif (sum(rev > 0) != sum(rev < 0)) + # unequal number of positive and negative eigenvalues + x_ha_err = 2; + else + # compute positive Riccati equation solution + u = d * u; + Xinf = u(n+1:2*n,1:n) / u(1:n,1:n); + if (!all(all(finite(Xinf)))) + x_ha_err = 3; + elseif (norm(Xinf-Xinf') >= 10*ptol) + # solution not symmetric + x_ha_err = 4; + else + # positive semidefinite? + # force symmetry (faster, avoids some convergence problems) + Xinf = (Xinf + Xinf')/2; + rev = eig(Xinf); + if (any(rev <= -ptol)) + x_ha_err = 5; + endif + endif + endif +endif
--- a/scripts/control/impulse.m +++ b/scripts/control/impulse.m @@ -36,7 +36,6 @@ # Written by Kai P. Mueller October 2, 1997 # based on lsim.m of Scottedward Hodel # modified by -# $Revision: 2.0.0.2 $ if((nargin < 1) || (nargin > 4)) usage("[y, u] = impulse(sys[, inp, tstop, n])");
--- a/scripts/control/is_abcd.m +++ b/scripts/control/is_abcd.m @@ -29,7 +29,6 @@ # Written by Kai P. Mueller November 4, 1997 # based on is_controllable.m of Scottedward Hodel # modified by - # $Revision: 2.0.0.2 $ retval = 0; switch (nargin)
--- a/scripts/control/is_controllable.m +++ b/scripts/control/is_controllable.m @@ -37,7 +37,6 @@ # Written by A. S. Hodel (scotte@eng.auburn.edu) August, 1993. # Updated by A. S. Hodel (scotte@eng.auburn.edu) Aubust, 1995 to use krylovb # Updated by John Ingram (ingraje@eng.auburn.edu) July, 1996 for packed systems -# $Revision: 1.16 $ deftol = 1; # assume default tolerance if(nargin < 1 | nargin > 3)
--- a/scripts/control/is_detectable.m +++ b/scripts/control/is_detectable.m @@ -31,7 +31,6 @@ # Written by A. S. Hodel (scotte@eng.auburn.edu) August 1993. # Updated by John Ingram (ingraje@eng.auburn.edu) July 1996. -# $Revision: 2.0.0.2 $ if( nargin < 1) usage("[retval,U] = is_detectable(a , c {, tol})");
--- a/scripts/control/is_dgkf.m +++ b/scripts/control/is_dgkf.m @@ -70,7 +70,6 @@ # Written by A. S. Hodel # Updated by John Ingram July 1996 to accept structured systems - # $Revision: 2.0.0.2 $ # # Revised by Kai P Mueller April 1998 to solve the general H_infinity # problem using unitary transformations Q (on w and z) @@ -92,7 +91,6 @@ # # This transformation together with the algorithm in [1] solves # the general problem (see [2] for example). - # $Revision: 2.0.0.2 $ if (nargin < 3) | (nargin > 4) usage("[retval,dgkf_struct] = is_dgkf(Asys,nu,ny{,tol})");
--- a/scripts/control/is_digital.m +++ b/scripts/control/is_digital.m @@ -22,7 +22,6 @@ # exits with an error of sys is a mixed (continuous and discrete) system # a s hodel July 1996 -# $Revision: 2.0.0.2 $ # SYS_INTERNAL accesses members of system structure # checked for sampled data system (mixed)
--- a/scripts/control/is_observable.m +++ b/scripts/control/is_observable.m @@ -31,7 +31,6 @@ # Written by A. S. Hodel (scotte@eng.auburn.edu) August 1993. # Updated by John Ingram (ingraje@eng.auburn.edu) July 1996. -# $Revision: 1.16 $ if( nargin < 1) usage("[retval,U] = is_observable(a , c {, tol})");
--- a/scripts/control/is_sample.m +++ b/scripts/control/is_sample.m @@ -22,7 +22,6 @@ # (real,scalar, > 0) # A. S. Hodel July 1995 -# $Revision: 2.0.0.2 $ out = (is_scalar(Ts) && (Ts == abs(Ts)) && (Ts != 0) );
new file mode 100644 --- /dev/null +++ b/scripts/control/is_signal_list.m @@ -0,0 +1,33 @@ +# Copyright (C) 1996,1998 A. Scottedward Hodel +# +# This file is part of Octave. +# +# Octave 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 2, or (at your option) any +# later version. +# +# Octave 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 Octave; see the file COPYING. If not, write to the Free +# Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. + +function flg = is_signal_list(mylist) +# function flg = is_signal_list(mylist) +# returns true if mylist is a list of individual strings. +# + +flg = is_list(mylist); +if(flg) + for ii=1:length(mylist) + if(!(isstr(nth(mylist,ii)) & rows(nth(mylist,ii)) ==1) ) + flg = 0; + endif + endfor +endif + +endfunction
--- a/scripts/control/is_siso.m +++ b/scripts/control/is_siso.m @@ -21,8 +21,6 @@ # return nonzero if the system sys is single-input, single-output. # a s hodel July 1996, 1998 -# $Revision: 2.0.0.2 $ -# SYS_INTERNAL accesses members of system structure if(nargin != 1) usage("SISO = is_siso(sys)");
--- a/scripts/control/is_stabilizable.m +++ b/scripts/control/is_stabilizable.m @@ -38,7 +38,6 @@ # Written by A. S. Hodel (scotte@eng.auburn.edu) August, 1993. # Updated by A. S. Hodel (scotte@eng.auburn.edu) Aubust, 1995 to use krylovb # Updated by John Ingram (ingraje@eng.auburn.edu) July, 1996 to accept systems -# $Revision: 2.0.0.2 $ if(nargin < 1) usage("[retval,U] = is_stabilizable(a {, b ,tol})"); elseif(is_struct(a))
--- a/scripts/control/is_stable.m +++ b/scripts/control/is_stable.m @@ -33,7 +33,6 @@ # Written by A. S. Hodel (scotte@eng.auburn.edu) August, 1993. # Updated by John Ingram (ingraje@eng.auburn.edu) July, 1996 for systems # Updated to simpler form by a.s.hodel 1998 -# $Revision: 2.0.0.2 $ if( (nargin < 1) | (nargin > 3) ) usage("is_stable(a {,tol,disc})"); elseif(is_struct(a))
--- a/scripts/control/jet707.m +++ b/scripts/control/jet707.m @@ -28,7 +28,6 @@ # Written by Kai P. Mueller September 28, 1997 # Updates - # $Revision: 2.0.0.2 $ if (nargin != 0) usage("outsys = jet707()")
--- a/scripts/control/lqe.m +++ b/scripts/control/lqe.m @@ -40,7 +40,6 @@ # e = closed loop poles of (A - K C) # Written by A. S. Hodel (scotte@eng.auburn.edu) August, 1993. -# $Revision: 1.16 $ if ( (nargin != 5) && (nargin != 6)) error ("lqe: invalid number of arguments");
--- a/scripts/control/lqg.m +++ b/scripts/control/lqg.m @@ -48,7 +48,6 @@ # Written by A. S. Hodel August 1995; revised for new system format # August 1996 -# $Revision: 2.0.0.2 $ sav_val = implicit_str_to_num_ok; implicit_str_to_num_ok = 1;
--- a/scripts/control/lqr.m +++ b/scripts/control/lqr.m @@ -47,7 +47,6 @@ # reference: Anderson and Moore, OPTIMAL CONTROL: LINEAR QUADRATIC METHODS, # Prentice-Hall, 1990, pp. 56-58 -# $Revision: 1.15 $ # Written by A. S. Hodel (scotte@eng.auburn.edu) August 1993.
--- a/scripts/control/lsim.m +++ b/scripts/control/lsim.m @@ -35,7 +35,6 @@ # Written by David Clem, A. S. Hodel July 1995 # modified by John Ingram for system format August 1996 -# $Revision: 2.0.0.2 $ if((nargin < 3)||(nargin > 4))
--- a/scripts/control/ltifr.m +++ b/scripts/control/ltifr.m @@ -32,7 +32,6 @@ # R. B. Tenison, D. Clem, A. S. Hodel, July 1995 # updated by John Ingram August 1996 for system format - # $Revision: 2.0.0.2 $ if ((nargin < 2) || (nargin > 3)) error("incorrect number of input arguments");
--- a/scripts/control/mb.m +++ b/scripts/control/mb.m @@ -1,8 +1,24 @@ +# Copyright (C) 1996,1998 A. Scottedward Hodel +# +# This file is part of Octave. +# +# Octave 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 2, or (at your option) any +# later version. +# +# Octave 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 Octave; see the file COPYING. If not, write to the Free +# Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. + # I think that this m-file can be deleted # a.s.hodel@eng.auburn.edu - 4 Dec. 1998 -# $Revision: 2.0.0.2 $ - Ap = [0 1;1960 0]; Bp = [0;-6261]; Cp = [1 0];
--- a/scripts/control/minfo.m +++ b/scripts/control/minfo.m @@ -33,7 +33,6 @@ # Written by R. Bruce Tenison July 29, 1994 # Modified by David Clem November 13, 1994 # Modified by A. S. Hodel July 1995 - # $Revision: 2.0.0.2 $ warning("minfo: obsolete. Use sys2ss, sys2tf, or sys2zp.");
--- a/scripts/control/moddemo.m +++ b/scripts/control/moddemo.m @@ -20,7 +20,6 @@ # Octave Controls toolbox demo: Model Manipulations demo # Written by David Clem August 15, 1994 -# $Revision: 2.0.0.2 $ # a s hodel: updated to reflect updated output order in ss2zp while (1)
--- a/scripts/control/nichols.m +++ b/scripts/control/nichols.m @@ -92,7 +92,6 @@ tistr = "(jw)"; endif xlabel("Phase (deg)"); - ylabel("Gain in dB"); if(is_siso(sys)) title(["Nichols plot of |[Y/U]",tistr,"|, u=", ... sysgetsignals(sys,"in",1,1), ", y=",sysgetsignals(sys,"out",1,1)]); @@ -101,7 +100,14 @@ printf("MIMO plot from\n%s\nto\n%s\n",outlist(inname," "), ... outlist(outname," ")); endif - md = 20*log10(mag); + if(max(mag) > 0) + ylabel("Gain in dB"); + md = 20*log10(mag); + else + ylabel("Gain |Y/U|") + md = mag; + endif + axvec = axis2dlim([vec(phase),vec(md)]); axis(axvec); plot(phase,md);
--- a/scripts/control/obsv.m +++ b/scripts/control/obsv.m @@ -39,7 +39,6 @@ # Written by Kai P. Mueller November 4, 1997 # modified by - # $Revision: 2.0.0.2 $ if (nargin == 2) a = sys;
--- a/scripts/control/ord2.m +++ b/scripts/control/ord2.m @@ -38,7 +38,6 @@ # Written by Kai P. Mueller September 28, 1997 # Updates - # $Revision: 2.0.0.2 $ if(nargin != 2 & nargin != 3) usage("outsys = ord2(nfreq, damp[, gain])")
--- a/scripts/control/outlist.m +++ b/scripts/control/outlist.m @@ -34,7 +34,6 @@ # prints the list to the screen, numbering each string in order. # A. S. Hodel Dec. 1995, 1998 -# $Revision: 2.0.0.2 $ #save for restore later save_empty = empty_list_elements_ok;
--- a/scripts/control/packedform.m +++ b/scripts/control/packedform.m @@ -16,7 +16,6 @@ # along with Octave; see the file COPYING. If not, write to the Free # Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. -# $Revision# save_var = page_screen_output; page_screen_output = 1; disp("Description of system data structure:")
--- a/scripts/control/packsys.m +++ b/scripts/control/packsys.m @@ -32,7 +32,6 @@ # Written by R. Bruce Tenison July 29, 1994 # Modified by David Clem November 13, 1994 # Modified by A. S. Hodel April 1995 - # $Revision: 2.0.0.2 $ warning("packsys is obsolete! Use ss2sys instead.");
--- a/scripts/control/parallel.m +++ b/scripts/control/parallel.m @@ -33,7 +33,6 @@ # Written by David Clem August 15, 1994 # completely rewritten Oct 1996 a s hodel # SYS_INTERNAL accesses members of system structure -# $Revision: 2.0.0.2 $ if(nargin != 2) usage("sysp = parallel(Asys,Bsys)"); @@ -48,10 +47,18 @@ if(mA != mB) error(["Asys has ",num2str(mA)," inputs, Bsys has ",num2str(mB)," inputs"]); endif + + # save signal names + Ain = sysgetsignals(Asys,"in"); + + # change signal names to avoid warning messages from sysgroup + Asys = syssetsignals(Asys,"in",sysdefioname(length(Ain),"Ain_u")); + Bsys = syssetsignals(Bsys,"in",sysdefioname(length(Ain),"Bin_u")); + sysp = sysgroup(Asys,Bsys); sysD = ss2sys([],[],[],[eye(mA);eye(mA)]); sysp = sysmult(sysp,sysD); - sysp = syssetsignals(sysp,"in",sysgetsignals(Asys,"in")); + sysp = syssetsignals(sysp,"in",Ain); endfunction
new file mode 100644 --- /dev/null +++ b/scripts/control/pinv.m @@ -0,0 +1,54 @@ +# Copyright (C) 1994 John W. Eaton +# +# This file is part of Octave. +# +# Octave 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 2, or (at your option) any +# later version. +# +# Octave 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 Octave; see the file COPYING. If not, write to the Free +# Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. + +function retval = pinv (X, tol) + +# usage: pinv (X, tol) +# +# Returns the pseudoinverse of X; singular values less than tol are +# ignored. +# +# If the second argument is omitted, it is assumed that +# +# tol = max (size (X)) * sigma_max (X) * eps, +# +# where sigma_max(X) is the maximal singular value of X. + +# Written by Kurt Hornik (hornik@neuro.tuwien.ac.at) March 1993. +# Dept of Probability Theory and Statistics TU Wien, Austria. + + if (nargin < 1 || nargin > 2) + error ("usage: pinv (X [, tol])"); + endif + + [U, S, V] = svd(X); + s = diag(S); + + if (nargin == 1) + tol = max (size (X)) * s (1) * eps; + endif + + r = sum (s > tol); + if (r == 0) + retval = zeros (X'); + else + D = diag (ones (r, 1) ./ s (1:r)); + retval = V (:, 1:r) * D * U (:, 1:r)'; + endif + +endfunction
--- a/scripts/control/place.m +++ b/scripts/control/place.m @@ -40,7 +40,6 @@ # # code adaped by A.S.Hodel (a.s.hodel@eng.auburn.edu) for use in controls # toolbox -# $Revision: 2.0.0.2 $ sav_val = empty_list_elements_ok; empty_list_elements_ok = 1;
--- a/scripts/control/polyout.m +++ b/scripts/control/polyout.m @@ -30,7 +30,6 @@ # Written by A. Scottedward Hodel (scotte@eng.auburn.edu) May 1995) # Nov 1998: Correctly handles complex coefficients -# $Revision: 2.0.0.2 $ if (nargin < 1 ) || (nargin > 2) || (nargout < 0 ) || (nargout > 1) usage("[y = ] polyout(c,[x])");
--- a/scripts/control/prompt.m +++ b/scripts/control/prompt.m @@ -23,7 +23,6 @@ # Written by David Clem August 15, 1994 # Modified A. S. Hodel June 1995 -# $Revision: 2.0.0.2 $ if(nargin > 1) usage("prompt([str])");
--- a/scripts/control/rldemo.m +++ b/scripts/control/rldemo.m @@ -20,7 +20,6 @@ # Octave Controls toolbox demo: Root Locus demo # Written by David Clem August 15, 1994 # Updated by John Ingram December 1996 -# $Revision: 2.0.0.2 $ while (1) clc
--- a/scripts/control/rlocus.m +++ b/scripts/control/rlocus.m @@ -42,7 +42,6 @@ # Written by Clem and Tenison # Updated by Kristi McGowan July 1996 for intelligent gain selection # Updated by John Ingram July 1996 for systems - # $Revision: 2.0.0.2 $ if (nargin < 1) | (nargin > 4) usage("rlocus(sys[,inc,mink,maxk])");
--- a/scripts/control/rotg.m +++ b/scripts/control/rotg.m @@ -22,7 +22,5 @@ # # NOTE: Use [c,s] = givens(a,b) instead. - # $Revision: 2.0.0.2 $ - [c,s] = givens(a,b); endfunction
--- a/scripts/control/run_cmd.m +++ b/scripts/control/run_cmd.m @@ -1,8 +1,24 @@ +# Copyright (C) 1996,1998 A. Scottedward Hodel +# +# This file is part of Octave. +# +# Octave 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 2, or (at your option) any +# later version. +# +# Octave 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 Octave; see the file COPYING. If not, write to the Free +# Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. + # run_cmd: short script used in demos # prints string cmd to the screen, then executes after a pause -# $Revision: 2.0.0.2 $ - disp(["Command: ",cmd]) puts("Press a key to execute command"); fflush(stdout);
--- a/scripts/control/series.m +++ b/scripts/control/series.m @@ -38,7 +38,6 @@ # Written by David Clem August 15, 1994 # If two arguments input, take care of mu system case -# $Revision: 2.0.0.2 $ warning("series is superseded by sysmult; use sysmult instead.")
--- a/scripts/control/sortcom.m +++ b/scripts/control/sortcom.m @@ -32,7 +32,6 @@ # idx: permutation vector: yy = xx(idx) # Written by A. S. Hodel June 1995 -# $Revision: 2.0.0.2 $ if( nargin < 1 | nargin > 2 ) usage("yy = sortcom(xx[,opt]"); @@ -46,34 +45,35 @@ error("sortcom: second argument must be a string"); endif endif - - if(strcmp(opt,"re")) datavec = real(xx); - elseif(strcmp(opt,"im")) datavec = imag(xx); - elseif(strcmp(opt,"mag")) datavec = abs(xx); + + if(isempty(xx)) + yy = idx = []; else - error(["sortcom: illegal option = ", opt]) - endif - - [datavec,idx] = sort(datavec); - yy= xx(idx); + if(strcmp(opt,"re")) datavec = real(xx); + elseif(strcmp(opt,"im")) datavec = imag(xx); + elseif(strcmp(opt,"mag")) datavec = abs(xx); + else error(["sortcom: illegal option = ", opt]) + endif - if(strcmp(opt,"re") | strcmp(opt,"mag")) - # sort so that complex conjugate pairs appear together + [datavec,idx] = sort(datavec); + yy= xx(idx); - ddiff = diff(datavec); - zidx = find(ddiff == 0); - - # sort common datavec values - if(!isempty(zidx)) - for iv=create_set(datavec(zidx)) - vidx = find(datavec == iv); - [vals,imidx] = sort(imag(yy(vidx))); - yy(vidx) = yy(vidx(imidx)); - idx(vidx) = idx(vidx(imidx)); - endfor + if(strcmp(opt,"re") | strcmp(opt,"mag")) + # sort so that complex conjugate pairs appear together + + ddiff = diff(datavec); + zidx = find(ddiff == 0); + + # sort common datavec values + if(!isempty(zidx)) + for iv=create_set(datavec(zidx)) + vidx = find(datavec == iv); + [vals,imidx] = sort(imag(yy(vidx))); + yy(vidx) = yy(vidx(imidx)); + idx(vidx) = idx(vidx(imidx)); + endfor + endif endif - - endif - + endif endfunction - +
--- a/scripts/control/ss2sys.m +++ b/scripts/control/ss2sys.m @@ -72,7 +72,6 @@ # # Written by John Ingram (ingraje@eng.auburn.edu) July 20, 1996 - # $Revision: 1.6 $ save_val = implicit_str_to_num_ok; # save for later implicit_str_to_num_ok = 1;
--- a/scripts/control/ss2tf.m +++ b/scripts/control/ss2tf.m @@ -34,7 +34,6 @@ # Written by R. Bruce Tenison (June 24, 1994) btenison@eng.auburn.edu # a s hodel: modified to allow for pure gain blocks Aug 1996 -# $Revision: 2.0.0.2 $ # Check args [n,m,p] = abcddim(a,b,c,d);
--- a/scripts/control/ss2zp.m +++ b/scripts/control/ss2zp.m @@ -27,7 +27,6 @@ # Written by David Clem August 15, 1994 # Hodel: changed order of output arguments to zer, pol, k. July 1996 # a s hodel: added argument checking, allow for pure gain blocks aug 1996 -# $Revision: 2.0.0.2 $ if(nargin != 4) usage("[zer,pol,k] = ss2zp(a,b,c,d)");
--- a/scripts/control/starp.m +++ b/scripts/control/starp.m @@ -46,7 +46,6 @@ # ny and/or nu may be negative (= negative feedback) # Written by Kai Mueller May 1998 -# $Revision: 2.0.0.2 $ if((nargin != 2) && (nargin != 4)) usage("[sys] = starp(P, K, ny, nu)");
--- a/scripts/control/step.m +++ b/scripts/control/step.m @@ -36,7 +36,6 @@ # Written by Kai P. Mueller September 30, 1997 # based on lsim.m of Scottedward Hodel # modified by -# $Revision: 2.0.0.2 $ if((nargin < 1) || (nargin > 4)) usage("[y, u] = step(sys[, inp, tstop, n])");
--- a/scripts/control/stepimp.m +++ b/scripts/control/stepimp.m @@ -31,7 +31,6 @@ # Written by Kai P. Mueller October 2, 1997 # based on lsim.m of Scottedward Hodel -# $Revision: 2.0.0.2 $ if (sitype == 1) IMPULSE = 0; elseif (sitype == 2) IMPULSE = 1;
--- a/scripts/control/swap.m +++ b/scripts/control/swap.m @@ -22,7 +22,6 @@ # A. S. Hodel July 24 1992 # Conversion to Octave R. Bruce Tenison July 4, 1994 - # $Revision: 2.0.0.2 $ a1 = b; b1 = a;
--- a/scripts/control/swapcols.m +++ b/scripts/control/swapcols.m @@ -22,7 +22,6 @@ # A. S. Hodel July 23, 1992 # Conversion to Octave R. Bruce Tenison July 4, 1994 - # $Revision: 2.0.0.2 $ m = length(A(1,:)); idx = m:-1:1;
--- a/scripts/control/swaprows.m +++ b/scripts/control/swaprows.m @@ -22,7 +22,6 @@ # A. S. Hodel July 23, 1992 # Conversion to Octave R. Bruce Tenison July 4, 1994 - # $Revision: 2.0.0.2 $ m = rows(A); idx = m:-1:1;
--- a/scripts/control/sys2fir.m +++ b/scripts/control/sys2fir.m @@ -20,7 +20,6 @@ # function [c,tsam,inname,outname] = sys2fir(sys) # extract fir system from system data structure -# $Revision: 2.0.0.2 $ # a s hodel July 1996 # let sys2tf do most of the work
--- a/scripts/control/sys2ss.m +++ b/scripts/control/sys2ss.m @@ -41,7 +41,6 @@ # Written by David Clem August 19, 1994 # Updates by John Ingram July 14, 1996 - # $Revision: 2.0.0.2 $ if(nargin != 1) usage("[a,b,c,d,tsam,n,nz,stname,inname,outname,yd] = sys2ss(sys)")
--- a/scripts/control/sys2tf.m +++ b/scripts/control/sys2tf.m @@ -31,7 +31,6 @@ # Written by R. Bruce Tenison (June 24, 1994) btenison@eng.auburn.edu # modified to make sys2tf by A. S. Hodel Aug 1995 # modified again for updated system format by John Ingram July 1996 -# $Revision: 2.0.0.2 $ if(nargin != 1) usage("[num,den,tsam,inname,outname] = sys2tf(Asys)");
--- a/scripts/control/sys2zp.m +++ b/scripts/control/sys2zp.m @@ -29,7 +29,6 @@ # inname, outname: input/output signal names (strings) # Created by John Ingram July 15 1996 -# $Revision: 2.0.0.2 $ if(nargin != 1) usage("[zer,pol,k,tsam,inname,outname] = sys2zp(sys)");
--- a/scripts/control/sysadd.m +++ b/scripts/control/sysadd.m @@ -36,7 +36,6 @@ # -------- # Written by John Ingram July 1996 -# $Revision: 2.0.0.2 $ save_val = implicit_str_to_num_ok; # save for later implicit_str_to_num_ok = 1; @@ -82,15 +81,17 @@ Gsys = sysupdate(Gsys,"ss"); Hsys = sysupdate(Hsys,"ss"); + # change signal names to avoid warning messages from sysgroup + Gsys = syssetsignals(Gsys,"in",sysdefioname(length(Gin),"Gin_u")); + Gsys = syssetsignals(Gsys,"out",sysdefioname(length(Gout),"Gout_u")); + Hsys = syssetsignals(Hsys,"in",sysdefioname(length(Hin),"Hin_u")); + Hsys = syssetsignals(Hsys,"out",sysdefioname(length(Hout),"Hout_u")); + sys = sysgroup(Gsys,Hsys); eyin = eye(mg); eyout = eye(pg); - - inname = Gin; - outname = Gout; - - sys = sysscale(sys,[eyout eyout],[eyin;eyin],outname,inname); + sys = sysscale(sys,[eyout eyout],[eyin;eyin],Gout,Gin); endfunction
--- a/scripts/control/sysappend.m +++ b/scripts/control/sysappend.m @@ -48,7 +48,6 @@ # sys = discrete: yd = ones(1,rows(c)) # written by John Ingram August 1996 - # $Revision: 2.0.0.2 $ sav_implicit_str_to_num_ok = implicit_str_to_num_ok; # save for later sav_empty_list_elements_ok = empty_list_elements_ok;
--- a/scripts/control/syschnames.m +++ b/scripts/control/syschnames.m @@ -40,7 +40,6 @@ # (or yd values, where appropriate) # Written by John Ingram August 1996; updated by A. S. Hodel 1998 -# $Revision: 2.0.0.2 $ retsys = syssetsignals(sys,opt,names,list);
--- a/scripts/control/syschtsam.m +++ b/scripts/control/syschtsam.m @@ -23,7 +23,6 @@ # This function changes the sampling time (tsam) of the system. # Written by John Ingram August 1996 -# $Revision: 2.0.0.2 $ if (nargin != 2) usage("retsys = syschtsam(sys,tsam)");
--- a/scripts/control/sysconnect.m +++ b/scripts/control/sysconnect.m @@ -50,7 +50,6 @@ # A. S. Hodel August 1995 # modified by John Ingram July 1996 -# $Revision: 2.0.0.2 $ save_val = implicit_str_to_num_ok; # save for later implicit_str_to_num_ok = 1;
--- a/scripts/control/syscont.m +++ b/scripts/control/syscont.m @@ -29,7 +29,6 @@ # returns csys empty if no continuous/continous path exists # Written by John Ingram August 1996 -# $Revision: 2.0.0.2 $ save_val = implicit_str_to_num_ok; # save for later save_empty = empty_list_elements_ok;
--- a/scripts/control/sysdefstname.m +++ b/scripts/control/sysdefstname.m @@ -21,8 +21,6 @@ # return default state names given n, nz # used internally, minimal argument checking -# $Revision: 2.0.0.2 $ - stname = list(); if(n > 0) for ii = 1:n
--- a/scripts/control/sysdisc.m +++ b/scripts/control/sysdisc.m @@ -26,8 +26,6 @@ # outputs, respectively. # -# $Revision: 2.0.0.2 $ - save_val = implicit_str_to_num_ok; # save for later save_empty = empty_list_elements_ok; empty_list_elements_ok = implicit_str_to_num_ok = 1;
--- a/scripts/control/sysdup.m +++ b/scripts/control/sysdup.m @@ -43,7 +43,6 @@ # A. S. Hodel August 1995 # modified by John Ingram July 1996 -# $Revision: 2.0.0.2 $ save_val = implicit_str_to_num_ok; # save for later implicit_str_to_num_ok = 1;
--- a/scripts/control/sysgettsam.m +++ b/scripts/control/sysgettsam.m @@ -1,9 +1,25 @@ +# Copyright (C) 1996,1998 A. Scottedward Hodel +# +# This file is part of Octave. +# +# Octave 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 2, or (at your option) any +# later version. +# +# Octave 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 Octave; see the file COPYING. If not, write to the Free +# Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. + function T = sysgettsam(sys) # T = sysgettsam(sys) # return the sampling time of the system -# $Revision: 2.0.0.2 $ - if(!is_struct(sys)) usage("T = sysgettsam(sys)"); endif
--- a/scripts/control/sysgroup.m +++ b/scripts/control/sysgroup.m @@ -40,7 +40,6 @@ # A. S. Hodel August 1995 # modified by John Ingram July 1996 -# $Revision: 2.0.0.2 $ save_val = implicit_str_to_num_ok; # save for later implicit_str_to_num_ok = 1;
--- a/scripts/control/sysgroupn.m +++ b/scripts/control/sysgroupn.m @@ -28,22 +28,21 @@ # # used internally in sysgroup -# $Revision: 2.0.0.2 $ - # check for duplicate names l = length(names); ii = 1; - while(ii < l-1) + while(ii <= l-1) st1 = nth(names,ii); jj = ii+1; - while ( jj < l) + while ( jj <= l) st2 = nth(names,jj); if(strcmp(st1,st2)) suffix = ["_",num2str(jj)]; - warning("sysgroup: %s name(%d) = %s name(%d); appending suffix %s to %d", ... - kind,ii,kind,jj,suffix,jj); + warning("sysgroup: %s name(%d) = %s name(%d) = %s", ... + kind,ii,kind,jj,st1); strval = sprintf("%s%s",st2,suffix); names(jj) = strval; + warning("sysgroup: changed %s name %d to %s",kind,jj,strval); # restart the check (just to be sure there's no further duplications) ii = 0; jj = l; endif
new file mode 100644 --- /dev/null +++ b/scripts/control/sysidx.m @@ -0,0 +1,50 @@ +function idxvec = sysidx(sys,sigtype,signamelist) +# idxvec = sysidx(sys,sigtype,signamelist) +# return indices of signals with specified signal names +# inputs: +# sys: OCST system data structure +# sigtype: signal type to be selected: "in", "out", "st" +# signamelist: list of desired signal names +# outputs: +# idxvec: vector of signal indices (appropriate for use with sysprune) + +if(nargin != 3) usage("idxvec = sysidx(sys,sigtype,signamelist)"); +elseif(!is_struct(sys)) error("sys must be a system data structure"); +elseif(!isstr(sigtype)) error("sigtype must be a string"); +elseif(rows(sigtype) != 1) + error("sigtype (%d x %d) must be a single string", ... + rows(sigtype),columns(sigtype)); +elseif(!is_signal_list(signamelist)) + error("signamelist must be a list of strings"); +endif + +sigtype_list = list("input","output","state"); +sigtnum = 0; +for idx = 1:length(sigtype_list) + thistype = nth(sigtype_list,idx); + if(strcmp(sigtype, thistype(1:length(sigtype)) )) sigtnum = idx; endif +endfor +if(sigtnum == 0) error("Illegal sigtype value = %s\n",sigtype); endif + +syssiglist = sysgetsignals(sys,sigtype); + +for idx = 1:length(signamelist) + signame = nth(signamelist,idx); + idxvec(idx) = 0; + for jdx = 1:sysdimensions(sys,sigtype) + if(strcmp(signame,sysgetsignals(sys,sigtype,jdx,1))) + if(idxvec(idx) != 0) + warning("Duplicate system input %s (%d,%d)\n", ... + sysgetsignals(sys,sigtype,jdx,1),jdx,idxvec(idx)); + else + idxvec(idx) = jdx; + endif + endif + endfor + if(idxvec(idx) == 0) + error("Did not find %s %s",nth(sigtype_list,sigtnum),signame); + endif +endfor + + +endfunction
--- a/scripts/control/sysmult.m +++ b/scripts/control/sysmult.m @@ -36,7 +36,6 @@ # does not recognize discrete inputs) # Written by John Ingram July 1996 -# $Revision: 2.0.0.2 $ save_val = implicit_str_to_num_ok; # save for later implicit_str_to_num_ok = 1; @@ -79,6 +78,10 @@ endif endif + # change signal names to avoid spurious warnings from sysgroup + Asys = syssetsignals(Asys,"in",sysdefioname(Am,"A_sysmult_tmp_name")); + Bsys = syssetsignals(Bsys,"out",sysdefioname(Bp,"B_sysmult_tmp_name")); + sys = sysgroup(Asys,Bsys); # connect outputs of B to inputs of A
--- a/scripts/control/sysout.m +++ b/scripts/control/sysout.m @@ -28,7 +28,6 @@ # "all": all of the above # Written by A S Hodel: 1995-1996 -# $Revision: 2.0.0.2 $ # save for restoring at end of routine save_val = implicit_str_to_num_ok;
--- a/scripts/control/sysprune.m +++ b/scripts/control/sysprune.m @@ -1,4 +1,4 @@ -# Copyright (C) 1996 A. Scottedward Hodel +# Copyright (C) 1996,1998 A. Scottedward Hodel # # This file is part of Octave. # @@ -16,21 +16,24 @@ # along with Octave; see the file COPYING. If not, write to the Free # Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. -function sys = sysprune(sys,output_list,input_list) -# function retsys = sysprune(Asys,output_list,input_list) -# Extract specified inputs/outputs from a system +function sys = sysprune(sys,output_idx,input_idx,state_idx) +# function retsys = sysprune(Asys,output_idx,input_idx{,state_idx}) +# Extract/permute specified inputs, outputs, and/or states of a system. # # inputs: # Asys: system data structure -# output_list,input_list: list of connections indices; the new -# system has outputs y(output_list(ii)) and inputs u(input_list(ii)). +# output_idx,input_idx: list of connections indices; the new +# system has outputs y(output_idx(ii)) and inputs u(input_idx(ii)). # May select as [] (empty matrix) to specify all outputs/inputs. +# state_idx: optional argument; list of indices of states to keep +# in the reduced model. May omit this argument or pass [] (empty +# matrix) to keep all states in the returned model # # output: retsys: resulting system: # ____________________ # | | # u1 ----->| |----> y1 -# (input_list) | Asys | (output_list) +# (input_idx) | Asys | (output_idx) # | | # u2 (deleted) |---->| |----| y2 (deleted) # | | @@ -38,58 +41,109 @@ # A. S. Hodel August 1995 # Updated by John Ingram 7-15-96 -# $Revision: 2.0.0.2 $ - if( nargin != 3 ) - usage("retsys = sysprune(sys,output_list,input_list)"); + if( nargin < 3 | nargin > 4 ) + usage("retsys = sysprune(sys,output_idx,input_idx{,state_idx})"); + elseif(nargin < 4) + state_idx = []; endif # default: no action [nn,nz,mm,pp] = sysdimensions(sys); - if(isempty(output_list)) - outputlist = 1:pp; - endif - if(isempty(input_list)) - input_list = 1:mm; - endif + if(isempty(output_idx)) output_idx = 1:pp; endif + if(isempty(input_idx)) input_idx = 1:mm; endif + if(isempty(state_idx)) state_idx = 1:(nn+nz); endif # check dimensions - if( !( - (is_vector(output_list) | isempty(output_list)) - & (is_vector(input_list) | isempty(input_list)) - )) - error("sysprune: output_list and input_list must be vectors"); - else - lo = length(output_list); - li = length(input_list); - - if(is_duplicate_entry(output_list) || is_duplicate_entry(input_list) ) - error("sysprune: duplicate entry in input and/or output list"); + if( !(is_vector(output_idx) | isempty(output_idx) ) ) + if(!is_matrix(output_idx)) + error("sysprune: bad argument passed for output_idx"); + else + error("sysprune: output_idx (%d x %d) must be a vector or empty", ... + rows(output_idx),columns(output_idx)); + endif + elseif(is_duplicate_entry(output_idx)) + error("sysprune: duplicate entries found in output_idx"); + endif + + if( !(is_vector(input_idx) | isempty(input_idx) ) ) + if(!is_matrix(input_idx)) + error("sysprune: bad argument passed for input_idx"); + else + error("sysprune: input_idx (%d x %d) must be a vector or empty", ... + rows(input_idx),columns(input_idx)); + endif + elseif(is_duplicate_entry(input_idx)) + error("sysprune: duplicate entries found in input_idx"); + endif + + if( !(is_vector(state_idx) | isempty(state_idx) ) ) + if(!is_matrix(state_idx)) + error("sysprune: bad argument passed for state_idx"); + else + error("sysprune: state_idx (%d x %d) must be a vector or empty", ... + rows(state_idx),columns(state_idx)); + endif + elseif(nn+nz > 0) + if(is_duplicate_entry(state_idx)) + error("sysprune: duplicate entries found in state_idx"); endif endif + + lo = length(output_idx); + li = length(input_idx); + lst = length(state_idx); - m = mm; - p = pp; - if( !is_struct(sys)) error("Asys must be a system data structure (see ss2sys, tf2sys, or zp2sys)") - elseif(p < lo) - error([num2str(lo)," output_list entries, system has only ", ... - num2str(p)," outputs"]); - elseif(m < li) - error([num2str(li)," input_list entries, system has only ", ... - num2str(m)," inputs"]); + elseif(pp < lo) + error([num2str(lo)," output_idx entries, system has only ", ... + num2str(pp)," outputs"]); + elseif(mm < li) + error([num2str(li)," input_idx entries, system has only ", ... + num2str(mm)," inputs"]); + elseif(nn+nz < lst) + error([num2str(lst)," state_idx entries, system has only ", ... + num2str(nn+nz)," states"]); endif [aa,bb,cc,dd,tsam,nn,nz,stnam,innam,outnam,yd] = sys2ss(sys); - bb = bb(:,input_list); - cc = cc(output_list,:); - dd = dd(output_list,input_list); - yd = yd(output_list); + + # check for legal state permutation + if(nn & nz) + c_idx = find(state_idx <= nn); + if(!isempty(c_idx)) max_c = max(c_idx); + else max_c = 0; endif + d_idx = find(state_idx > nn); + if(!isempty(d_idx)) min_d = min(d_idx); + else min_d = nn+nz; endif + if(max_c > min_d) + warning("sysprune: state_idx(%d)=%d (discrete) preceeds", ... + min_d,state_idx(min_d)); + warning(" state_idx(%d)=%d (continuous)",... + max_c,state_idx(max_c)); + warning("sysprune: sys has %d continuous states, %d discrete states", ... + nn,nz); + error("continuous/discrete state partition not preserved ; see ss2sys"); + endif + endif - # this part needs rewritten for the list structure of signal names - innam = innam(input_list); - outnam = outnam(output_list); - - sys = ss2sys(aa,bb,cc,dd,tsam,nn,nz,stnam,innam,outnam,find(yd)); + idx = input_idx; + odx = output_idx; + if(isempty(state_idx)) + idx = []; + odx = []; + endif + aa = aa(state_idx,state_idx); + bb = bb(state_idx,idx); + cc = cc(odx,state_idx); + dd = dd(output_idx,input_idx); + yd = yd(output_idx); + + innam = innam(input_idx); + outnam = outnam(output_idx); + stnam = stnam(state_idx); + nn1 = length(find(state_idx <= nn)); + nz1 = length(find(state_idx > nn)); + sys = ss2sys(aa,bb,cc,dd,tsam,nn1,nz1,stnam,innam,outnam,find(yd)); endfunction
--- a/scripts/control/sysreorder.m +++ b/scripts/control/sysreorder.m @@ -28,7 +28,6 @@ # to use this at home! # A. S. Hodel, Aug 1995 -# $Revision: 2.0.0.2 $ #disp('sysreorder: entry')
--- a/scripts/control/sysrepdemo.m +++ b/scripts/control/sysrepdemo.m @@ -23,8 +23,6 @@ # Written by A. S. Hodel June 1995 # Revised Aug 1995 for system data structure format -# $Revision: 2.0.0.2 $ - save_val = page_screen_output; page_screen_output = 1;
--- a/scripts/control/sysscale.m +++ b/scripts/control/sysscale.m @@ -39,7 +39,6 @@ # A. S. Hodel August 1995 # modified by John Ingram 7-15-96 -# $Revision: 2.0.0.2 $ if( (nargin < 3) || (nargin > 5) ) usage("retsys = sysscale(Asys,output_list,input_list{,inname,outname})");
new file mode 100644 --- /dev/null +++ b/scripts/control/syssetsignals.m @@ -0,0 +1,168 @@ +# Copyright (C) 1996,1998 A. Scottedward Hodel +# +# This file is part of Octave. +# +# Octave 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 2, or (at your option) any +# later version. +# +# Octave 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 Octave; see the file COPYING. If not, write to the Free +# Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. + +function retsys = syssetsignals(sys,opt,names,sig_idx) +# retsys = syssetsignals(sys,opt,names{,sig_idx}) +# change the names of selected inputs, outputs and states. +# inputs: +# sys: system data structure +# opt: []: change default name (output) +# "out": change selected output names +# "in": change selected input names +# "st": change selected state names +# "yd": change selected outputs from discrete to continuous or +# from continuous to discrete. +# names: opt = "out", "in", or "st": string or or list of strings +# opt = "yd": desired output continuous/discrete flag. +# names(ii) = 0: output ii is continuous +# names(ii) = 1: output ii is discrete +# list: vector of indices of outputs, yd, inputs, or +# states whose respective names should be changed. +# Default: replace entire signal list/vector value +# outputs: +# retsys=sys with appropriate signal names changed +# (or yd values, where appropriate) + +# Written by John Ingram August 1996 + + save_val = implicit_str_to_num_ok; # save for later + implicit_str_to_num_ok = 1; + + if (nargin < 3 | nargin > 4) + usage("retsys=syssetsignals(sys,opt,names{,sig_idx})"); + elseif (!is_struct(sys)) + error("sys must be a system data structure"); + elseif (isempty(opt)) + opt = "out"; + elseif( ! isstr(opt) ) + error("opt must be a string"); + elseif( ! (strcmp(opt,"out") + strcmp(opt,"yd") + ... + strcmp(opt,"in") + strcmp(opt,"st") ) ) + error("opt must be one of [], ""out"", ""yd"", ""in"", or ""st"""); + elseif(nargin == 4) + if(min(size(sig_idx)) > 1) + disp("syssetsignals: sig_idx=") + disp(sig_idx); + error("sig_idx must be a vector") + endif + endif + + sig_vals = sysgetsignals(sys,opt); + + # make sure it's in state space form if state names are given + if(strcmp(opt,"st")) sys = sysupdate(sys,"ss"); endif + + if(strcmp(opt,"yd") == 0) + # it's a signal name list we're changing + if(!is_list(names)) + names = list(names); + endif + if(!is_signal_list(names)) + if(isstr(nth(names,1))) + warning("syssetsignals(opt=%s): converting string matrix \"names\" to a list of strings",opt); + tmpstr = nth(names,1); + for ii=1:rows(tmpstr) + names(ii) = deblank(tmpstr(ii,:)); + endfor + else + names + error("parameter \"names\" must be a list of strings"); + endif + endif + nsigs = length(sig_vals); + + if(nargin == 3) + # replace all signal names + if(length(names) != nsigs) + error("opt=%s, sig_idx omitted: names(len=%d) should have %d entries ", ... + opt,length(names),nsigs); + endif + sig_idx = 1:nsigs; + elseif(length(names) != length(sig_idx)) + # replace specified signal names + error("opt=%s, sig_idx(len=%d), names(len=%d) mismatch",opt, ... + length(sig_idx), length(names)); + endif + + for ii=1:length(sig_idx) + jj = sig_idx(ii); + if(jj < 1 | jj > nsigs | jj != floor(jj+0.5)) + error("opt=%s, sig_idx(%d)=%d, %e: must be an integer between 1 and %d", ... + opt, ii, jj, jj, nsigs); + endif + sig_vals(jj) = nth(names,ii); + endfor + + else + nsigs = length(sig_vals); + if(!is_vector(names)) + error("syssetsignals: opt=yd, names(%dx%d) must be a vector", ... + rows(names), columns(names)); + endif + if(nargin == 3) + if(length(names) != nsigs) + error("opt=yd, sig_idx omitted: names(%d) should be length(%d)", ... + length(names), nsigs); + endif + sig_idx = 1:nsigs; + elseif(length(names) != length(sig_idx)) + error("opt=yd: length(names)=%d, length(sig_idx)=%d",length(names), ... + length(sig_idx) ); + endif + + badidx = find(names != 0 & names != 1); + if(! isempty(badidx) ) + for ii=1:length(badidx) + warning("syssetsignals: opt=yd: names(%d)=%e, must be 0 or 1", ... + badidx(ii), names(badidx(ii)) ); + endfor + error("opt=yd: illegal values in names"); + endif + + for ii=1:length(sig_idx) + jj = sig_idx(ii); + if(jj < 1 | jj > nsigs | jj != floor(jj)) + error("sig_idx(%d)=%d, %e: must be an integer between 1 and %d", ... + ii,jj, jj, nsigs); + endif + sig_vals(jj) = names(ii); + endfor + if(any(sig_vals == 1) & sysgettsam(sys) == 0) + warning("Setting system sampling time to 1"); + printf("syssetsignals: original system sampling time=0 but output(s)\n"); + disp(find(sig_vals==1)) + printf("are digital\n"); + sys = syschtsam(sys,1); + endif + + endif + + if(strcmp(opt,"st")) + sys.stname = sig_vals; + elseif(strcmp(opt,"in")) + sys.inname = sig_vals; + elseif(strcmp(opt,"out")) + sys.outname = sig_vals; + elseif(strcmp(opt,"yd")) + sys.yd = sig_vals; + endif + + retsys = sys; + implicit_str_to_num_ok = save_val; # restore value + +endfunction
--- a/scripts/control/syssub.m +++ b/scripts/control/syssub.m @@ -36,7 +36,6 @@ # -------- # Written by John Ingram July 1996 -# $Revision: 2.0.0.2 $ save_val = implicit_str_to_num_ok; # save for later implicit_str_to_num_ok = 1; @@ -82,13 +81,16 @@ Gsys = sysupdate(Gsys,"ss"); Hsys = sysupdate(Hsys,"ss"); + # change signal names to avoid warning messages from sysgroup + Gsys = syssetsignals(Gsys,"in",sysdefioname(length(Gin),"Gin_u")); + Gsys = syssetsignals(Gsys,"out",sysdefioname(length(Gout),"Gout_u")); + Hsys = syssetsignals(Hsys,"in",sysdefioname(length(Hin),"Hin_u")); + Hsys = syssetsignals(Hsys,"out",sysdefioname(length(Hout),"Hout_u")); + sys = sysgroup(Gsys,Hsys); eyin = eye(mg); eyout = eye(pg); - inname = Gin; - outname = Gout; - - sys = sysscale(sys,[eyout -eyout],[eyin;eyin],outname,inname); + sys = sysscale(sys,[eyout -eyout],[eyin;eyin],Gout,Gin); endfunction
--- a/scripts/control/sysupdate.m +++ b/scripts/control/sysupdate.m @@ -34,7 +34,6 @@ # see also: tf2sys, ss2sys, zp2sys, sysout, sys2ss, sys2tf, sys2zp # Written by John Ingram 7-9-96 -# $Revision: 2.0.0.2 $ # check for correct number of inputs if (nargin != 2)
--- a/scripts/control/tf2ss.m +++ b/scripts/control/tf2ss.m @@ -39,7 +39,6 @@ # Written by R. Bruce Tenison (June 22, 1994) btenison@eng.auburn.edu # mod A S Hodel July, Aug 1995 - # $Revision: 2.0.0.2 $ if(nargin != 2) error("tf2ss: wrong number of input arguments") elseif(isempty(num)) error("tf2ss: empty numerator");
--- a/scripts/control/tf2sys.m +++ b/scripts/control/tf2sys.m @@ -29,7 +29,6 @@ # Written by R. Bruce Tenison July 29, 1994 # Name changed to TF2SYS July 1995 # updated for new system data structure format July 1996 - # $Revision: 2.0.0.2 $ save_val = implicit_str_to_num_ok; implicit_str_to_num_ok = 1;
--- a/scripts/control/tf2sysl.m +++ b/scripts/control/tf2sysl.m @@ -22,7 +22,6 @@ # used internally in tf2sys # strip leading zero coefficients to get the true polynomial length -# $Revision: 2.0.0.2 $ while( (length(vec) > 1) & (vec(1) == 0) ) vec = vec(2:length(vec));
--- a/scripts/control/tf2zp.m +++ b/scripts/control/tf2zp.m @@ -23,7 +23,6 @@ # defined by num/den. K is a gain associated with the system zeros. # Written by A. S. Hodel, etc. -# $Revision: 2.0.0.2 $ if(nargin == 2) if(length(den) > 1) pol = roots(den);
--- a/scripts/control/tfout.m +++ b/scripts/control/tfout.m @@ -28,7 +28,6 @@ # filter, polyderiv, polyinteg, polyout # Written by A. Scottedward Hodel (scotte@eng.auburn.edu) June 1995) -# $Revision: 2.0.0.2 $ save_val = implicit_str_to_num_ok; save_empty = empty_list_elements_ok;
--- a/scripts/control/tzero.m +++ b/scripts/control/tzero.m @@ -39,7 +39,6 @@ # R. Bruce Tenison July 4, 1994 # A. S. Hodel Aug 1995: allow for MIMO and system data structures - # $Revision: 1.18 $ # get A,B,C,D and Asys variables, regardless of initial form if(nargin == 4)
--- a/scripts/control/tzero2.m +++ b/scripts/control/tzero2.m @@ -27,7 +27,6 @@ # Needs to incorporate mvzero algorithm to isolate finite zeros. # Written by A. S. Hodel (scotte@eng.auburn.edu) August 1993. -# $Revision: 2.0.0.2 $ if (nargin == 4) bal = "B";
--- a/scripts/control/ugain.m +++ b/scripts/control/ugain.m @@ -28,7 +28,6 @@ # Written by Kai P. Mueller April, 1998 # Updates - # $Revision: 2.0.0.2 $ if((nargin != 1) || (nargout > 1)) usage("outsys = ugain(n)")
--- a/scripts/control/unpacksys.m +++ b/scripts/control/unpacksys.m @@ -21,7 +21,6 @@ # Obsolete. Use sys2ss instead. # Written by David Clem August 19, 1994 - # $Revision: 2.0.0.2 $ warning("unpacksys obsolete; calling sys2ss"); [a,b,c,d] = sys2ss(syst);
--- a/scripts/control/wgt1o.m +++ b/scripts/control/wgt1o.m @@ -30,7 +30,6 @@ # fc = Corner frequency (in Hz, *not* in rad/sec) # Written by Kai P. Mueller September 30, 1997 -# $Revision: 2.0.0.2 $ if (nargin != 3) usage("wsys = wgt1o(vl, vh, fc)");
--- a/scripts/control/zgfmul.m +++ b/scripts/control/zgfmul.m @@ -30,7 +30,6 @@ # A. S. Hodel July 24 1992 # Conversion to Octave July 3, 1994 - # $Revision: 2.0.0.2 $ [n,m] = size(b); [p,m1] = size(c);
--- a/scripts/control/zgfslv.m +++ b/scripts/control/zgfslv.m @@ -22,7 +22,6 @@ # Written by A. Scotte Hodel # Converted to Octave by R Bruce Tenison, July 3, 1994 - # $Revision: 2.0.0.2 $ nmp = n+m+p; gam1 = (2*n)+m+p; gam2 = n+p; gam3 = n+m;
--- a/scripts/control/zginit.m +++ b/scripts/control/zginit.m @@ -29,7 +29,6 @@ # A. S. Hodel July 24 1992 # Conversion to Octave by R. Bruce Tenison, July 3, 1994 - # $Revision: 2.0.0.2 $ [nn,mm] = size(b); [pp,mm] = size(d);
--- a/scripts/control/zgpbal.m +++ b/scripts/control/zgpbal.m @@ -38,7 +38,6 @@ # A. S. Hodel July 24 1992 # Conversion to Octave by R. Bruce Tenison July 3, 1994 - # $Revision: 2.0.0.2 $ if( (nargin != 1) | (!is_struct(Asys))) usage("retsys = zgpbal(Asys)");
--- a/scripts/control/zgreduce.m +++ b/scripts/control/zgreduce.m @@ -23,7 +23,6 @@ # # used internally in tzero; minimal argument checking performed -#$Revision: 2.0.0.2 $ # SYS_INTERNAL accesses members of system data structure is_digital(Asys); # make sure it's pure digital/continuous
--- a/scripts/control/zgrownorm.m +++ b/scripts/control/zgrownorm.m @@ -22,7 +22,6 @@ # returns nonz = number of rows of mat whose two norm exceeds meps # zer = number of rows of mat whose two norm is less than meps -# $Revision: 2.0.0.2 $ rownorm = []; for ii=1:rows(mat)
--- a/scripts/control/zgscal.m +++ b/scripts/control/zgscal.m @@ -28,7 +28,7 @@ # A. S. Hodel July 24 1992 # Conversion to Octave R. Bruce Tenison July 3, 1994 - # $Revision: 2.0.0.2 $ + #************************************************************************** #initialize parameters: @@ -43,6 +43,41 @@ rk1 = z; k = 0; + # construct balancing least squares problem + F = eye(nmp); + for kk=1:nmp + F(1:nmp,kk) = zgfmul(a,b,c,d,F(:,kk)); + endfor + + [U,H,k1] = krylov(F,z,nmp,1e-12); + if(!is_square(H)) + if(columns(H) != k1) + error("zgscal(tzero): k1=%d, columns(H)=%d",k1,columns(H)); + elseif(rows(H) != k1+1) + error("zgscal: k1=%d, rows(H) = %d",k1,rows(H)); + elseif ( norm(H(k1+1,:)) > 1e-12*norm(H,"inf") ) + zgscal_last_row_of_H = H(k1+1,:) + error("zgscal: last row of H nonzero (norm(H)=%e)",norm(H,"inf")) + endif + H = H(1:k1,1:k1); + U = U(:,1:k1); + endif + + # tridiagonal H can still be rank deficient, so do permuted qr + # factorization + [qq,rr,pp] = qr(H); # H = qq*rr*pp' + nn = rank(rr); + qq = qq(:,1:nn); + rr = rr(1:nn,:); # rr may not be square, but "\" does least + xx = U*pp*(rr\qq'*(U'*z)); # squares solution, so this works + #xx1 = pinv(F)*z; + #zgscal_x_xx1_err = [xx,xx1,xx-xx1] + return; + + # the rest of this is left from the original zgscal; + # I've had some numerical problems with the GCG algorithm, + # so for now I'm solving it with the krylov routine. + #initialize residual error norm rnorm = norm(rk1,1);
--- a/scripts/control/zgsgiv.m +++ b/scripts/control/zgsgiv.m @@ -23,7 +23,6 @@ # A. S. Hodel July 29, 1992 # Convertion to Octave by R. Bruce Tenison July 3, 1994 - # $Revision: 2.0.0.2 $ t1 = c*a + s*b; t2 = -s*a + c*b;
--- a/scripts/control/zgshsr.m +++ b/scripts/control/zgshsr.m @@ -24,7 +24,6 @@ # A. S. Hodel July 24, 1992 # Conversion to Octave by R. Bruce Tenison July 3, 1994 - # $Revision: 2.0.0.2 $ if(!is_vector(y)) error(sprintf("y(%dx%d) must be a vector",rows(y),columns(y)));
--- a/scripts/control/zp2ss.m +++ b/scripts/control/zp2ss.m @@ -39,7 +39,6 @@ # k is a gain that is associated with the zero vector. # Written by David Clem August 15, 1994 -# $Revision: 2.0.0.2 $ sav_val = empty_list_elements_ok; empty_list_elements_ok = 1;
--- a/scripts/control/zp2ssg2.m +++ b/scripts/control/zp2ssg2.m @@ -23,7 +23,6 @@ # extract 2 values from rvals (if possible) and construct # a polynomial with those roots. -# $Revision: 2.0.0.2 $ # A. S. Hodel Aug 1996 # locate imaginary roots (if any)