Mercurial > hg > octave-terminal
changeset 7126:4a375de63f66
[project @ 2007-11-08 03:44:14 by jwe]
line wrap: on
line diff
--- a/scripts/control/base/__bodquist__.m +++ b/scripts/control/base/__bodquist__.m @@ -61,36 +61,36 @@ endif ## check each argument to see if it's in the correct form - if (!isstruct(sys)) - error("sys must be a system data structure"); + if (! isstruct (sys)) + error ("sys must be a system data structure"); endif ## let __freqresp__ determine w if it's not already given - USEW = freqchkw(w); + USEW = freqchkw (w); ## get initial dimensions (revised below if sysprune is called) - [nn,nz,mm,pp ] = sysdimensions(sys); + [nn, nz, mm, pp] = sysdimensions (sys); ## check for an output vector and to see whether it`s correct - if (!isempty(outputs)) - if (isempty(inputs)) + if (! isempty (outputs)) + if (isempty (inputs)) inputs = 1:mm; # use all inputs - warning([rname,": outputs specified but not inputs"]); - elseif(is_signal_list(inputs) | ischar(inputs)) - inputs = sysidx(sys,"in",inputs); + warning ("%s: outputs specified but not inputs", rname); + elseif (is_signal_list (inputs) || ischar (inputs)) + inputs = sysidx (sys, "in", inputs); endif - if(is_signal_list(outputs) | ischar(outputs)) - outputs = sysidx(sys,"out",outputs); + if (is_signal_list (outputs) || ischar (outputs)) + outputs = sysidx (sys, "out", outputs); end - sys = sysprune(sys,outputs,inputs); - [nn,nz,mm,pp ] = sysdimensions(sys); + 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 = zp(zer,pol,k,tsam,inname,outname); + if (is_siso (sys) && strcmp (sysgettype (sys), "ss")) + [zer, pol, k, tsam, inname, outname] = sys2zp (sys); + sys = zp (zer, pol, k, tsam, inname, outname); endif ## get system frequency response @@ -98,46 +98,49 @@ phase = arg(f)*180.0/pi; - if(!USEW) + 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 = find(pd > dphase); + while (pcnt) + pd = abs (diff (phase)); # phase variation + pdbig = find (pd > dphase); + + ## relative variation + lp = length (f); + lp1 = lp-1; - lp = length(f); lp1 = lp-1; # relative variation - fd = abs(diff(f)); - fm = max(abs([f(1:lp1); f(2:lp)])); - fdbig = find(fd > fm/10); + fd = abs (diff (f)); + fm = max (abs ([f(1:lp1); f(2:lp)])); + fdbig = find (fd > fm/10); - bigpts = union(fdbig, pdbig); + bigpts = union (fdbig, pdbig); - if(isempty(bigpts) ) + 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)); + 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)); + 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 + wnew = wnew(find (wnew != 0)); + wnew = sort (wnew); + wnew = create_set (wnew); + if (isempty (wnew)) # all small crossovers pcnt = 0; else ## get new freq resp points, combine with old, and sort. @@ -153,12 +156,12 @@ endif ## ensure unique frequency values - [w,idx] = sort(w); + [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_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);
--- a/scripts/control/base/__freqresp__.m +++ b/scripts/control/base/__freqresp__.m @@ -50,82 +50,80 @@ ## SYS_INTERNAL accesses members of system data structure ## Check Args - if ((nargin < 2) || (nargin > 4)) + if (nargin < 2 || nargin > 4) print_usage (); - elseif (USEW & (nargin < 3) ) + 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"); + elseif (USEW && isempty (w)) + warning ("USEW = 1 but w is empty; setting USEW=0"); USEW = 0; endif - DIGITAL = is_digital(sys); + DIGITAL = is_digital (sys); ## compute default w if needed - if(!USEW) - if(is_siso(sys)) - sys = sysupdate(sys,"zp"); - [zer,pol] = sys2zp(sys); + if (! USEW) + if (is_siso (sys)) + sys = sysupdate (sys, "zp"); + [zer, pol] = sys2zp (sys); else - zer = tzero(sys); - pol = eig(sys2ss(sys)); + 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); + [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 + 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)); + if (DIGITAL) + jw = exp (i*w*sysgettsam(sys)); else jw = i*w; endif - [nn,nz,mm,pp] = sysdimensions(sys); + [nn, nz, mm, pp] = sysdimensions (sys); ## now compute the frequency response - divide by zero yields a warning - if (strcmp(sysgettype(sys),"zp")) + 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)); + [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) + for ii = (l1+1):length(pol) ff = ff ./ (jw - pol(ii)); endfor ff = ff*sysk; - - elseif (strcmp(sysgettype(sys),"tf")) + elseif (strcmp (sysgettype (sys), "tf")) ## transfer function form - [num,den] = sys2tf(sys); - ff = polyval(num,jw)./polyval(den,jw); + [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); + [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; + 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); + [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); + 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)); + w = reshape (w, 1, length (w)); + ff = reshape (ff, 1, length (ff)); endfunction
--- a/scripts/control/base/__stepimp__.m +++ b/scripts/control/base/__stepimp__.m @@ -269,5 +269,5 @@ y = []; t = []; endif - ## printf("##STEPIMP-DEBUG: gratulations, successfull completion.\n"); + endfunction
--- a/scripts/control/base/are.m +++ b/scripts/control/base/are.m @@ -70,10 +70,11 @@ if (nargin == 3 || nargin == 4) if (nargin == 4) - if (! (strcmp (opt, "N") || strcmp (opt, "P") ... - || strcmp (opt, "S") || strcmp (opt, "B") ... - || strcmp (opt, "n") || strcmp (opt, "p") ... - || strcmp (opt, "s") || strcmp (opt, "b"))) + if (! (ischar (opt) + && (strcmp (opt, "N") || strcmp (opt, "P") + || strcmp (opt, "S") || strcmp (opt, "B") + || strcmp (opt, "n") || strcmp (opt, "p") + || strcmp (opt, "s") || strcmp (opt, "b")))) warning ("are: opt has an invalid value; setting to B"); opt = "B"; endif
--- a/scripts/control/base/ctrb.m +++ b/scripts/control/base/ctrb.m @@ -47,23 +47,23 @@ if (nargin == 2) a = sys; - elseif (nargin == 1 && isstruct(sys)) - sysupdate(sys,"ss"); - [a,b] = sys2ss(sys); + elseif (nargin == 1 && isstruct (sys)) + sysupdate (sys, "ss"); + [a, b] = sys2ss (sys); else print_usage (); endif - if (!is_abcd(a,b)) + if (! is_abcd (a, b)) Qs = []; else ## no need to check dimensions, we trust is_abcd(). - [na, ma] = size(a); + [na, ma] = size (a); ## using imb avoids name conflict with the "mb" function - [inb, imb] = size(b); - Qs = zeros(na, ma*imb); + [inb, imb] = size (b); + Qs = zeros (na, ma*imb); for i = 1:na - Qs(:, (i-1)*imb+1:i*imb) = b; + Qs(:,(i-1)*imb+1:i*imb) = b; b = a * b; endfor endif
--- a/scripts/control/base/damp.m +++ b/scripts/control/base/damp.m @@ -35,51 +35,53 @@ ## assume a continuous system DIGITAL = 0; - if(nargin < 1 || nargin > 2) + + if (nargin < 1 || nargin > 2) print_usage (); endif - if(isstruct(p)) + + if (isstruct (p)) if (nargin != 1) error("damp: when p is a system, tsamp parameter is not allowed."); endif - [aa, b, c, d, t_samp] = sys2ss(p); - DIGITAL = is_digital(p); + [aa, b, c, d, t_samp] = sys2ss (p); + DIGITAL = is_digital (p); else aa = p; if (nargin == 2) - DIGITAL = 1; - t_samp = tsam; + DIGITAL = 1; + t_samp = tsam; endif endif - if (!issquare(aa)) - error("damp: Matrix p is not square.") + if (! issquare (aa)) + error ("damp: Matrix p is not square.") endif if (DIGITAL && t_samp <= 0.0) - error("damp: Sampling time tsam must not be <= 0.") + error ("damp: Sampling time tsam must not be <= 0.") endif ## all checks done. - e = eig(aa); - [n, m] = size(aa); + e = eig (aa); + [n, m] = size (aa); if (DIGITAL) - printf(" (discrete system with sampling time %f)\n", t_samp); + printf (" (discrete system with sampling time %f)\n", t_samp); endif - printf("............... Eigenvalue ........... Damping Frequency\n"); - printf("--------[re]---------[im]--------[abs]----------------------[Hz]\n"); + printf ("............... Eigenvalue ........... Damping Frequency\n"); + printf ("--------[re]---------[im]--------[abs]----------------------[Hz]\n"); for i = 1:n pole = e(i); cpole = pole; if (DIGITAL) - cpole = log(pole) / t_samp; + cpole = log (pole) / t_samp; endif - d0 = -cos(atan2(imag(cpole), real(cpole))); - f0 = 0.5 / pi * abs(cpole); - if (abs(imag(cpole)) < eps) - printf("%12f --- %12f %10f %12f\n", - real(pole), abs(pole), d0, f0); + d0 = -cos (atan2 (imag (cpole), real (cpole))); + f0 = 0.5 / pi * abs (cpole); + if (abs (imag (cpole)) < eps) + printf ("%12f --- %12f %10f %12f\n", + real (pole), abs (pole), d0, f0); else - printf("%12f %12f %12f %10f %12f\n", - real(pole), imag(pole), abs(pole), d0, f0); + printf ("%12f %12f %12f %10f %12f\n", + real (pole), imag (pole), abs (pole), d0, f0); endif endfor
--- a/scripts/control/base/dare.m +++ b/scripts/control/base/dare.m @@ -78,9 +78,11 @@ function x = dare (a, b, q, r, opt) - if (nargin == 4 | nargin == 5) + if (nargin == 4 || nargin == 5) if (nargin == 5) - if (opt != "N" || opt != "P" || opt != "S" || opt != "B") + if (! (ischar (opt) + && (strcmp (opt, "N") || strcmp (opt, "P") + || strcmp (opt, "S") || strcmp (opt, "B")))) warning ("dare: opt has an invalid value -- setting to B"); opt = "B"; endif @@ -88,22 +90,20 @@ opt = "B"; endif - if ((p = issquare (q)) == 0) q = q'*q; endif ##Checking positive definiteness - if (isdefinite(r)<=0) - error("dare: r not positive definite"); + if (isdefinite (r) <= 0) + error ("dare: r not positive definite"); end - if (isdefinite(q)<0) - error("dare: q not positive semidefinite"); + if (isdefinite (q) < 0) + error ("dare: q not positive semidefinite"); end - ## Check r dimensions. - [n,m] = size(b); + [n, m] = size (b); if ((m1 = issquare (r)) == 0) error ("dare: r is not square"); elseif (m1 != m) @@ -112,8 +112,11 @@ s1 = [a, zeros(n) ; -q, eye(n)]; s2 = [eye(n), (b/r)*b' ; zeros(n), a']; - [c,d,s1,s2] = balance(s1,s2,opt); - [aa,bb,u,lam] = qz(s1,s2,"S"); + + [c, d, s1, s2] = balance (s1, s2, opt); + + [aa, bb, u, lam] = qz (s1, s2, "S"); + u = d*u; n1 = n+1; n2 = 2*n;
--- a/scripts/control/base/dcgain.m +++ b/scripts/control/base/dcgain.m @@ -30,26 +30,30 @@ function gm = dcgain (sys, tol) - if((nargin < 1) || (nargin > 2) || (nargout > 1)) + if (nargin < 1 || nargin > 2 || nargout > 1) print_usage (); endif - if(!isstruct(sys)) - error("dcgain: first argument is not a system data structure.") + if (! isstruct (sys)) + error ("dcgain: first argument is not a system data structure.") endif - sys = sysupdate(sys, "ss"); - [aa,bb,cc,dd] = sys2ss(sys); - if (is_digital(sys)) aa = aa - eye(size(aa)); endif - if (nargin == 1) tol = 1.0e-10; endif - r = rank(aa, tol); - if (r < rows(aa)) + sys = sysupdate (sys, "ss"); + [aa, bb, cc, dd] = sys2ss (sys); + if (is_digital (sys)) + aa = aa - eye (size (aa)); + endif + if (nargin == 1) + tol = 1.0e-10; + endif + r = rank (aa, tol); + if (r < rows (aa)) gm = []; else gm = -cc / aa * bb + dd; endif - if(!is_stable(sys)) - [nn,nz,mm,pp] = sysdimensions(sys); - warning("dcgain: unstable system; dimensions [nc=%d,nz=%d,mm=%d,pp=%d]", ... - nn,nz,mm,pp); + if (! is_stable (sys)) + [nn, nz, mm, pp] = sysdimensions (sys); + warning ("dcgain: unstable system; dimensions: nc=%d, nz=%d, mm=%d, pp=%d", + nn, nz, mm, pp); endif endfunction
--- a/scripts/control/base/dgram.m +++ b/scripts/control/base/dgram.m @@ -62,7 +62,7 @@ function m = dgram (a, b) if (nargin != 2) - print_usage(); + print_usage (); endif ## let dlyap do the error checking...
--- a/scripts/control/base/dlqr.m +++ b/scripts/control/base/dlqr.m @@ -128,16 +128,17 @@ qo = q; endif - ## Checking stabilizability and detectability (dimensions are checked inside these calls) + ## Checking stabilizability and detectability (dimensions are checked + ## inside these calls). tol = 200*eps; - if (is_stabilizable (ao, b,tol,1) == 0) + if (is_stabilizable (ao, b, tol, 1) == 0) error ("dlqr: (a,b) not stabilizable"); endif - dflag = is_detectable (ao, qo, tol,1); - if ( dflag == 0) + dflag = is_detectable (ao, qo, tol, 1); + if (dflag == 0) warning ("dlqr: (a,q) not detectable"); - elseif ( dflag == -1) - error("dlqr: (a,q) has non minimal modes near unit circle"); + elseif (dflag == -1) + error ("dlqr: (a,q) has non minimal modes near unit circle"); end ## Compute the Riccati solution @@ -145,6 +146,4 @@ k = (r+b'*p*b)\(b'*p*a + s'); e = eig (a - b*k); - endfunction -
--- a/scripts/control/base/dre.m +++ b/scripts/control/base/dre.m @@ -97,72 +97,79 @@ function [tvals, Plist] = dre (sys, Q, R, Qf, t0, tf, Ptol, maxits) - if(nargin < 6 | nargin > 8 | nargout != 2) + if (nargin < 6 || nargin > 8 || nargout != 2) print_usage (); - elseif(!isstruct(sys)) - error("sys must be a system data structure") - elseif(is_digital(sys)) - error("sys must be a continuous time system") - elseif(!ismatrix(Q) | !ismatrix(R) | !ismatrix(Qf)) - error("Q, R, and Qf must be matrices."); - elseif(!isscalar(t0) | !isscalar(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(!isscalar(Ptol)) error("Ptol must be a scalar"); - elseif(Ptol <= 0) error("Ptol must be positive"); + elseif (! isstruct (sys)) + error ("sys must be a system data structure") + elseif (is_digital (sys)) + error ("sys must be a continuous time system") + elseif (! ismatrix (Q) || ! ismatrix (R) || ! ismatrix (Qf)) + error ("Q, R, and Qf must be matrices"); + elseif (! isscalar (t0) || ! isscalar (tf)) + error ("t0 and tf must be scalars") + elseif (t0 >= tf) + error ("t0=%e >= tf=%e", t0, tf); + elseif (nargin < 7) + Ptol = 0.1; + elseif (! isscalar (Ptol)) + error ("Ptol must be a scalar"); + elseif (Ptol <= 0) + error ("Ptol must be positive"); endif - if(nargin < 8) maxits = 10; - elseif(!isscalar(maxits)) error("maxits must be a scalar"); - elseif(maxits <= 0) error("maxits must be positive"); + if (nargin < 8) + maxits = 10; + elseif (! isscalar (maxits)) + error ("maxits must be a scalar"); + elseif (maxits <= 0) + error ("maxits must be positive"); endif - maxits = ceil(maxits); + maxits = ceil (maxits); - [aa,bb] = sys2ss(sys); - nn = sysdimensions(sys,"cst"); - mm = sysdimensions(sys,"in"); - pp = sysdimensions(sys,"out"); + [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); + 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); + 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) + while (! done) done = 1; # assume this pass will do the job ## sort time values in reverse order - tvals = -sort(-tvals); - tvlen = length(tvals); + 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 ; Plist{ii-1} ]; + for ii = 2:tvlen + uv_i_minus_1 = [In; Plist{ii-1}]; delta_t = tvals(ii-1) - tvals(ii); - uv = expm(-H*delta_t)*uv_i_minus_1; + 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(Plist{ii} - Plist{ii-1})/norm(Plist{ii}); - maxerr = max(maxerr,Perr); - if(Perr > Ptol) - new_t = mean(tvals([ii,ii-1])); + Perr = norm (Plist{ii} - Plist{ii-1})/norm(Plist{ii}); + maxerr = max (maxerr,Perr); + if (Perr > Ptol) + new_t = mean (tvals([ii,ii-1])); tvals = [tvals, new_t]; done = 0; endif @@ -170,11 +177,11 @@ ## check number of iterations maxits = maxits - 1; - done = done+(maxits==0); + 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); + if (maxerr > Ptol) + warning ("dre: exiting with %d points, max rel chg. = %e, Ptol = %e", + tvlen, maxerr, Ptol); tvals = tvals(1:length(Plist)); endif
--- a/scripts/control/base/impulse.m +++ b/scripts/control/base/impulse.m @@ -56,11 +56,7 @@ function [y, t] = impulse (sys, inp, tstop, n) - if ((nargin < 1) || (nargin > 4)) - print_usage (); - endif - - if (nargout > 2) + if (nargin < 1 || nargin > 4 || nargout > 2) print_usage (); endif
--- a/scripts/control/base/lqe.m +++ b/scripts/control/base/lqe.m @@ -92,7 +92,7 @@ function [k, p, e] = lqe (a, g, c, sigw, sigv, zz) - if ( (nargin != 5) && (nargin != 6)) + if (nargin != 5 && nargin != 6) error ("lqe: invalid number of arguments"); endif
--- a/scripts/control/base/lqg.m +++ b/scripts/control/base/lqg.m @@ -70,89 +70,98 @@ function [K, Q1, P1, Ee, Er] = lqg (sys, Sigw, Sigv, Q, R, input_list) - if ( (nargin < 5) | (nargin > 6)) + if (nargin < 5 || nargin > 6) print_usage (); - - elseif(!isstruct(sys) ) - error("sys must be in system data structure"); + elseif (! isstruct (sys)) + error ("sys must be in system data structure"); endif - DIG = is_digital(sys); - [A,B,C,D,tsam,n,nz,stname,inname,outname] = sys2ss(sys); - [n,nz,nin,nout] = sysdimensions(sys); - if(nargin == 5) + DIG = is_digital (sys); + [A, B, C, D, tsam, n, nz, stname, inname, outname] = sys2ss (sys); + [n, nz, nin, nout] = sysdimensions (sys); + if (nargin == 5) ## construct default input_list input_list = (columns(Sigw)+1):nin; endif - if( !(n+nz) ) - error(["lqg: 0 states in system"]); + if (! (n+nz)) + error("lqg: 0 states in system"); - elseif(nin != columns(Sigw)+ columns(R)) - error(["lqg: sys has ",num2str(nin)," inputs, dim(Sigw)=", ... - num2str(columns(Sigw)),", dim(u)=",num2str(columns(R))]) + elseif (nin != columns (Sigw) + columns (R)) + error ("lqg: sys has %d inputs, dim(Sigw)=%d, dim(u)=%d", + nin, columns (Sigw), columns (R)); - elseif(nout != columns(Sigv)) - error(["lqg: sys has ",num2str(nout)," outputs, dim(Sigv)=", ... - num2str(columns(Sigv)),")"]) + elseif (nout != columns (Sigv)) + error ("lqg: sys has %d outputs, dim(Sigv)=%d", nout, columns (Sigv)); endif ## check for names of signals - if(is_signal_list(input_list) | ischar(input_list)) - input_list = sysidx(sys,"in",input_list); + if (is_signal_list (input_list) || ischar (input_list)) + input_list = sysidx (sys, "in", input_list); endif - if(length(input_list) != columns(R)) - error(["lqg: length(input_list)=",num2str(length(input_list)), ... - ", columns(R)=", num2str(columns(R))]); + if (length (input_list) != columns (R)) + error ("lqg: length(input_list)=%d, columns(R)=%d", + length (input_list), columns (R)); + endif + + if (! issquare (Sigw)) + error ("lqg: Sigw is not square"); endif - varname = {"Sigw","Sigv","Q","R"}; - for kk=1:length(varname); - eval(sprintf("chk = issquare(%s);",varname{kk})); - if(! chk ) error("lqg: %s is not square",varname{kk}); endif - endfor + if (! issquare (Sigv)) + error ("lqg: Sigv is not square"); + endif + + if (! issquare (Q)) + error ("lqg: Q is not square"); + endif + + if (! issquare (R)) + error ("lqg: Q is not square"); + endif ## permute (if need be) - if(nargin == 6) - all_inputs = sysreorder(nin,input_list); + if (nargin == 6) + all_inputs = sysreorder (nin, input_list); B = B(:,all_inputs); - inname = inname(all_inputs); + inname = inname (all_inputs); endif ## put parameters into correct variables - m1 = columns(Sigw); + m1 = columns (Sigw); m2 = m1+1; G = B(:,1:m1); B = B(:,m2:nin); ## now we can just do the design; call dlqr and dlqe, since all matrices ## are not given in Cholesky factor form (as in h2syn case) - if(DIG) - [Ks, P1, Er] = dlqr(A,B,Q,R); - [Ke, Q1, jnk, Ee] = dlqe(A,G,C,Sigw,Sigv); + if (DIG) + [Ks, P1, Er] = dlqr (A, B, Q, R); + [Ke, Q1, jnk, Ee] = dlqe (A, G, C, Sigw, Sigv); else - [Ks, P1, Er] = lqr(A,B,Q,R); - [Ke, Q1, Ee] = lqe(A,G,C,Sigw,Sigv); + [Ks, P1, Er] = lqr (A, B, Q, R); + [Ke, Q1, Ee] = lqe (A, G, C, Sigw, Sigv); endif Ac = A - Ke*C - B*Ks; Bc = Ke; Cc = -Ks; - Dc = zeros(rows(Cc),columns(Bc)); + Dc = zeros (rows (Cc), columns (Bc)); ## fix state names - stname1 = strappend(stname,"_e"); + stname1 = strappend (stname, "_e"); ## fix controller output names - outname1 = strappend(inname(m2:nin),"_K"); + outname1 = strappend (inname(m2:nin), "_K"); ## fix controller input names - inname1 = strappend(outname,"_K"); + inname1 = strappend (outname, "_K"); - if(DIG) - K = ss(Ac,Bc,Cc,Dc,tsam,n,nz,stname1,inname1,outname1,1:rows(Cc)); + if (DIG) + K = ss (Ac, Bc, Cc, Dc, tsam, n, nz, stname1, inname1, outname1, + 1:rows(Cc)); else - K = ss(Ac,Bc,Cc,Dc,tsam,n,nz,stname,inname1,outname1); + K = ss (Ac, Bc, Cc, Dc, tsam, n, nz, stname, inname1, outname1); endif endfunction
--- a/scripts/control/base/lqr.m +++ b/scripts/control/base/lqr.m @@ -119,7 +119,7 @@ ## disp("lqr: entry"); - if ((nargin != 4) && (nargin != 5)) + if (nargin != 4 && nargin != 5) error ("lqr: invalid number of arguments"); endif @@ -135,19 +135,19 @@ endif ## Check q. - if ( ((n1 = issquare (q)) == 0) || (n1 != n)) + if ((n1 = issquare (q)) == 0 || n1 != n) error ("lqr: q must be square and conformal with a"); endif ## Check r. - if ( ((m1 = issquare(r)) == 0) || (m1 != m)) + if ((m1 = issquare(r)) == 0 || m1 != m) error ("lqr: r must be square and conformal with column dimension of b"); endif ## Check if n is there. if (nargin == 5) [n1, m1] = size (s); - if ( (n1 != n) || (m1 != m)) + if (n1 != n || m1 != m) error ("lqr: z must be identically dimensioned with b"); endif @@ -162,7 +162,7 @@ ## Check that q, (r) are symmetric, positive (semi)definite - if (issymmetric (q) && issymmetric (r) ... + if (issymmetric (q) && issymmetric (r) && all (eig (q) >= 0) && all (eig (r) > 0)) p = are (ao, (b/r)*b', qo); k = r\(b'*p + s'); @@ -171,5 +171,4 @@ error ("lqr: q (r) must be symmetric positive (semi) definite"); endif - ## disp("lqr: exit"); endfunction
--- a/scripts/control/base/lsim.m +++ b/scripts/control/base/lsim.m @@ -40,38 +40,40 @@ function [y, x] = lsim (sys, u, t, x0) - if((nargin < 3)||(nargin > 4)) + if (nargin < 3 || nargin > 4) print_usage (); endif - if(!isstruct(sys)) - error("sys must be in system data structure"); + if (! isstruct (sys)) + error ("sys must be in system data structure"); endif - sys = sysupdate(sys,"ss"); + sys = sysupdate (sys,"ss"); - [ncstates, ndstates, nin, nout] = sysdimensions(sys); - [a,b,c,d] = sys2ss(sys); + [ncstates, ndstates, nin, nout] = sysdimensions (sys); + [a, b, c, d] = sys2ss (sys); - if (nargin == 3) x0 = zeros(columns(a),1); endif + if (nargin == 3) + x0 = zeros (columns (a), 1); + endif - if(rows(u) ~= length(t)) - error("lsim: There should be an input value (row) for each time instant"); + if (rows (u) != length (t)) + error ("lsim: There should be an input value (row) for each time instant"); endif - if(columns(u) ~= columns(d)) - error("lsim: U and d should have the same number of inputs"); + if (columns (u) != columns (d)) + error ("lsim: U and d should have the same number of inputs"); endif - if(columns(x0) > 1) - error("lsim: Initial condition vector should have only one column"); + if (columns (x0) > 1) + error ("lsim: Initial condition vector should have only one column"); endif - if(rows(x0) > rows(a)) - error("lsim: Initial condition vector is too large"); + if (rows (x0) > rows p(a)) + error ("lsim: Initial condition vector is too large"); endif Ts = 0; t(2)-t(1); u=u'; - n = max(size(t)); + n = max (size (t)); for ii = 1:(n-1) @@ -81,8 +83,8 @@ if (abs (t(ii+1) - t(ii) - Ts) > 10 * eps) Ts = t(ii+1) - t(ii); ## [F,G] = c2d(a,b,Ts); - dsys = c2d(sys, Ts); - [F,G] = sys2ss(dsys); + dsys = c2d (sys, Ts); + [F, G] = sys2ss (dsys); endif x(:,ii) = x0; @@ -93,9 +95,10 @@ x(:,n) = x0; y = c*x + d*u; - if(nargout == 0) - plot(t,y); - y=[]; - x=[]; + if (nargout == 0) + plot (t, y); + y = []; + x = []; endif + endfunction
--- a/scripts/control/base/ltifr.m +++ b/scripts/control/base/ltifr.m @@ -59,53 +59,56 @@ function out = ltifr (a, b, w) - if ((nargin < 2) || (nargin > 3)) + if (nargin < 2 || nargin > 3) error("incorrect number of input arguments"); endif if (nargin == 2) sys = a; w = b; - if(!isstruct(sys)) - error("two arguments: 1st must be a system data structure"); + if(! isstruct (sys)) + error ("two arguments: 1st must be a system data structure"); endif - if (!isvector(w)) - error("w must be a vector"); + if (! isvector (w)) + error ("w must be a vector"); endif - [nn,nz,mm,pp] = sysdimensions(sys); - if(mm != 1) error("sys has %d > 1 inputs",mm); endif + [nn, nz, mm, pp] = sysdimensions (sys); + if (mm != 1) + error("sys has %d > 1 inputs", mm); + endif - [a,b] = sys2ss(sys); + [a, b] = sys2ss (sys); else - if (columns(a) != rows(b)), - error("ltifr: A(%dx%d), B(%dx%d) not compatibly dimensioned", ... - rows(a), columns(a), rows(b), columns(b)); + if (columns (a) != rows (b)), + error ("ltifr: A(%dx%d), B(%dx%d) not compatibly dimensioned", + rows (a), columns(a), rows(b), columns(b)); endif - if(columns(b) != 1) - error("ltifr: b(%dx%d) must be a single column vector", ... - rows(b),columns(b)); + if (columns (b) != 1) + error ("ltifr: b(%dx%d) must be a single column vector", + rows(b), columns(b)); endif - if (!issquare(a)) - error("ltifr: A(%dx$d) must be square.",rows(a),columns(a)) + if (! issquare (a)) + error ("ltifr: A(%dx$d) must be square", rows(a), columns(a)) endif endif - if (!isvector(w)) - error("w must be a vector"); + if (! isvector (w)) + error ("w must be a vector"); endif - ey = eye(size(a)); - lw = length(w); - out = ones(columns(a),lw); + ey = eye (size (a)); + lw = length (w); + out = ones (columns (a), lw); - for ii=1:lw, + for ii = 1:lw, out(:,ii) = (w(ii)*ey-a)\b; endfor + endfunction
--- a/scripts/control/base/nichols.m +++ b/scripts/control/base/nichols.m @@ -100,7 +100,7 @@ [f, w, sys] = __bodquist__ (sys, w, outputs, inputs, "nichols"); - [stname,inname,outname] = sysgetsignals (sys); + [stname, inname, outname] = sysgetsignals (sys); systsam = sysgettsam (sys); ## Get the magnitude and phase of f.
--- a/scripts/control/base/nyquist.m +++ b/scripts/control/base/nyquist.m @@ -168,7 +168,7 @@ th = atan2 (real (df), imag (df)) * 180 / pi; ## get angle at fmax - if (fi == length(f)) + if (fi == length (f)) fi = fi-1; endif thm = th(fi);
--- a/scripts/control/base/obsv.m +++ b/scripts/control/base/obsv.m @@ -54,19 +54,19 @@ if (nargin == 2) a = sys; elseif (nargin == 1 && isstruct(sys)) - sysupdate(sys,"ss"); - [a,b,c] = sys2ss(sys); + sysupdate (sys, "ss"); + [a, b, c] = sys2ss (sys); else print_usage (); endif - if (!is_abcd(a,c')) + if (! is_abcd (a, c')) Qb = []; else ## no need to check dimensions, we trust is_abcd(). - [na, ma] = size(a); - [nc, mc] = size(c); - Qb = zeros(na*nc, ma); + [na, ma] = size (a); + [nc, mc] = size (c); + Qb = zeros (na*nc, ma); for i = 1:na Qb((i-1)*nc+1:i*nc, :) = c; c = c * a;
--- a/scripts/control/base/place.m +++ b/scripts/control/base/place.m @@ -51,7 +51,7 @@ ## check arguments - if(! isstruct (sys)) + if (! isstruct (sys)) error ("sys must be in system data structure format (see ss)"); endif sys = sysupdate (sys, "ss"); # make sure it has state space form up to date @@ -66,7 +66,7 @@ is_digital (sys); [n, nz, m, p] = sysdimensions (sys); nx = n+nz; # already checked that it's not a mixed system. - if(m != 1) + if (m != 1) error ("sys has %d inputs; need only 1", m); endif
--- a/scripts/control/base/rlocus.m +++ b/scripts/control/base/rlocus.m @@ -66,22 +66,22 @@ endif ## Convert the input to a transfer function if necessary - [num,den] = sys2tf(sys); # extract numerator/denom polyomials - lnum = length(num); - lden = length(den); - # equalize length of num, den polynomials - if(lden < 2) - error("system has no poles"); - elseif(lnum < lden) + [num, den] = sys2tf (sys); # extract numerator/denom polyomials + lnum = length (num); + lden = length (den); + ## equalize length of num, den polynomials + if (lden < 2) + error ("system has no poles"); + elseif (lnum < lden) num = [zeros(1,lden-lnum), num]; # so that derivative is shortened by one endif - olpol = roots(den); - olzer = roots(num); - nas = lden -lnum; # number of asymptotes + olpol = roots (den); + olzer = roots (num); + nas = lden - lnum; # number of asymptotes maxk = 0; - if(nas > 0) - cas = ( sum(olpol) - sum(olzer) )/nas; + if (nas > 0) + cas = (sum (olpol) - sum (olzer)) / nas; angles = (2*[1:nas]-1)*pi/nas; # printf("rlocus: there are %d asymptotes centered at %f\n", nas, cas); else @@ -91,33 +91,33 @@ # compute real axis break points and corresponding gains - dnum=polyderiv(num); - dden=polyderiv(den); - brkp = conv(den, dnum) - conv(num, dden); - real_ax_pts = roots(brkp); - real_ax_pts = real_ax_pts(find(imag(real_ax_pts) == 0)); - k_break = -polyval(den,real_ax_pts) ./ polyval(num, real_ax_pts); - idx = find(k_break >= 0); + dnum = polyderiv (num); + dden = polyderiv (den); + brkp = conv (den, dnum) - conv (num, dden); + real_ax_pts = roots (brkp); + real_ax_pts = real_ax_pts(find (imag (real_ax_pts) == 0)); + k_break = -polyval (den, real_ax_pts) ./ polyval (num, real_ax_pts); + idx = find (k_break >= 0); k_break = k_break(idx); real_ax_pts = real_ax_pts(idx); - if(!isempty(k_break)) - maxk = max(max(k_break),maxk); + if (! isempty (k_break)) + maxk = max (max (k_break), maxk); endif - if(nas == 0) - maxk = max(1, 2*maxk); % get at least some root locus + if (nas == 0) + maxk = max (1, 2*maxk); # get at least some root locus else - # get distance from breakpoints, poles, and zeros to center of asymptotes - dmax = 3*max(abs( [vec(olzer); vec(olpol); vec(real_ax_pts)] - cas )); - if(dmax == 0) + ## get distance from breakpoints, poles, and zeros to center of asymptotes + dmax = 3*max (abs ([vec(olzer); vec(olpol); vec(real_ax_pts)] - cas)); + if (dmax == 0) dmax = 1; endif # get gain for dmax along each asymptote, adjust maxk if necessary - svals = cas + dmax*exp(j*angles); - kvals = -polyval(den,svals) ./ polyval(num, svals); - maxk = max(maxk, max(real(kvals))); - end + svals = cas + dmax * exp (j*angles); + kvals = -polyval (den, svals) ./ polyval (num, svals); + maxk = max (maxk, max (real (kvals))); + endif ## check for input arguments: if (nargin > 2) @@ -129,8 +129,8 @@ maxk = max_k; endif if (nargin > 1) - if(increment <= 0) - error("increment must be positive"); + if (increment <= 0) + error ("increment must be positive"); else ngain = (maxk-mink)/increment; endif @@ -139,48 +139,49 @@ endif ## vector of gains - ngain = max(30,ngain); - gvec = linspace(mink,maxk,ngain); - if(length(k_break)) - gvec = sort([gvec, vec(k_break)']); + ngain = max (30, ngain); + gvec = linspace (mink, maxk, ngain); + if (length (k_break)) + gvec = sort ([gvec, vec(k_break)']); endif ## Find the open loop zeros and the initial poles - rlzer = roots(num); + rlzer = roots (num); ## update num to be the same length as den - lnum = length(num); - if(lnum < lden) + lnum = length (num); + if (lnum < lden) num = [zeros(1,lden - lnum),num]; endif ## compute preliminary pole sets - nroots = lden-1; - for ii=1:ngain + nroots = lden - 1; + for ii = 1:ngain gain = gvec(ii); - rlpol(1:nroots,ii) = vec(sortcom(roots(den + gain*num))); + rlpol(1:nroots,ii) = vec(sortcom (roots (den + gain*num))); endfor ## set smoothing tolerance - smtolx = 0.01*( max(max(real(rlpol))) - min(min(real(rlpol)))); - smtoly = 0.01*( max(max(imag(rlpol))) - min(min(imag(rlpol)))); - smtol = max(smtolx, smtoly); - rlpol = sort_roots(rlpol,smtolx, smtoly); % sort according to nearest-neighbor + smtolx = 0.01*(max (max (real (rlpol))) - min (min (real (rlpol)))); + smtoly = 0.01*(max (max (imag (rlpol))) - min (min (imag (rlpol)))); + smtol = max (smtolx, smtoly); + ## sort according to nearest-neighbor + rlpol = sort_roots (rlpol, smtolx, smtoly); - done=(nargin == 4); # perform a smoothness check - while((!done) & ngain < 1000) + done = (nargin == 4); # perform a smoothness check + while (! done && ngain < 1000) done = 1 ; # assume done - dp = abs(diff(rlpol'))'; - maxdp = max(dp); + dp = abs (diff (rlpol'))'; + maxdp = max (dp); ## search for poles whose neighbors are distant - if(lden == 2) - idx = find(dp > smtol); + if (lden == 2) + idx = find (dp > smtol); else - idx = find(maxdp > smtol); + idx = find (maxdp > smtol); endif - for ii=1:length(idx) + for ii = 1:length(idx) i1 = idx(ii); g1 = gvec(i1); p1 = rlpol(:,i1); @@ -190,39 +191,40 @@ p2 = rlpol(:,i2); ## isolate poles in p1, p2 - if( max(abs(p2-p1)) > smtol) - newg = linspace(g1,g2,5); + if (max (abs (p2-p1)) > smtol) + newg = linspace (g1, g2, 5); newg = newg(2:4); - gvec = [gvec,newg]; + gvec = [gvec,newg]; done = 0; # need to process new gains endif endfor ## process new gain values - ngain1 = length(gvec); - for ii=(ngain+1):ngain1 + ngain1 = length (gvec); + for ii = (ngain+1):ngain1 gain = gvec(ii); - rlpol(1:nroots,ii) = vec(sortcom(roots(den + gain*num))); + rlpol(1:nroots,ii) = vec(sortcom (roots (den + gain*num))); endfor - [gvec,idx] = sort(gvec); + [gvec, idx] = sort (gvec); rlpol = rlpol(:,idx); - ngain = length(gvec); - rlpol = sort_roots(rlpol,smtolx, smtoly); % sort according to nearest-neighbor + ngain = length (gvec); + ## sort according to nearest-neighbor + rlpol = sort_roots (rlpol, smtolx, smtoly); endwhile rldata = rlpol; ## Plot the data - if(nargout == 0) + if (nargout == 0) rlpolv = vec(rlpol); - axdata = [real(rlpolv),imag(rlpolv); real(olzer), imag(olzer)]; - axlim = axis2dlim(axdata); + axdata = [real(rlpolv), imag(rlpolv); real(olzer), imag(olzer)]; + axlim = axis2dlim (axdata); rldata = [real(rlpolv), imag(rlpolv) ]; - [stn,inname,outname] = sysgetsignals(sys); + [stn, inname, outname] = sysgetsignals (sys); ## build plot command args pole by pole - n_rlpol = rows(rlpol); + n_rlpol = rows (rlpol); nelts = n_rlpol+1; if (! isempty (rlzer)) nelts++; @@ -236,7 +238,7 @@ kk = 0; # asymptotes first if (n_A > 0) - len_A = 2*max(abs(axlim)); + len_A = 2*max (abs (axlim)); sigma_A = (sum(olpol) - sum(olzer))/n_A; for i_A=0:n_A-1 phi_A = pi*(2*i_A + 1)/n_A; @@ -250,7 +252,7 @@ endfor endif # locus next - for ii=1:rows(rlpol) + for ii = 1:rows(rlpol) args{1,++kk} = real (rlpol (ii,:)); args{2,kk} = imag (rlpol (ii,:)); if (ii == 1) @@ -260,22 +262,22 @@ endif endfor # poles and zeros last - args{1,++kk} = real(olpol); - args{2,kk} = imag(olpol); + args{1,++kk} = real (olpol); + args{2,kk} = imag (olpol); args{3,kk} = "rx;open loop poles;"; - if (! isempty(rlzer)) - args{1,++kk} = real(rlzer); - args{2,kk} = imag(rlzer); + if (! isempty (rlzer)) + args{1,++kk} = real (rlzer); + args{2,kk} = imag (rlzer); args{3,kk} = "go;zeros;"; endif set (gcf,"visible","off"); hplt = plot (args{:}); set (hplt(kk--), "markersize", 2); - if (! isempty(rlzer)) - set(hplt(kk--), "markersize", 2); + if (! isempty (rlzer)) + set (hplt(kk--), "markersize", 2); endif - for ii=1:rows(rlpol) + for ii = 1:rows(rlpol) set (hplt(kk--), "linewidth", 2); endfor legend ("boxon", 2); @@ -284,42 +286,42 @@ xlabel (sprintf ("Root locus from %s to %s, gain=[%f,%f]: Real axis", inname{1}, outname{1}, gvec(1), gvec(ngain))); ylabel ("Imag. axis"); - set (gcf,"visible","on"); + set (gcf (), "visible","on"); rldata = []; endif endfunction function rlpol = sort_roots (rlpol,tolx, toly) # no point sorting of you've only got one pole! - if(rows(rlpol) == 1) - return + if (rows (rlpol) == 1) + return; endif # reorder entries in each column of rlpol to be by their nearest-neighbors - dp = diff(rlpol')'; - drp = max(real(dp)); - dip = max(imag(dp)); - idx = find( drp > tolx | dip > toly ); - if(isempty(idx) ) - return + dp = diff (rlpol')'; + drp = max (real (dp)); + dip = max (imag (dp)); + idx = find (drp > tolx | dip > toly); + if (isempty (idx)) + return; endif - [np,ng] = size(rlpol); # num poles, num gains + [np, ng] = size (rlpol); # num poles, num gains for jj = idx vals = rlpol(:,[jj,jj+1]); jdx = (jj+1):ng; - for ii=1:rows(rlpol-1) + for ii = 1:rows(rlpol-1) rdx = ii:np; - dval = abs(rlpol(rdx,jj+1)-rlpol(ii,jj)); - mindist = min(dval); - sidx = min( find ( dval == mindist)) + ii - 1; - if( sidx != ii) - c1 = norm(diff(vals')); - [vals(ii,2), vals(sidx,2)] = swap( vals(ii,2), vals(sidx,2)); - c2 = norm(diff(vals')); - if(c1 > c2 ) - # perform the swap - [rlpol(ii,jdx), rlpol(sidx,jdx)] = swap( rlpol(ii,jdx), rlpol(sidx,jdx)); + dval = abs (rlpol(rdx,jj+1)-rlpol(ii,jj)); + mindist = min (dval); + sidx = min (find (dval == mindist)) + ii - 1; + if (sidx != ii) + c1 = norm (diff(vals')); + [vals(ii,2), vals(sidx,2)] = swap (vals(ii,2), vals(sidx,2)); + c2 = norm (diff (vals')); + if (c1 > c2) + ## perform the swap + [rlpol(ii,jdx), rlpol(sidx,jdx)] = swap (rlpol(ii,jdx), rlpol(sidx,jdx)); vals = rlpol(:,[jj,jj+1]); endif endif
--- a/scripts/control/base/step.m +++ b/scripts/control/base/step.m @@ -57,11 +57,7 @@ function [y, t] = step (sys, inp, tstop, n) - if ((nargin < 1) || (nargin > 4)) - print_usage (); - endif - - if (nargout > 2) + if (nargin < 1 || nargin > 4 || nargout > 2) print_usage (); endif
--- a/scripts/control/base/tzero.m +++ b/scripts/control/base/tzero.m @@ -70,23 +70,23 @@ function [zer, gain] = tzero (A, B, C, D) ## get A,B,C,D and Asys variables, regardless of initial form - if(nargin == 4) - Asys = ss(A,B,C,D); - elseif( (nargin == 1) && (! isstruct(A))) - print_usage (); - elseif(nargin != 1) + if (nargin == 4) + Asys = ss (A, B, C, D); + elseif (nargin == 1 && ! isstruct (A)) + error ("tzero: expecting argument to be system structure"); + elseif (nargin != 1) print_usage (); else Asys = A; - [A,B,C,D] = sys2ss(Asys); + [A, B, C, D] = sys2ss (Asys); endif Ao = Asys; # save for leading coefficient - siso = is_siso(Asys); - digital = is_digital(Asys); # check if it's mixed or not + siso = is_siso (Asys); + digital = is_digital (Asys); # check if it's mixed or not ## see if it's a gain block - if(isempty(A)) + if (isempty (A)) zer = []; gain = D; return; @@ -95,44 +95,50 @@ ## First, balance the system via the zero computation generalized eigenvalue ## problem balancing method (Hodel and Tiller, Linear Alg. Appl., 1992) - Asys = __zgpbal__ (Asys); [A,B,C,D] = sys2ss(Asys); # balance coefficients + ## balance coefficients + Asys = __zgpbal__ (Asys); + [A, B, C, D] = sys2ss (Asys); meps = 2*eps*norm ([A, B; C, D], "fro"); - Asys = zgreduce(Asys,meps); [A, B, C, D] = sys2ss(Asys); # ENVD algorithm - if(!isempty(A)) + ## ENVD algorithm + Asys = zgreduce (Asys, meps); + [A, B, C, D] = sys2ss (Asys); + if (! isempty (A)) ## repeat with dual system - Asys = ss(A', C', B', D'); Asys = zgreduce(Asys,meps); + Asys = ss (A', C', B', D'); + Asys = zgreduce (Asys, meps); ## transform back - [A,B,C,D] = sys2ss(Asys); Asys = ss(A', C', B', D'); + [A, B, C, D] = sys2ss (Asys); + Asys = ss (A', C', B', D'); endif zer = []; # assume none - [A,B,C,D] = sys2ss(Asys); - if( !isempty(C) ) - [W,r,Pi] = qr([C, D]'); - [nonz,ztmp] = zgrownorm(r,meps); - if(nonz) + [A, B, C, D] = sys2ss (Asys); + if (! isempty (C)) + [W, r, Pi] = qr ([C, D]'); + [nonz, ztmp] = zgrownorm (r, meps); + if (nonz) ## We can now solve the generalized eigenvalue problem. - [pp,mm] = size(D); - nn = rows(A); + [pp, mm] = size (D); + nn = rows (A); Afm = [A , B ; C, D] * W'; Bfm = [eye(nn), zeros(nn,mm); zeros(pp,nn+mm)]*W'; jdx = (mm+1):(mm+nn); Af = Afm(1:nn,jdx); Bf = Bfm(1:nn,jdx); - zer = qz(Af,Bf); + zer = qz (Af, Bf); endif endif - mz = length(zer); - [A,B,C,D] = sys2ss(Ao); # recover original system + mz = length (zer); + [A, B, C, D] = sys2ss (Ao); # recover original system ## compute leading coefficient - if ( (nargout == 2) && siso) - n = rows(A); - if ( mz == n) + if (nargout == 2 && siso) + n = rows (A); + if (mz == n) gain = D; - elseif ( mz < n ) + elseif (mz < n) gain = C*(A^(n-1-mz))*B; endif else