Mercurial > hg > octave-lojdl
changeset 3437:3c609681b2bf
[project @ 2000-01-14 06:00:00 by jwe]
author | jwe |
---|---|
date | Fri, 14 Jan 2000 06:00:04 +0000 |
parents | 655add875c71 |
children | 2e06c3941943 |
files | scripts/control/base/__bodquist__.m scripts/control/base/__freqresp__.m scripts/control/base/__stepimp__.m scripts/control/system/__abcddims__.m scripts/control/system/__syschnamesl__.m scripts/control/system/__syscont_disc__.m scripts/control/system/__sysdefioname__.m scripts/control/system/__sysdefstname__.m scripts/control/system/__sysgroupn__.m scripts/control/system/__tf2sysl__.m scripts/control/system/__zp2ssg2__.m scripts/control/util/__outlist__.m scripts/control/util/__zgpbal__.m |
diffstat | 13 files changed, 1224 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
new file mode 100644 --- /dev/null +++ b/scripts/control/base/__bodquist__.m @@ -0,0 +1,156 @@ +## Copyright (C) 1996, 1998 Auburn University. All rights reserved. +## +## 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, 59 Temple Place, Suite 330, Boston, MA 02111 USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{f}, @var{w}] =} bodquist (@var{sys}, @var{w}, @var{out_idx}, @var{in_idx}) +## used internally by bode, nyquist; compute system frequency response. +## +## @strong{Inputs} +## @table @var +## @item sys +## input system structure +## @item w +## range of frequencies; empty if user wants default +## @item out_idx +## list of outputs; empty if user wants all +## @item in_idx +## list of inputs; empty if user wants all +## @item rname +## name of routine that called bodquist ("bode" or "nyquist") +## @end table +## @strong{Outputs} +## @table @var +## @item w +## list of frequencies +## @item f +## frequency response of sys; @math{f(ii) = f(omega(ii))} +## @end table +## @strong{Note} bodquist could easily be incorporated into a Nichols +## plot function; this is in a "to do" list. +## +## Both bode and nyquist share the same introduction, so the common parts are +## in bodquist. It contains the part that finds the number of arguments, +## determines whether or not the system is SISO, and computes the frequency +## response. Only the way the response is plotted is different between the +## two functions. +## @end deftypefn + +function [f, w] = bodquist (sys, w, outputs, inputs, rname) + + ## check number of input arguments given + if (nargin != 5) + usage("[f,w] = bodquist(sys,w,outputs,inputs,rname)"); + endif + + ## check each argument to see if it's in the correct form + if (!is_struct(sys)) + error("sys must be a system data structure"); + endif + + ## let freqresp determine w if it's not already given + USEW = freqchkw(w); + + ## get initial dimensions (revised below if sysprune is called) + [nn,nz,mm,pp ] = sysdimensions(sys); + + ## check for an output vector and to see whether it`s correct + if (!isempty(outputs)) + if (isempty(inputs)) + inputs = 1:mm; # use all inputs + warning([rname,": outputs specified but not inputs"]); + endif + sys = sysprune(sys,outputs,inputs); + [nn,nz,mm,pp ] = sysdimensions(sys); + endif + + ## for speed in computation, convert local copy of + ## SISO state space systems to zero-pole form + if( is_siso(sys) & strcmp( sysgettype(sys), "ss") ) + [zer,pol,k,tsam,inname,outname] = sys2zp(sys); + sys = zp2sys(zer,pol,k,tsam,inname,outname); + endif + + ## get system frequency response + [f,w] = freqresp(sys,USEW,w); + + phase = arg(f)*180.0/pi; + + if(!USEW) + ## smooth plots + pcnt = 5; # max number of refinement steps + dphase = 5; # desired max change in phase + dmag = 0.2; # desired max change in magnitude + while(pcnt) + pd = abs(diff(phase)); # phase variation + pdbig = vec(find(pd > dphase)); + + lp = length(f); lp1 = lp-1; # relative variation + fd = abs(diff(f)); + fm = max(abs([f(1:lp1); f(2:lp)])); + fdbig = vec(find(fd > fm/10)); + + bigpts = union(fdbig, pdbig); + + if(isempty(bigpts) ) + pcnt = 0; + else + pcnt = pcnt - 1; + wnew = []; + crossover_points = find ( phase(1:lp1).*phase(2:lp) < 0); + pd(crossover_points) = abs(359.99+dphase - pd(crossover_points)); + np_pts = max(3,ceil(pd/dphase)+2); # phase points + nm_pts = max(3,ceil(log(fd./fm)/log(dmag))+2); # magnitude points + npts = min(5,max(np_pts, nm_pts)); + + w1 = log10(w(1:lp1)); + w2 = log10(w(2:lp)); + for ii=bigpts + if(npts(ii)) + wtmp = logspace(w1(ii),w2(ii),npts(ii)); + wseg(ii,1:(npts(ii)-2)) = wtmp(2:(npts(ii)-1)); + endif + endfor + wnew = vec(wseg)'; # make a row vector + wnew = wnew(find(wnew != 0)); + wnew = sort(wnew); + wnew = create_set(wnew); + if(isempty(wnew)) # all small crossovers + pcnt = 0; + else + [fnew,wnew] = freqresp(sys,1,wnew); # get new freq resp points + w = [w,wnew]; # combine with old freq resp + f = [f,fnew]; + [w,idx] = sort(w); # sort into order + f = f(idx); + phase = arg(f)*180.0/pi; + endif + endif + endwhile + endif + + ## ensure unique frequency values + [w,idx] = sort(w); + f = f(idx); + + w_diff = diff(w); + w_dup = find(w_diff == 0); + w_idx = complement(w_dup,1:length(w)); + w = w(w_idx); + f = f(w_idx); + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/control/base/__freqresp__.m @@ -0,0 +1,133 @@ +## Copyright (C) 1996, 1998 Auburn University. All rights reserved. +## +## 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, 59 Temple Place, Suite 330, Boston, MA 02111 USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{out} =} freqresp (@var{sys}, @var{USEW}@{,@var{w}@}); +## Frequency response function - used internally by @code{bode}, @code{nyquist}. +## minimal argument checking; "do not attempt to do this at home" +## +## @strong{Inputs} +## @table @var +## @item sys +## system data structure +## @item USEW +## returned by @code{freqchkw} +## @item optional +## must be present if @var{USEW} is true (nonzero) +## @end table +## @strong{Outputs} +## @table @var +## @item @var{out} +## vector of finite @math{G(j*w)} entries (or @math{||G(j*w)||} for MIMO) +## @item w +## vector of corresponding frequencies +## @end table +## @end deftypefn + +## Author: R. Bruce Tenison <btenison@eng.auburn.edu> +## Created: July 11, 1994 + +function [ff, w] = freqresp (sys, USEW, w); + + ## SYS_INTERNAL accesses members of system data structure + + save_val = empty_list_elements_ok; + empty_list_elements_ok = 1; + + ## Check Args + if( (nargin < 2) || (nargin > 4) ) + usage ("[ff,w] = freqresp(sys,USEW{,w})"); + elseif( USEW & (nargin < 3) ) + error("USEW=1 but w was not passed."); + elseif( USEW & isempty(w)) + warning("USEW=1 but w is empty; setting USEW=0"); + USEW=0; + endif + + DIGITAL = is_digital(sys); + + ## compute default w if needed + if(!USEW) + if(is_siso(sys)) + sys = sysupdate(sys,"zp"); + [zer,pol] = sys2zp(sys); + else + zer = tzero(sys); + pol = eig(sys2ss(sys)); + endif + + ## get default frequency range + [wmin,wmax] = bode_bounds(zer,pol,DIGITAL,sysgettsam(sys)); + w = logspace(wmin,wmax,50); + else + w = reshape(w,1,length(w)); # make sure it's a row vector + endif + + ## now get complex values of s or z + if(DIGITAL) + jw = exp(i*w*sysgettsam(sys)); + else + jw = i*w; + endif + + [nn,nz,mm,pp] = sysdimensions(sys); + + ## now compute the frequency response - divide by zero yields a warning + if (strcmp(sysgettype(sys),"zp")) + ## zero-pole form (preferred) + [zer,pol,sysk] = sys2zp(sys); + ff = ones(size(jw)); + l1 = min(length(zer)*(1-isempty(zer)),length(pol)*(1-isempty(pol))); + for ii=1:l1 + ff = ff .* (jw - zer(ii)) ./ (jw - pol(ii)); + endfor + + ## require proper transfer function, so now just get poles. + for ii=(l1+1):length(pol) + ff = ff ./ (jw - pol(ii)); + endfor + ff = ff*sysk; + + elseif (strcmp(sysgettype(sys),"tf")) + ## transfer function form + [num,den] = sys2tf(sys); + ff = polyval(num,jw)./polyval(den,jw); + elseif (mm==pp) + ## The system is square; do state-space form bode plot + [sysa,sysb,sysc,sysd,tsam,sysn,sysnz] = sys2ss(sys); + n = sysn + sysnz; + for ii=1:length(jw); + ff(ii) = det(sysc*((jw(ii).*eye(n)-sysa)\sysb)+sysd); + endfor; + else + ## Must be state space... bode + [sysa,sysb,sysc,sysd,tsam,sysn,sysnz] = sys2ss(sys); + n = sysn + sysnz; + for ii=1:length(jw); + ff(ii) = norm(sysc*((jw(ii)*eye(n)-sysa)\sysb)+sysd); + endfor + + endif + + w = reshape(w,1,length(w)); + ff = reshape(ff,1,length(ff)); + + ## restore global variable + empty_list_elements_ok = save_val; +endfunction +
new file mode 100644 --- /dev/null +++ b/scripts/control/base/__stepimp__.m @@ -0,0 +1,279 @@ +## Copyright (C) 1996, 1998 Auburn University. All rights reserved. +## +## 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, 59 Temple Place, Suite 330, Boston, MA 02111 USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{y}, @var{t}] = } stepimp (@var{sitype}, @var{sys} [, @var{inp}, @var{tstop}, @var{n}]) +## Impulse or step response for a linear system. +## The system can be discrete or multivariable (or both). +## This m-file contains the "common code" of step and impulse. +## +## Produces a plot or the response data for system sys. +## +## Limited argument checking; "do not attempt to do this at home". +## Used internally in @code{impulse}, @code{step}. Use @code{step} +## or @code{impulse} instead. +## @end deftypefn +## @seealso{step and impulse} + +## Author: Kai P. Mueller <mueller@ifr.ing.tu-bs.de> +## Created: October 2, 1997 +## based on lsim.m of Scottedward Hodel + +function [y, t] = stepimp (sitype, sys, inp, tstop, n) + + if (sitype == 1) IMPULSE = 0; + elseif (sitype == 2) IMPULSE = 1; + else error("stepimp: invalid sitype argument.") + endif + sys = sysupdate(sys,"ss"); + + USE_DEF = 0; # default tstop and n if we have to give up + N_MIN = 50; # minimum number of points + N_MAX = 2000; # maximum number of points + T_DEF = 10.0; # default simulation time + + ## collect useful information about the system + [ncstates,ndstates,NIN,NOUT] = sysdimensions(sys); + TSAMPLE = sysgettsam(sys); + + if (nargin < 3) inp = 1; + elseif (inp < 1 | inp > NIN) error("Argument inp out of range") + endif + + DIGITAL = is_digital(sys); + if (DIGITAL) + NSTATES = ndstates; + if (TSAMPLE < eps) + error("stepimp: sampling time of discrete system too small.") + endif + else NSTATES = ncstates; endif + if (NSTATES < 1) + error("step: pure gain block (n_states < 1), step response is trivial"); + endif + if (nargin < 5) + ## we have to compute the time when the system reaches steady state + ## and the step size + ev = eig(sys2ss(sys)); + if (DIGITAL) + ## perform bilinear transformation on poles in z + for i = 1:NSTATES + pole = ev(i); + if (abs(pole + 1) < 1.0e-10) + ev(i) = 0; + else + ev(i) = 2 / TSAMPLE * (pole - 1) / (pole + 1); + endif + endfor + endif + ## remove poles near zero from eigenvalue array ev + nk = NSTATES; + for i = 1:NSTATES + if (abs(ev(i)) < 1.0e-10) + ev(i) = 0; + nk = nk - 1; + endif + endfor + if (nk == 0) + USE_DEF = 1; + ## printf("##STEPIMP-DEBUG: using defaults.\n"); + else + ev = ev(find(ev)); + x = max(abs(ev)); + t_step = 0.2 * pi / x; + x = min(abs(real(ev))); + t_sim = 5.0 / x; + ## round up + yy = 10^(ceil(log10(t_sim)) - 1); + t_sim = yy * ceil(t_sim / yy); + ## printf("##STEPIMP-DEBUG: nk=%d t_step=%f t_sim=%f\n", + ## nk, t_step, t_sim); + endif + endif + + if (DIGITAL) + ## ---- sampled system + if (nargin == 5) + n = round(n); + if (n < 2) + error("stepimp: n must not be less than 2.") + endif + else + if (nargin == 4) + ## n is unknown + elseif (nargin >= 1) + ## tstop and n are unknown + if (USE_DEF) + tstop = (N_MIN - 1) * TSAMPLE; + else + tstop = t_sim; + endif + endif + n = floor(tstop / TSAMPLE) + 1; + if (n < 2) n = 2; endif + if (n > N_MAX) + n = N_MAX; + printf("Hint: number of samples limited to %d by default.\n", \ + N_MAX); + printf(" ==> increase \"n\" parameter for longer simulations.\n"); + endif + endif + tstop = (n - 1) * TSAMPLE; + t_step = TSAMPLE; + else + ## ---- continuous system + if (nargin == 5) + n = round(n); + if (n < 2) + error("step: n must not be less than 2.") + endif + t_step = tstop / (n - 1); + else + if (nargin == 4) + ## only n in unknown + if (USE_DEF) + n = N_MIN; + t_step = tstop / (n - 1); + else + n = floor(tstop / t_step) + 1; + endif + else + ## tstop and n are unknown + if (USE_DEF) + tstop = T_DEF; + n = N_MIN; + t_step = tstop / (n - 1); + else + tstop = t_sim; + n = floor(tstop / t_step) + 1; + endif + endif + if (n < N_MIN) + n = N_MIN; + t_step = tstop / (n - 1); + endif + if (n > N_MAX) + tstop = (n - 1) * t_step; + t_step = tstop / (N_MAX - 1); + n = N_MAX; + endif + endif + tstop = (n - 1) * t_step; + [jnk,B] = sys2ss(sys); + B = B(:,inp); + sys = c2d(sys, t_step); + endif + ## printf("##STEPIMP-DEBUG: t_step=%f n=%d tstop=%f\n", t_step, n, tstop); + + F = sys.a; + G = sys.b(:,inp); + C = sys.c; + D = sys.d(:,inp); + y = zeros(NOUT, n); + t = linspace(0, tstop, n); + + if (IMPULSE) + if (!DIGITAL && (D'*D > 0)) + error("impulse: D matrix is nonzero, impulse response infinite.") + endif + if (DIGITAL) + y(:,1) = D / t_step; + x = G / t_step; + else + x = B; + y(:,1) = C * x; + x = F * x; + endif + for i = 2:n + y(:,i) = C * x; + x = F * x; + endfor + else + x = zeros(NSTATES, 1); + for i = 1:n + y(:,i) = C * x + D; + x = F * x + G; + endfor + endif + + if(nargout == 0) + ## Plot the information + oneplot(); + gset nogrid + gset nologscale + gset autoscale + gset nokey + clearplot(); + if (gnuplot_has_multiplot) + if (IMPULSE) + gm = zeros(NOUT, 1); + tt = "impulse"; + else + ssys = ss2sys(F, G, C, D, t_step); + gm = dcgain(ssys); + tt = "step"; + endif + ncols = floor(sqrt(NOUT)); + nrows = ceil(NOUT / ncols); + for i = 1:NOUT + subplot(nrows, ncols, i); + title(sprintf("%s: | %s -> %s", tt,sysgetsignals(sys,"in",inp,1), ... + sysgetsignals(sys,"out",i,1))); + if (DIGITAL) + [ts, ys] = stairs(t, y(i,:)); + ts = ts(1:2*n-2)'; ys = ys(1:2*n-2)'; + if (length(gm) > 0) + yy = [ys; gm(i)*ones(size(ts))]; + else + yy = ys; + endif + grid("on"); + xlabel("time [s]"); + ylabel("y(t)"); + plot(ts, yy); + else + if (length(gm) > 0) + yy = [y(i,:); gm(i)*ones(size(t))]; + else + yy = y(i,:); + endif + grid("on"); + xlabel("time [s]"); + ylabel("y(t)"); + plot(t, yy); + endif + endfor + ## leave gnuplot in multiplot mode is bad style + oneplot(); + else + ## plot everything in one diagram + title([tt, " response | ", sysgetsignals(sys,"in",inp,1), ... + " -> all outputs"]); + if (DIGITAL) + stairs(t, y(i,:)); + else + grid("on"); + xlabel("time [s]"); + ylabel("y(t)"); + plot(t, y(i,:)); + endif + endif + y=[]; + t=[]; + endif + ## printf("##STEPIMP-DEBUG: gratulations, successfull completion.\n"); +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/control/system/__abcddims__.m @@ -0,0 +1,38 @@ +## Copyright (C) 1997 Auburn University. All rights reserved. +## +## 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, 59 Temple Place, Suite 330, Boston, MA 02111 USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{y}, @var{my}, @var{ny}] =} abcddims (@var{x}) +## +## Used internally in @code{abcddim}. If @var{x} is a zero-size matrix, +## both dimensions are set to 0 in @var{y}. +## @var{my} and @var{ny} are the row and column dimensions of the result. +## @end deftypefn + +## Author: A. S. Hodel <a.s.hodel@eng.auburn.edu> +## Created: February 1997 + +function [y, my, ny] = abcddims (x) + + y = x; + if(isempty(y)) + y = []; + endif + [my,ny] = size(y); + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/control/system/__syschnamesl__.m @@ -0,0 +1,116 @@ +## Copyright (C) 1996 Auburn University. All rights reserved. +## +## 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, 59 Temple Place, Suite 330, Boston, MA 02111 USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{old_names} =} syschnamesl (@var{olist}, @var{old_names}, @var{inames}, @var{listname}) +## used internally in syschnames +## item olist: index list +## old_names: original list names +## inames: new names +## listname: name of index list +## +## combines the two string lists old_names and inames +## @end deftypefn + +function old_names = syschnamesl (olist, old_names, inames, listname) + + probstr = []; + if( max(olist) > rows(old_names) ) + probstr = ["index list value(s) exceed(s) number of signals (", ... + num2str(rows(old_names)),")"]; + + elseif( length(olist) > rows(inames) ) + probstr = ["index list dimension exceeds number of replacement names (", ... + num2str(rows(inames)),")"]; + + elseif(isempty(olist)) + probstr = []; # do nothing, no changes + + elseif(min(size(olist)) != 1 ) + probstr = "index list must be either a vector or an empty matrix"; + + elseif(max(olist) > rows(old_names)) + probstr = ["max(",listname,")=",num2str(max(olist))," > ", ... + num2str(rows(old_names)),", too big"]; + + elseif(min(olist) < 1) + probstr = ["min(",listname,")=",num2str(min(olist))," < 1, too small"]; + + else + if( length(olist) == 1) + len_in = columns(inames); + len_out = columns(old_names); + + if (len_in < len_out) + inames(1,(len_in+1):(len_out)) = zeros(1,(len_out - len_in)); + endif + + old_names(olist,1:length(inames)) = inames; + elseif(length(olist) > 1) + for ii=1:length(olist) + mystr = inames(ii,:); + len_my = columns(mystr); + len_out = columns(old_names); + + if (len_my < len_out) + mystr(1,(len_my+1):(len_out)) = " "*ones(1,(len_out - len_my)); + len_my = len_out; + endif + + old_names(olist(ii),1:len_my) = mystr; + endfor + endif + endif + if(!isempty(probstr)) + ## the following lines are NOT debugging code! + disp("Problem in syschnames: old names are") + outlist(old_names," ") + disp("new names are") + outlist(inames," ") + disp("list indices are") + disp(olist) + error(sprintf("syschnames: \"%s\" dim=(%d x %d)--\n\t%s\n", ... + listname, rows(olist), columns(olist),probstr)); + endif + + ## change zeros to blanks + if( find(old_names == 0) ) + ## disp("syschnamesl: old_names contains zeros ") + ## old_names + ## disp("/syschnamesl"); + + [ii,jj] = find(old_names == 0); + for idx=1:length(ii) + old_names(ii(idx),jj(idx)) = " "; + endfor + + ## disp("syschnamesl: old_names fixed zeros ") + ## old_names + ## disp("/syschnamesl"); + endif + + ## just in case it's not a string anymore + if( !isstr(old_names) ) + old_names = setstr(old_names); + endif + + ## disp("syschnamesl: exit, old_names=") + ## old_names + ## disp("/syschnamesl: exiting") + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/control/system/__syscont_disc__.m @@ -0,0 +1,54 @@ +## Copyright (C) 1996, 1998 Auburn University. All rights reserved. +## +## 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, 59 Temple Place, Suite 330, Boston, MA 02111 USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{n_tot}, @var{st_c}, @var{st_d}, @var{y_c}, @var{y_d}] =} syscont_disc(@var{sys}) +## Used internally in syscont and sysdisc. +## +## @strong{Inputs} +## @var{ sys} is a system data structure. +## +## @strong{Outputs} +## @table @var +## @item n_tot +## total number of states +## @item st_c +## vector of continuous state indices (empty if none) +## @item st_d +## vector of discrete state indices (empty if none) +## @item y_c +## vector of continuous output indices +## @item y_d +## vector of discrete output indices +## @end table +## @end deftypefn + +## Author: A. S. Hodel <a.s.hodel@eng.auburn.edu> +## Created: February 1997 + +function [n_tot, st_c, st_d, y_c, y_d] = syscont_disc (sys) + + ## get ranges for discrete/continuous states and outputs + [nn,nz,mm,pp,yd] = sysdimensions(sys); + n_tot = nn + nz; + st_c = 1:(nn); + st_d = nn + (1:nz); + y_c = find(yd == 0); # y_c, y_d will be empty if there are none. + y_d = find(yd == 1); + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/control/system/__sysdefioname__.m @@ -0,0 +1,57 @@ +## Copyright (C) 1996, 1998 Auburn University. All rights reserved. +## +## 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, 59 Temple Place, Suite 330, Boston, MA 02111 USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{ioname} =} sysdefioname (@var{n},@var{str} @{,@var{m}@}) +## return default input or output names given @var{n}, @var{str}, @var{m}. +## @var{n} is the final value, @var{str} is the string prefix, and @var{m} +## is start value +## +## used internally, minimal argument checking +## +## @strong{Example} @code{ioname = sysdefioname(5,"u",3)} +## returns the list: +## @example +## ioname = +## ( +## [1] = u_3 +## [2] = u_4 +## [3] = u_5 +## ) +## @end example +## @end deftypefn + +function ioname = sysdefioname (n, str, m) + + if (nargin < 2 | nargin > 3) + usage("ioname = sysdefioname(n,str[,m])"); + endif + + if (nargin == 2) m = min(1,n); endif + + ioname = list(); + jj = 1; + if(n > 0 & m > 0 & m <= n) + for ii = m:n + ioname(ii+1-m) = sprintf("%s_%d",str,ii); + endfor + elseif(m > n) + error("str=%s; start value m=%d > final value n=%d",str,m,n); + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/control/system/__sysdefstname__.m @@ -0,0 +1,42 @@ +## Copyright (C) 1996 Auburn University. All rights reserved. +## +## 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, 59 Temple Place, Suite 330, Boston, MA 02111 USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{stname} =} sysdefstname (@var{n}, @var{nz}) +## return default state names given @var{n}, @var{nz} +## +## used internally, minimal argument checking +## @end deftypefn + +function stname = sysdefstname (n, nz) + + stname = list(); + if(n > 0) + for ii = 1:n + stname(ii) = sprintf("x_%d",ii); + endfor + endif + + ## Set default names for discrete states + if(nz > 0) + for ii = (n+1):(n+nz) + stname(ii) = sprintf("xd_%d",ii); + endfor + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/control/system/__sysgroupn__.m @@ -0,0 +1,57 @@ +## Copyright (C) 1996, 1998 Auburn University. All rights reserved. +## +## 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, 59 Temple Place, Suite 330, Boston, MA 02111 USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{names} =} sysgroupn (@var{names}) +## names = sysgroupn(names) +## Locate and mark duplicate names +## inputs: +## names: list of signal names +## kind: kind of signal name (used for diagnostic message purposes only) +## outputs: +## returns names with unique suffixes added; diagnostic warning +## message is printed to inform the user of the new signal name +## +## used internally in sysgroup and elsewhere. +## @end deftypefn + +function names = sysgroupn (names, kind) + + ## check for duplicate names + l = length(names); + ii = 1; + while(ii <= l-1) + st1 = nth(names,ii); + jj = ii+1; + while ( jj <= l) + st2 = nth(names,jj); + if(strcmp(st1,st2)) + suffix = ["_",num2str(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 + jj = jj+1; + endwhile + ii = ii+1; + endwhile +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/control/system/__tf2sysl__.m @@ -0,0 +1,35 @@ +## Copyright (C) 1996 Auburn University. All rights reserved. +## +## 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, 59 Temple Place, Suite 330, Boston, MA 02111 USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{vec} =} tf2sysl (@var{vec}) +## used internally in @ref{tf2sys}. +## strip leading zero coefficients to get the true polynomial length +## @end deftypefn + +function vec = tf2sysl (vec) + + while( (length(vec) > 1) & (vec(1) == 0) ) + vec = vec(2:length(vec)); + endwhile + + if(vec(1) == 0) + warning("tf2sys: polynomial has no nonzero coefficients!") + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/control/system/__zp2ssg2__.m @@ -0,0 +1,68 @@ +## Copyright (C) 1996, 1998 Auburn University. All rights reserved. +## +## 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, 59 Temple Place, Suite 330, Boston, MA 02111 USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{poly}, @var{rvals}] =} zp2ssg2 (@var{rvals}) +## Used internally in @code{zp2ss} +## Extract 2 values from @var{rvals} (if possible) and construct +## a polynomial with those roots. +## @end deftypefn + +## Author: A. S. Hodel <a.s.hodel@eng.auburn.edu> +## Created: August 1996 + +function [poly, rvals] = zp2ssg2 (rvals) + + ## locate imaginary roots (if any) + cidx = find(imag(rvals)); + + if(!isempty(cidx)) + ## select first complex root, omit from cidx + r1i = cidx(1); r1 = rvals(r1i); cidx = complement(r1i,cidx); + + ## locate conjugate root (must be in cidx list, just in case there's + ## roundoff) + err = abs(rvals(cidx) - r1'); + minerr = min(err); + c2i = find(err == minerr); + r2i = cidx(c2i); + r2 = rvals(r2i); + cidx = complement(r2i,cidx); + + ## don't check for divide by zero, since 0 is not complex. + if(abs(r2 - r1')/abs(r1) > 1e-12) + error(sprintf("r1=(%f,%f); r2=(%f,%f), not conjugates.", ... + real(r1),imag(r1),real(r2),imag(r2))); + endif + + ## complex conjugate pair + poly = [1, -2*real(r1), real(r1)^2+imag(r1)^2]; + else + ## select two roots (they're all real) + r1 = rvals(1); + r2 = rvals(2); + poly = [1, -(r1+r2), (r1*r2)]; + r1i = 1; r2i = 2; + endif + + ## remove roots used + idx = complement([r1i, r2i],1:length(rvals)); + rvals = rvals(idx); + +endfunction +
new file mode 100644 --- /dev/null +++ b/scripts/control/util/__outlist__.m @@ -0,0 +1,80 @@ +## Copyright (C) 1996, 1998 Auburn University. All rights reserved. +## +## 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, 59 Temple Place, Suite 330, Boston, MA 02111 USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} outlist (@var{lmat}@{, @var{tabchar}, @var{yd}, @var{ilist} @}) +## Prints an enumerated list of strings. +## internal use only; minimal argument checking performed +## +## @strong{Inputs} +## @table @var +## @item lmat +## list of strings +## @item tabchar +## tab character (default: none) +## @item yd +## indices of strings to append with the string "(discrete)" +## (used by @var{sysout}; minimal checking of this argument) +## @math{yd = [] } indicates all outputs are continuous +## @item ilist +## index numbers to print with names. +## +## default: @code{1:rows(lmat)} +## @end table +## +## @strong{Outputs} +## prints the list to the screen, numbering each string in order. +## @end deftypefn + +## Author: A. S. Hodel <a.s.hodel@eng.auburn.edu> +## Created: December 1995 + +function str_val = outlist (name_list, tabchar, yd, ilist) + + ## save for restore later + save_empty = empty_list_elements_ok; + empty_list_elements_ok = 1; + + if( nargin < 1 | nargin > 4 ) + usage("str_val = outlist(x[,tabchar,yd,ilist])"); + endif + + m = length(name_list); + if(nargin < 4) ilist = 1:m; endif + if(nargin ==1) + empty_list_elements_ok = 1; + tabchar = ""; + endif + + if(nargin < 3) yd = zeros(1,m); + elseif(isempty(yd)) yd = zeros(1,m); endif + + str_val = ""; + dstr = list(""," (discrete)"); + if((m >= 1) && (is_list(name_list))) + for ii=1:m + str_val = sprintf("%s%s%d: %s%s\n",str_val,tabchar, ilist(ii), ... + nth(name_list,ii),nth(dstr,yd(ii)+1)); + endfor + else + str_val = sprintf("%sNone",tabchar); + endif + + empty_list_elements_ok = save_empty; + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/control/util/__zgpbal__.m @@ -0,0 +1,109 @@ +## Copyright (C) 1996 Auburn University. All rights reserved. +## +## 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, 59 Temple Place, Suite 330, Boston, MA 02111 USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{retsys} =} zgpbal (@var{Asys}) +## +## used internally in @code{tzero}; minimal argument checking performed +## +## implementation of zero computation generalized eigenvalue problem +## balancing method (Hodel and Tiller, Allerton Conference, 1991) +## Based on Ward's balancing algorithm (SIAM J. Sci Stat. Comput., 1981) +## +## zgpbal computes a state/input/output weighting that attempts to +## reduced the range of the magnitudes of the nonzero elements of [a,b,c,d] +## The weighting uses scalar multiplication by powers of 2, so no roundoff +## will occur. +## +## zgpbal should be followed by zgpred +## @end deftypefn + +## References: +## ZGEP: Hodel, "Computation of Zeros with Balancing," 1992, submitted to LAA +## Generalized CG: Golub and Van Loan, "Matrix Computations, 2nd ed" 1989 + +## Author: A. S. Hodel <a.s.hodel@eng.auburn.edu> +## Created: July 24, 1992 +## Conversion to Octave by R. Bruce Tenison July 3, 1994 + +function retsys = zgpbal (Asys) + + if( (nargin != 1) | (!is_struct(Asys))) + usage("retsys = zgpbal(Asys)"); + endif + + Asys = sysupdate(Asys,"ss"); + [a,b,c,d] = sys2ss(Asys); + + [nn,mm,pp] = abcddim(a,b,c,d); + + np1 = nn+1; + nmp = nn+mm+pp; + + ## set up log vector zz, incidence matrix ff + zz = zginit(a,b,c,d); + + ## disp("zgpbal: zginit returns") + ## zz + ## disp("/zgpbal") + + if (norm(zz)) + ## generalized conjugate gradient approach + xx = zgscal(a,b,c,d,zz,nn,mm,pp); + + for i=1:nmp + xx(i) = floor(xx(i)+0.5); + xx(i) = 2.0^xx(i); + endfor + + ## now scale a + ## block 1: a = sigma a inv(sigma) + for i=1:nn + a(i,1:nn) = a(i,1:nn)*xx(i); + a(1:nn,i) = a(1:nn,i)/xx(i); + endfor + ## block 2: b= sigma a phi + for j=1:mm + j1 = j+nn; + b(1:nn,j) = b(1:nn,j)*xx(j1); + endfor + for i=1:nn + b(i,1:mm) = b(i,1:mm)*xx(i); + endfor + for i=1:pp + i1 = i+nn+mm; + ## block 3: c = psi C inv(sigma) + c(i,1:nn) = c(i,1:nn)*xx(i1); + endfor + for j=1:nn + c(1:pp,j) = c(1:pp,j)/xx(j); + endfor + ## block 4: d = psi D phi + for j=1:mm + j1 = j+nn; + d(1:pp,j) = d(1:pp,j)*xx(j1); + endfor + for i=1:pp + i1 = i + nn + mm; + d(i,1:mm) = d(i,1:mm)*xx(i1); + endfor + endif + + retsys = ss2sys(a,b,c,d); +endfunction +