Mercurial > hg > octave-max
changeset 3432:e39d90787668
[project @ 2000-01-14 04:22:59 by jwe]
| author | jwe |
|---|---|
| date | Fri, 14 Jan 2000 04:28:06 +0000 |
| parents | 99ab64f4a09d |
| children | 0f6b72798944 |
| files | scripts/control/base/ctrb.m scripts/control/base/damp.m scripts/control/base/dcgain.m scripts/control/base/lsim.m scripts/control/base/ltifr.m scripts/control/base/place.m scripts/control/base/pzmap.m scripts/control/base/starp.m scripts/control/base/tzero.m scripts/control/base/tzero2.m scripts/control/ctrb.m scripts/control/damp.m scripts/control/dcgain.m scripts/control/dhinfdemo.m scripts/control/hinf/dhinfdemo.m scripts/control/hinf/wgt1o.m scripts/control/jet707.m scripts/control/lsim.m scripts/control/ltifr.m scripts/control/outlist.m scripts/control/place.m scripts/control/pzmap.m scripts/control/starp.m scripts/control/system/jet707.m scripts/control/system/ugain.m scripts/control/tzero.m scripts/control/tzero2.m scripts/control/ugain.m scripts/control/util/outlist.m scripts/control/wgt1o.m |
| diffstat | 30 files changed, 1307 insertions(+), 1307 deletions(-) [+] |
line wrap: on
line diff
new file mode 100644 --- /dev/null +++ b/scripts/control/base/ctrb.m @@ -0,0 +1,63 @@ +## Copyright (C) 1997 Kai P. 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, 59 Temple Place, Suite 330, Boston, MA 02111 USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{Qs} =} ctrb(@var{sys} @{, @var{b}@}) +## @deftypefnx {Function File} {@var{Qs} =} ctrb(@var{A}, @var{B}) +## Build controllability matrix +## @example +## 2 n-1 +## Qs = [ B AB A B ... A B ] +## @end example +## +## of a system data structure or the pair (@var{A}, @var{B}). +## +## @strong{Note} @code{ctrb} forms the controllability matrix. +## The numerical properties of @code{is_controllable} +## are much better for controllability tests. +## @end deftypefn + +## Author: Kai P. Mueller <mueller@ifr.ing.tu-bs.de> +## Created: November 4, 1997 +## based on is_controllable.m of Scottedward Hodel + +function Qs = ctrb (sys, b) + + if (nargin == 2) + a = sys; + elseif (nargin == 1 && is_struct(sys)) + sysupdate(sys,"ss"); + [a,b] = sys2ss(sys); + else + usage("ctrb(sys [, b])") + endif + + if (!is_abcd(a,b)) + Qs = []; + else + ## no need to check dimensions, we trust is_abcd(). + [na, ma] = size(a); + ## using imb avoids name conflict with the "mb" function + [inb, imb] = size(b); + Qs = zeros(na, ma*imb); + for i = 1:na + Qs(:, (i-1)*imb+1:i*imb) = b; + b = a * b; + endfor + endif +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/control/base/damp.m @@ -0,0 +1,85 @@ +## Copyright (C) 1993, 1994, 1995 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, 59 Temple Place, Suite 330, Boston, MA 02111 USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} damp(@var{p}@{, @var{tsam}@}) +## Displays eigenvalues, natural frequencies and damping ratios +## of the eigenvalues of a matrix @var{p} or the @var{A}-matrix of a +## system @var{p}, respectively. +## If @var{p} is a system, @var{tsam} must not be specified. +## If @var{p} is a matrix and @var{tsam} is specified, eigenvalues +## of @var{p} are assumed to be in @var{z}-domain. +## @end deftypefn +## @seealso{eig} + +## Author: Kai P. Mueller <mueller@ifr.ing.tu-bs.de> +## Created: September 29, 1997. + +function damp (p, tsam) + + ## assume a continuous system + DIGITAL = 0; + if(nargin < 1 || nargin > 2) + usage("damp(p,[ tsamp])") + endif + if(is_struct(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); + else + aa = p; + if (nargin == 2) + DIGITAL = 1; + t_samp = tsam; + endif + endif + if (!is_square(aa)) + error("damp: Matrix p is not square.") + endif + if (DIGITAL && t_samp <= 0.0) + error("damp: Sampling time tsam must not be <= 0.") + endif + + ## all checks done. + e = eig(aa); + [n, m] = size(aa); + if (DIGITAL) + printf(" (discrete system with sampling time %f)\n", t_samp); + endif + 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; + 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); + else + printf("%12f %12f %12f %10f %12f\n", + real(pole), imag(pole), abs(pole), d0, f0); + endif + endfor + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/control/base/dcgain.m @@ -0,0 +1,54 @@ +## Copyright (C) 1993, 1994, 1995 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, 59 Temple Place, Suite 330, Boston, MA 02111 USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{gm} =} dcgain (@var{sys}@{, tol@}) +## Returns dc-gain matrix. If dc-gain is infinite +## an empty matrix is returned. +## The argument @var{tol} is an optional tolerance for the condition +## number of @var{A}-Matrix in @var{sys} (default @var{tol} = 1.0e-10) +## @end deftypefn + +## Author: Kai P. Mueller <mueller@ifr.ing.tu-bs.de> +## Created: October 1, 1997 + +function gm = dcgain (sys, tol) + + if((nargin < 1) || (nargin > 2) || (nargout > 1)) + usage("[gm, ok] = dcgain(sys[, tol])"); + endif + if(!is_struct(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)) + 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); + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/control/base/lsim.m @@ -0,0 +1,99 @@ +## 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} {} lsim (@var{sys}, @var{u}, @var{t}@{,@var{x0}@}) +## Produce output for a linear simulation of a system +## +## Produces a plot for the output of the system, sys. +## +## U is an array that contains the system's inputs. Each column in u +## corresponds to a different time step. Each row in u corresponds to a +## different input. T is an array that contains the time index of the +## system. T should be regularly spaced. If initial conditions are required +## on the system, the x0 vector should be added to the argument list. +## +## When the lsim function is invoked with output parameters: +## [y,x] = lsim(sys,u,t,[x0]) +## a plot is not displayed, however, the data is returned in y = system output +## and x = system states. +## @end deftypefn + +## Author: David Clem +## Author: A. S. Hodel <a.s.hodel@eng.auburn.edu> +## Created: July 1995 +## modified by John Ingram for system format August 1996 + +function [y, x] = lsim (sys, u, t, x0) + + if((nargin < 3)||(nargin > 4)) + usage("[y,x] = lsim(sys,u,t[,x0])"); + endif + + if(!is_struct(sys)) + error("sys must be in system data structure"); + endif + + sys = sysupdate(sys,"ss"); + + [ncstates, ndstates, nin, nout] = sysdimensions(sys); + [a,b,c,d] = sys2ss(sys); + + 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"); + endif + 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"); + endif + if(rows(x0) > rows(a)) + error("lsim: Initial condition vector is too large"); + endif + + Ts = 0; + t(2)-t(1); + u=u'; + n = max(size(t)); + for ii = 1:(n-1) + + ## check if step size changed + if (t(ii+1) - t(ii) != Ts) + Ts = t(ii+1) - t(ii); + ## [F,G] = c2d(a,b,Ts); + dsys = c2d(sys, Ts); + [F,G] = sys2ss(dsys); + endif + + x(:,ii) = x0; + x0 = F*x0 + G*u(:,ii); + endfor + + ## pick up last point + x(:,n) = x0; + + y = c*x + d*u; + if(nargout == 0) + plot(t,y); + y=[]; + x=[]; + endif +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/control/base/ltifr.m @@ -0,0 +1,99 @@ +## 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{out} =} ltifr (@var{A}, @var{B}, @var{w}) +## @deftypefnx {Function File} {@var{out} =} ltifr (@var{sys}, @var{w}) +## Linear time invariant frequency response of single input systems +## @strong{Inputs} +## @table @var +## @item A +## @itemx B +## coefficient matrices of @math{dx/dt = A x + B u} +## @item sys +## system data structure +## @item w +## vector of frequencies +## @end table +## @strong{Outputs} +## @var{out} +## @example +## -1 +## G(s) = (jw I-A) B +## @end example +## for complex frequencies @math{s = jw}. +## @end deftypefn + +## Author: R. Bruce Tenison <btenison@eng.auburn.edu> +## Author: David Clem +## Author: A. S. Hodel <a.s.hodel@eng.auburn.edu> +## Created: July 1995 +## updated by John Ingram August 1996 for system format + +function out = ltifr (a, b, w) + + if ((nargin < 2) || (nargin > 3)) + error("incorrect number of input arguments"); + endif + + if (nargin == 2) + sys = a; + w = b; + if(!is_struct(sys)) + error("two arguments: 1st must be a system data structure"); + endif + + if (!is_vector(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 + + [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)); + endif + + if(columns(b) != 1) + error("ltifr: b(%dx%d) must be a single column vector", ... + rows(b),columns(b)); + endif + + if (!is_square(a)) + error("ltifr: A(%dx$d) must be square.",rows(a),columns(a)) + endif + + endif + + if (!is_vector(w)) + error("w must be a vector"); + endif + + ey = eye(size(a)); + lw = length(w); + out = ones(columns(a),lw); + + for ii=1:lw, + out(:,ii) = (w(ii)*ey-a)\b; + endfor +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/control/base/place.m @@ -0,0 +1,126 @@ +## Copyright (C) 1997 Jose Daniel Munoz Frias +## +## 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{K} =} place (@var{sys}, @var{P}) +## Computes the matrix K such that if the state +## is feedback with gain K, then the eigenvalues of the closed loop +## system (i.e. A-BK) are those specified in the vector P. +## +## Version: Beta (May-1997): If you have any comments, please let me know. +## (see the file place.m for my address) +## @end deftypefn + +## Author: Jose Daniel Munoz Frias + +## Universidad Pontificia Comillas +## ICAIdea +## Alberto Aguilera, 23 +## 28015 Madrid, Spain +## +## E-Mail: daniel@dea.icai.upco.es +## +## Phone: 34-1-5422800 Fax: 34-1-5596569 +## +## Algorithm taken from "The Control Handbook", IEEE press pp. 209-212 +## +## code adaped by A.S.Hodel (a.s.hodel@eng.auburn.edu) for use in controls +## toolbox + +function K = place (sys, P) + + sav_val = empty_list_elements_ok; + empty_list_elements_ok = 1; + + ## check arguments + + if(!is_struct(sys)) + error("sys must be in system data structure format (see ss2sys)"); + endif + sys = sysupdate(sys,"ss"); # make sure it has state space form up to date + if(!is_controllable(sys)) + error("sys is not controllable."); + elseif( min(size(P)) != 1) + error("P must be a vector") + else + P = reshape(P,length(P),1); # make P a column vector + endif + ## system must be purely continuous or discrete + is_digital(sys); + [n,nz,m,p] = sysdimensions(sys); + nx = n+nz; # already checked that it's not a mixed system. + if(m != 1) + error(["sys has ", num2str(m)," inputs; need only 1"]); + endif + + ## takes the A and B matrix from the system representation + [A,B]=sys2ss(sys); + sp = length(P); + if(nx == 0) + error("place: A matrix is empty (0x0)"); + elseif(nx != length(P)) + error(["A=(",num2str(nx),"x",num2str(nx),", P has ", num2str(length(P)), ... + "entries."]) + endif + + ## arguments appear to be compatible; let's give it a try! + ## The second step is the calculation of the characteristic polynomial ofA + PC=poly(A); + + ## Third step: Calculate the transformation matrix T that transforms the state + ## equation in the controllable canonical form. + + ## first we must calculate the controllability matrix M: + M=B; + AA=A; + for n = 2:nx + M(:,n)=AA*B; + AA=AA*A; + endfor + + ## second, construct the matrix W + PCO=PC(nx:-1:1); + PC1=PCO; # Matrix to shift and create W row by row + + for n = 1:nx + W(n,:) = PC1; + PC1=[PCO(n+1:nx),zeros(1,n)]; + endfor + + T=M*W; + + ## finaly the matrix K is calculated + PD = poly(P); # The desired characteristic polynomial + PD = PD(nx+1:-1:2); + PC = PC(nx+1:-1:2); + + K = (PD-PC)/T; + + ## Check if the eigenvalues of (A-BK) are the same specified in P + Pcalc = eig(A-B*K); + + Pcalc = sortcom(Pcalc); + P = sortcom(P); + + if(max( (abs(Pcalc)-abs(P))./abs(P) ) > 0.1) + disp("Place: Pole placed at more than 10% relative error from specified"); + endif + + empty_list_elements_ok = sav_val; +endfunction +
new file mode 100644 --- /dev/null +++ b/scripts/control/base/pzmap.m @@ -0,0 +1,86 @@ +## 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{zer}, @var{pol}]=} pzmap (@var{sys}) +## Plots the zeros and poles of a system in the complex plane. +## @strong{Inputs} +## @var{sys} system data structure +## +## @strong{Outputs} +## if omitted, the poles and zeros are plotted on the screen. +## otherwise, pol, zer are returned as the system poles and zeros. +## (see sys2zp for a preferable function call) +## @end deftypefn + +function [zer, pol]=pzmap (sys) + + save_emp = empty_list_elements_ok; + + empty_list_elements_ok = 1; + + if(nargin != 1) + usage("pzmap(sys) or [zer,pol] = pzmap(sys)"); + elseif (!is_struct(sys)); + error("sys must be in system format"); + endif + + [zer,pol] = sys2zp(sys); + + ## force to column vectors, split into real, imaginary parts + zerdata = poldata = []; + if(length(zer)) + zer = reshape(zer,length(zer),1); + zerdata = [real(zer(:,1)), imag(zer(:,1))]; + endif + if(length(pol)) + pol = reshape(pol,length(pol),1); + poldata = [real(pol(:,1)), imag(pol(:,1))]; + endif + + ## determine continuous or discrete plane + vars = "sz"; + varstr = vars(is_digital(sys) + 1); + + ## Plot the data + gset nologscale xy; + if(is_siso(sys)) + title(sprintf("Pole-zero map from %s to %s", ... + sysgetsignals(sys,"in",1,1), sysgetsignals(sys,"out",1,1) )); + endif + xlabel(["Re(",varstr,")"]); + ylabel(["Im(",varstr,")"]); + grid; + + ## compute axis limits + axis(axis2dlim([zerdata;poldata])); + grid + ## finally, plot the data + if(length(zer) == 0) + plot(poldata(:,1), poldata(:,2),"@12 ;poles (no zeros);"); + elseif(length(pol) == 0) + plot(zerdata(:,1), zerdata(:,2),"@31 ;zeros (no poles);"); + else + plot(zerdata(:,1), zerdata(:,2),"@31 ;zeros;", ... + poldata(:,1), poldata(:,2),"@12 ;poles;"); + endif + replot + + empty_list_elements_ok = save_emp; + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/control/base/starp.m @@ -0,0 +1,126 @@ +## 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{outputs} =} starp (@var{inputs}) +## @format +## +## sys = starp(P, K, ny, nu) +## +## Redheffer star product or upper/lower LFT, respectively. +## +## +## +-------+ +## --------->| |---------> +## | P | +## +--->| |---+ ny +## | +-------+ | +## +-------------------+ +## | | +## +----------------+ | +## | | +## | +-------+ | +## +--->| |------+ nu +## | K | +## --------->| |---------> +## +-------+ +## +## If ny and nu "consume" all inputs and outputs of K then the result +## is a lower fractional transformation. If ny and nu "consume" all +## inputs and outputs of P then the result is an upper fractional +## transformation. +## +## ny and/or nu may be negative (= negative feedback) +## @end format +## @end deftypefn + +## Author: Kai P. Mueller <mueller@ifr.ing.tu-bs.de> +## Created: May 1998 + +function sys = starp (P, K, ny, nu); + + if((nargin != 2) && (nargin != 4)) + usage("sys = starp(P, K, ny, nu)"); + endif + if (!is_struct(P)) + error("---> P must be in system data structure"); + endif + if (!is_struct(K)) + error("---> K must be in system data structure"); + endif + + P = sysupdate(P, "ss"); + [n, nz, mp, pp] = sysdimensions(P); + np = n + nz; + K = sysupdate(K, "ss"); + [n, nz, mk, pk] = sysdimensions(K); + nk = n + nz; + ny_sign = 1; + nu_sign = 1; + if (nargin == 2) + ## perform a LFT of P and K (upper or lower) + ny = min([pp, mk]); + nu = min([pk, mp]); + else + if (ny < 0) + ny = -ny; + ny_sign = -1; + endif + if (nu < 0) + nu = -nu; + nu_sign = -1; + endif + endif + if (ny > pp) + error("---> P has not enough outputs."); + endif + if (nu > mp) + error("---> P has not enough inputs."); + endif + if (ny > mk) + error("---> K has not enough inputs."); + endif + if (nu > pk) + error("---> K has not enough outputs."); + endif + nwp = mp - nu; + nzp = pp - ny; + nwk = mk - ny; + nzk = pk - nu; + if ((nwp + nwk) < 1) + error("---> no inputs left for star product."); + endif + if ((nzp + nzk) < 1) + error("---> no outputs left for star product."); + endif + + ## checks done, form sys + if (nzp) Olst = [1:nzp]; endif + if (nzk) Olst = [Olst, pp+nu+1:pp+pk]; endif + if (nwp) Ilst = [1:nwp]; endif + if (nwk) Ilst = [Ilst, mp+ny+1:mp+mk]; endif + Clst = zeros(ny+nu,2); + for ii = 1:nu + Clst(ii,:) = [nwp+ii, nu_sign*(pp+ii)]; + endfor + for ii = 1:ny + Clst(nu+ii,:) = [mp+ii, ny_sign*(nzp+ii)]; + endfor + sys = buildssic(Clst,[],Olst,Ilst,P,K); + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/control/base/tzero.m @@ -0,0 +1,125 @@ +## 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} {} tzero (@var{a}, @var{b}, @var{c}, @var{d}@{, @var{opt}@}) +## @deftypefnx {Function File} {} tzero (@var{sys}@{,@var{opt}@}) +## Compute transmission zeros of a continuous +## @example +## . +## x = Ax + Bu +## y = Cx + Du +## @end example +## or discrete +## @example +## x(k+1) = A x(k) + B u(k) +## y(k) = C x(k) + D u(k) +## @end example +## system. +## @strong{Outputs} +## @table @var +## @item zer +## transmission zeros of the system +## @item gain +## leading coefficient (pole-zero form) of SISO transfer function +## returns gain=0 if system is multivariable +## @end table +## @strong{References} +## @enumerate +## @item Emami-Naeini and Van Dooren, Automatica, 1982. +## @item Hodel, "Computation of Zeros with Balancing," 1992 Lin. Alg. Appl. +## @end enumerate +## @end deftypefn + +## Author: R. Bruce Tenison <btenison@eng.auburn.edu> +## Created: July 4, 1994 +## A. S. Hodel Aug 1995: allow for MIMO and system data structures + +function [zer, gain] = tzero (A, B, C, D) + + ## get A,B,C,D and Asys variables, regardless of initial form + if(nargin == 4) + Asys = ss2sys(A,B,C,D); + elseif( (nargin == 1) && (! is_struct(A))) + usage("[zer,gain] = tzero(A,B,C,D) or zer = tzero(Asys)"); + elseif(nargin != 1) + usage("[zer,gain] = tzero(A,B,C,D) or zer = tzero(Asys)"); + else + Asys = A; + [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 + + ## see if it's a gain block + if(isempty(A)) + zer = []; + gain = D; + return; + endif + + ## 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 + meps = 2*eps*norm ([A, B; C, D], "fro"); + Asys = zgreduce(Asys,meps); [A, B, C, D] = sys2ss(Asys); # ENVD algorithm + if(!isempty(A)) + ## repeat with dual system + Asys = ss2sys(A', C', B', D'); Asys = zgreduce(Asys,meps); + + ## transform back + [A,B,C,D] = sys2ss(Asys); Asys = ss2sys(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) + ## We can now solve the generalized eigenvalue problem. + [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); + endif + endif + + 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) + gain = D; + elseif ( mz < n ) + gain = C*(A^(n-1-mz))*B; + endif + else + gain = []; + endif +endfunction +
new file mode 100644 --- /dev/null +++ b/scripts/control/base/tzero2.m @@ -0,0 +1,66 @@ +## Copyright (C) 1993 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{zr} =} tzero2 (@var{a}, @var{b}, @var{c}, @var{d}, @var{bal}) +## Compute the transmission zeros of a, b, c, d. +## +## bal = balancing option (see balance); default is "B". +## +## Needs to incorporate @code{mvzero} algorithm to isolate finite zeros; use +## @code{tzero} instead. +## @end deftypefn + +## Author: A. S. Hodel <a.s.hodel@eng.auburn.edu> +## Created: August 1993 + +function zr = tzero2 (a, b, c, d, bal) + + if (nargin == 4) + bal = "B"; + elseif (nargin != 5) + error ("tzero: invalid number of arguments"); + endif + + [n, m, p] = abcddim (a, b, c, d); + + if (n > 0 && m > 0 && p > 0) + if (m != p) + fprintf (stderr, "tzero: number of inputs,outputs differ. squaring up"); + if (p > m) + fprintf (stderr, " by padding b and d with zeros."); + b = [b, (zeros (n, p-m))]; + d = [d, (zeros (p, p-m))]; + m = p; + else + fprintf (stderr, " by padding c and d with zeros."); + c = [c; (zeros (m-p, n))]; + d = [d; (zeros (m-p, m))]; + p = m; + endif + fprintf (stderr, "This is a kludge. Try again with SISO system."); + endif + ab = [-a, -b; c, d]; + bb = [(eye (n)), (zeros (n, m)); (zeros (p, n)), (zeros (p, m))]; + [ab,bb] = balance (ab, bb); + zr = -qz (ab, bb); + else + error ("tzero: a, b, c, d not compatible. exiting"); + endif + +endfunction
deleted file mode 100644 --- a/scripts/control/ctrb.m +++ /dev/null @@ -1,63 +0,0 @@ -## Copyright (C) 1997 Kai P. 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, 59 Temple Place, Suite 330, Boston, MA 02111 USA. - -## -*- texinfo -*- -## @deftypefn {Function File} {@var{Qs} =} ctrb(@var{sys} @{, @var{b}@}) -## @deftypefnx {Function File} {@var{Qs} =} ctrb(@var{A}, @var{B}) -## Build controllability matrix -## @example -## 2 n-1 -## Qs = [ B AB A B ... A B ] -## @end example -## -## of a system data structure or the pair (@var{A}, @var{B}). -## -## @strong{Note} @code{ctrb} forms the controllability matrix. -## The numerical properties of @code{is_controllable} -## are much better for controllability tests. -## @end deftypefn - -## Author: Kai P. Mueller <mueller@ifr.ing.tu-bs.de> -## Created: November 4, 1997 -## based on is_controllable.m of Scottedward Hodel - -function Qs = ctrb (sys, b) - - if (nargin == 2) - a = sys; - elseif (nargin == 1 && is_struct(sys)) - sysupdate(sys,"ss"); - [a,b] = sys2ss(sys); - else - usage("ctrb(sys [, b])") - endif - - if (!is_abcd(a,b)) - Qs = []; - else - ## no need to check dimensions, we trust is_abcd(). - [na, ma] = size(a); - ## using imb avoids name conflict with the "mb" function - [inb, imb] = size(b); - Qs = zeros(na, ma*imb); - for i = 1:na - Qs(:, (i-1)*imb+1:i*imb) = b; - b = a * b; - endfor - endif -endfunction
deleted file mode 100644 --- a/scripts/control/damp.m +++ /dev/null @@ -1,85 +0,0 @@ -## Copyright (C) 1993, 1994, 1995 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, 59 Temple Place, Suite 330, Boston, MA 02111 USA. - -## -*- texinfo -*- -## @deftypefn {Function File} {} damp(@var{p}@{, @var{tsam}@}) -## Displays eigenvalues, natural frequencies and damping ratios -## of the eigenvalues of a matrix @var{p} or the @var{A}-matrix of a -## system @var{p}, respectively. -## If @var{p} is a system, @var{tsam} must not be specified. -## If @var{p} is a matrix and @var{tsam} is specified, eigenvalues -## of @var{p} are assumed to be in @var{z}-domain. -## @end deftypefn -## @seealso{eig} - -## Author: Kai P. Mueller <mueller@ifr.ing.tu-bs.de> -## Created: September 29, 1997. - -function damp (p, tsam) - - ## assume a continuous system - DIGITAL = 0; - if(nargin < 1 || nargin > 2) - usage("damp(p,[ tsamp])") - endif - if(is_struct(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); - else - aa = p; - if (nargin == 2) - DIGITAL = 1; - t_samp = tsam; - endif - endif - if (!is_square(aa)) - error("damp: Matrix p is not square.") - endif - if (DIGITAL && t_samp <= 0.0) - error("damp: Sampling time tsam must not be <= 0.") - endif - - ## all checks done. - e = eig(aa); - [n, m] = size(aa); - if (DIGITAL) - printf(" (discrete system with sampling time %f)\n", t_samp); - endif - 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; - 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); - else - printf("%12f %12f %12f %10f %12f\n", - real(pole), imag(pole), abs(pole), d0, f0); - endif - endfor - -endfunction
deleted file mode 100644 --- a/scripts/control/dcgain.m +++ /dev/null @@ -1,54 +0,0 @@ -## Copyright (C) 1993, 1994, 1995 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, 59 Temple Place, Suite 330, Boston, MA 02111 USA. - -## -*- texinfo -*- -## @deftypefn {Function File} {@var{gm} =} dcgain (@var{sys}@{, tol@}) -## Returns dc-gain matrix. If dc-gain is infinite -## an empty matrix is returned. -## The argument @var{tol} is an optional tolerance for the condition -## number of @var{A}-Matrix in @var{sys} (default @var{tol} = 1.0e-10) -## @end deftypefn - -## Author: Kai P. Mueller <mueller@ifr.ing.tu-bs.de> -## Created: October 1, 1997 - -function gm = dcgain (sys, tol) - - if((nargin < 1) || (nargin > 2) || (nargout > 1)) - usage("[gm, ok] = dcgain(sys[, tol])"); - endif - if(!is_struct(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)) - 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); - endif - -endfunction
deleted file mode 100644 --- a/scripts/control/dhinfdemo.m +++ /dev/null @@ -1,143 +0,0 @@ -## Copyright (C) 1996, 1998 Kai P. 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, 59 Temple Place, Suite 330, Boston, MA 02111 USA. - -## ------------------------------------------------------------ -## dhinfdemo Design of a discrete H_infinity controller. -## This is not a true discrete design. The design -## is carried out in continuous time while the -## effect of sampling is described by a bilinear -## transformation of the sampled system. -## This method works quite well if the sampling -## period is "small" compared to the plant time -## constants. -## -## This is a script file for OCTAVE. -## ------------------------------------------------------------ -## -## continuous plant: -## 1 -## G(s) = -------------- -## (s + 2)(s + 1) -## -## discretised plant with ZOH (Sampling period = Ts = 1 second) -## -## 0.39958z + 0.14700 -## G(s) = -------------------------- -## (z - 0.36788)(z - 0.13533) -## -## +----+ -## -------------------->| W1 |---> v1 -## z | +----+ -## ----|-------------+ || T || => min. -## | | vz infty -## | +---+ v +----+ -## *--->| G |--->O--*-->| W2 |---> v2 -## | +---+ | +----+ -## | | -## | +---+ | -## -----| K |<------- -## +---+ -## -## W1 and W2 are the robustness and performancs weighting -## functions - -## K. Mueller, <mueller@ifr.ing.tu-bs.de> -## Technical University of Braunschweig, IfR - -echo off -disp(" "); -disp(" --------------------------------------------------"); -disp(" Discrete H_infinity optimal control for the plant:"); -disp(" "); -disp(" 0.39958z + 0.14700"); -disp(" G(s) = --------------------------"); -disp(" (z - 0.36788)(z - 0.13533)"); -disp(" --------------------------------------------------"); -disp(" "); - -disp("sampling time:") -cmd = "Ts = 1.0;"; -disp(cmd); -eval(cmd); -disp("weighting on actuator value u"); -cmd = "W1 = wgt1o(0.1, 200.0, 50.0);"; -disp(cmd); -eval(cmd); -disp("weighting on controlled variable y"); -cmd = "W2 = wgt1o(350.0, 0.05, 0.0002);"; -disp(cmd); -eval(cmd); -## omega axis (column vector) -ww = vec(logspace(-4.99, 3.99, 100)); - -disp("Create ZOH equivalent model of a continuous plant"); -cmd = "G = tf2sys(2,[1 3 2]); Gd = c2d(G, Ts);"; -run_cmd - -## w-plane (continuous representation of the sampled system) -disp("W-plane transform of discrete time system:"); -cmd = "Gw = d2c(Gd, \"bi\");"; -run_cmd - -disp(" "); -disp(" o building P..."); -## need One as the pseudo transfer function One = 1 -cmd = "One = ugain(1);"; -disp(cmd); -eval(cmd); -cmd = " psys = buildssic([1 4;2 4;3 1],[3],[2 3 5],[3 4],Gw,W1,W2,One);"; -run_cmd; -disp(" o controller design..."); -cmd = "[K, gfin, GWC] = hinfsyn(psys, 1, 1, 0.1, 10.0, 0.02);"; -run_cmd - -disp(" "); -fig_n = 1; -yn = input(" * Plot magnitudes of W1KS and W2S? [n]: ","S"); -if (length(yn) >= 1) - if ((yn(1) == "y") || (yn(1) == 'Y')) - disp(" o magnitudes of W1KS and W2S..."); - gwx = sysprune(GWC, 1, 1); - mag1 = bode(gwx, ww); - if (columns(mag1) > 1); mag1 = mag1'; endif - gwx = sysprune(GWC, 2, 1); - mag2 = bode(gwx, ww); - if (columns(mag2) > 1); mag2 = mag2'; endif - figure(fig_n) - fig_n = fig_n + 1; - gset grid - loglog(ww, [mag1 mag2]); - endif -endif - -Kd = c2d(K, "bi", Ts); -GG = buildssic([1 2; 2 1], [], [1 2], [-2], Gd, Kd); -disp(" o closed loop poles..."); -damp(GG); - -disp(" "); -yn = input(" * Plot closed loop step responses? [n]: ","S"); -if (length(yn) >= 1) - if ((yn(1) == "y") || (yn(1) == 'Y')) - disp(" o step responses of T and KS..."); - figure(fig_n) - step(GG, 1, 10); - endif -endif - -## --------- End of dhinfdemo/kpm
new file mode 100644 --- /dev/null +++ b/scripts/control/hinf/dhinfdemo.m @@ -0,0 +1,143 @@ +## Copyright (C) 1996, 1998 Kai P. 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, 59 Temple Place, Suite 330, Boston, MA 02111 USA. + +## ------------------------------------------------------------ +## dhinfdemo Design of a discrete H_infinity controller. +## This is not a true discrete design. The design +## is carried out in continuous time while the +## effect of sampling is described by a bilinear +## transformation of the sampled system. +## This method works quite well if the sampling +## period is "small" compared to the plant time +## constants. +## +## This is a script file for OCTAVE. +## ------------------------------------------------------------ +## +## continuous plant: +## 1 +## G(s) = -------------- +## (s + 2)(s + 1) +## +## discretised plant with ZOH (Sampling period = Ts = 1 second) +## +## 0.39958z + 0.14700 +## G(s) = -------------------------- +## (z - 0.36788)(z - 0.13533) +## +## +----+ +## -------------------->| W1 |---> v1 +## z | +----+ +## ----|-------------+ || T || => min. +## | | vz infty +## | +---+ v +----+ +## *--->| G |--->O--*-->| W2 |---> v2 +## | +---+ | +----+ +## | | +## | +---+ | +## -----| K |<------- +## +---+ +## +## W1 and W2 are the robustness and performancs weighting +## functions + +## K. Mueller, <mueller@ifr.ing.tu-bs.de> +## Technical University of Braunschweig, IfR + +echo off +disp(" "); +disp(" --------------------------------------------------"); +disp(" Discrete H_infinity optimal control for the plant:"); +disp(" "); +disp(" 0.39958z + 0.14700"); +disp(" G(s) = --------------------------"); +disp(" (z - 0.36788)(z - 0.13533)"); +disp(" --------------------------------------------------"); +disp(" "); + +disp("sampling time:") +cmd = "Ts = 1.0;"; +disp(cmd); +eval(cmd); +disp("weighting on actuator value u"); +cmd = "W1 = wgt1o(0.1, 200.0, 50.0);"; +disp(cmd); +eval(cmd); +disp("weighting on controlled variable y"); +cmd = "W2 = wgt1o(350.0, 0.05, 0.0002);"; +disp(cmd); +eval(cmd); +## omega axis (column vector) +ww = vec(logspace(-4.99, 3.99, 100)); + +disp("Create ZOH equivalent model of a continuous plant"); +cmd = "G = tf2sys(2,[1 3 2]); Gd = c2d(G, Ts);"; +run_cmd + +## w-plane (continuous representation of the sampled system) +disp("W-plane transform of discrete time system:"); +cmd = "Gw = d2c(Gd, \"bi\");"; +run_cmd + +disp(" "); +disp(" o building P..."); +## need One as the pseudo transfer function One = 1 +cmd = "One = ugain(1);"; +disp(cmd); +eval(cmd); +cmd = " psys = buildssic([1 4;2 4;3 1],[3],[2 3 5],[3 4],Gw,W1,W2,One);"; +run_cmd; +disp(" o controller design..."); +cmd = "[K, gfin, GWC] = hinfsyn(psys, 1, 1, 0.1, 10.0, 0.02);"; +run_cmd + +disp(" "); +fig_n = 1; +yn = input(" * Plot magnitudes of W1KS and W2S? [n]: ","S"); +if (length(yn) >= 1) + if ((yn(1) == "y") || (yn(1) == 'Y')) + disp(" o magnitudes of W1KS and W2S..."); + gwx = sysprune(GWC, 1, 1); + mag1 = bode(gwx, ww); + if (columns(mag1) > 1); mag1 = mag1'; endif + gwx = sysprune(GWC, 2, 1); + mag2 = bode(gwx, ww); + if (columns(mag2) > 1); mag2 = mag2'; endif + figure(fig_n) + fig_n = fig_n + 1; + gset grid + loglog(ww, [mag1 mag2]); + endif +endif + +Kd = c2d(K, "bi", Ts); +GG = buildssic([1 2; 2 1], [], [1 2], [-2], Gd, Kd); +disp(" o closed loop poles..."); +damp(GG); + +disp(" "); +yn = input(" * Plot closed loop step responses? [n]: ","S"); +if (length(yn) >= 1) + if ((yn(1) == "y") || (yn(1) == 'Y')) + disp(" o step responses of T and KS..."); + figure(fig_n) + step(GG, 1, 10); + endif +endif + +## --------- End of dhinfdemo/kpm
new file mode 100644 --- /dev/null +++ b/scripts/control/hinf/wgt1o.m @@ -0,0 +1,59 @@ +## Copyright (C) 1998 Kai P. 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, 59 Temple Place, Suite 330, Boston, MA 02111 USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{wsys} =} wgt1o (@var{vl}, @var{vh}, @var{fc}) +## State space description of a first order weighting function. +## +## Weighting function are needed by the H2/H_infinity design procedure. +## These function are part of thye augmented plant P (see hinfdemo +## for an applicattion example). +## +## vl = Gain at low frequencies +## +## vh = Gain at high frequencies +## +## fc = Corner frequency (in Hz, *not* in rad/sec) +## @end deftypefn + +## Author: Kai P. Mueller <mueller@ifr.ing.tu-bs.de> +## Created: September 30, 1997 + +function wsys = wgt1o (vl, vh, fc) + + if (nargin != 3) + usage("wsys = wgt1o(vl, vh, fc)"); + endif + + if(nargout > 1) + usage("wsys = wgt1o(vl, vh, fc)"); + endif + + if (vl == vh) + a = []; + b = []; + c = []; + else + a = [-2*pi*fc]; + b = [-2*pi*fc]; + c = [vh-vl]; + endif + d=[vh]; + + wsys = ss2sys(a,b,c,d); +endfunction
deleted file mode 100644 --- a/scripts/control/jet707.m +++ /dev/null @@ -1,57 +0,0 @@ -## Copyright (C) 1997 Kai P. 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, 59 Temple Place, Suite 330, Boston, MA 02111 USA. - -## -*- texinfo -*- -## @deftypefn {Function File} {@var{outsys} =} jet707 () -## Creates linearized state space model of a Boeing 707-321 aircraft -## at v=80m/s. (M = 0.26, Ga0 = -3 deg, alpha0 = 4 deg, kappa = 50 deg) -## System inputs: (1) thrust and (2) elevator angle -## System outputs: (1) airspeed and (2) pitch angle -## Ref: R. Brockhaus: Flugregelung (Flight Control), Springer, 1994 -## @end deftypefn -## @seealso{ord2} - -## Author: Kai P. Mueller <mueller@ifr.ing.tu-bs.de> -## Created: September 28, 1997 - -function outsys = jet707 () - - if (nargin != 0) - usage("outsys = jet707()") - endif - if (nargin > 1) - usage("outsys = jet707()") - endif - - a = [ -0.46E-01, 0.10681415316, 0.0, -0.17121680433; - -0.1675901504661613, -0.515, 1.0, 0.6420630320636088E-02; - 0.1543104215347786, -0.547945, -0.906, -0.1521689385990753E-02; - 0.0, 0.0, 1.0, 0.0 ]; - b = [ 0.1602300107479095, 0.2111848453E-02; - 0.8196877780963616E-02, -0.3025E-01; - 0.9173594317692437E-01, -0.75283075; - 0.0, 0.0 ]; - c = [ 1.0, 0.0, 0.0, 0.0; - 0.0, 0.0, 0.0, 1.0 ]; - d=zeros(2,2); - inam = ["thrust"; "rudder"]; - onam = ["speed"; "pitch"]; - snam = ["x1"; "x2"; "x3"; "x4"]; - outsys = ss2sys(a, b, c, d, 0.0, 4, 0, snam, inam, onam); - -endfunction
deleted file mode 100644 --- a/scripts/control/lsim.m +++ /dev/null @@ -1,99 +0,0 @@ -## 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} {} lsim (@var{sys}, @var{u}, @var{t}@{,@var{x0}@}) -## Produce output for a linear simulation of a system -## -## Produces a plot for the output of the system, sys. -## -## U is an array that contains the system's inputs. Each column in u -## corresponds to a different time step. Each row in u corresponds to a -## different input. T is an array that contains the time index of the -## system. T should be regularly spaced. If initial conditions are required -## on the system, the x0 vector should be added to the argument list. -## -## When the lsim function is invoked with output parameters: -## [y,x] = lsim(sys,u,t,[x0]) -## a plot is not displayed, however, the data is returned in y = system output -## and x = system states. -## @end deftypefn - -## Author: David Clem -## Author: A. S. Hodel <a.s.hodel@eng.auburn.edu> -## Created: July 1995 -## modified by John Ingram for system format August 1996 - -function [y, x] = lsim (sys, u, t, x0) - - if((nargin < 3)||(nargin > 4)) - usage("[y,x] = lsim(sys,u,t[,x0])"); - endif - - if(!is_struct(sys)) - error("sys must be in system data structure"); - endif - - sys = sysupdate(sys,"ss"); - - [ncstates, ndstates, nin, nout] = sysdimensions(sys); - [a,b,c,d] = sys2ss(sys); - - 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"); - endif - 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"); - endif - if(rows(x0) > rows(a)) - error("lsim: Initial condition vector is too large"); - endif - - Ts = 0; - t(2)-t(1); - u=u'; - n = max(size(t)); - for ii = 1:(n-1) - - ## check if step size changed - if (t(ii+1) - t(ii) != Ts) - Ts = t(ii+1) - t(ii); - ## [F,G] = c2d(a,b,Ts); - dsys = c2d(sys, Ts); - [F,G] = sys2ss(dsys); - endif - - x(:,ii) = x0; - x0 = F*x0 + G*u(:,ii); - endfor - - ## pick up last point - x(:,n) = x0; - - y = c*x + d*u; - if(nargout == 0) - plot(t,y); - y=[]; - x=[]; - endif -endfunction
deleted file mode 100644 --- a/scripts/control/ltifr.m +++ /dev/null @@ -1,99 +0,0 @@ -## 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{out} =} ltifr (@var{A}, @var{B}, @var{w}) -## @deftypefnx {Function File} {@var{out} =} ltifr (@var{sys}, @var{w}) -## Linear time invariant frequency response of single input systems -## @strong{Inputs} -## @table @var -## @item A -## @itemx B -## coefficient matrices of @math{dx/dt = A x + B u} -## @item sys -## system data structure -## @item w -## vector of frequencies -## @end table -## @strong{Outputs} -## @var{out} -## @example -## -1 -## G(s) = (jw I-A) B -## @end example -## for complex frequencies @math{s = jw}. -## @end deftypefn - -## Author: R. Bruce Tenison <btenison@eng.auburn.edu> -## Author: David Clem -## Author: A. S. Hodel <a.s.hodel@eng.auburn.edu> -## Created: July 1995 -## updated by John Ingram August 1996 for system format - -function out = ltifr (a, b, w) - - if ((nargin < 2) || (nargin > 3)) - error("incorrect number of input arguments"); - endif - - if (nargin == 2) - sys = a; - w = b; - if(!is_struct(sys)) - error("two arguments: 1st must be a system data structure"); - endif - - if (!is_vector(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 - - [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)); - endif - - if(columns(b) != 1) - error("ltifr: b(%dx%d) must be a single column vector", ... - rows(b),columns(b)); - endif - - if (!is_square(a)) - error("ltifr: A(%dx$d) must be square.",rows(a),columns(a)) - endif - - endif - - if (!is_vector(w)) - error("w must be a vector"); - endif - - ey = eye(size(a)); - lw = length(w); - out = ones(columns(a),lw); - - for ii=1:lw, - out(:,ii) = (w(ii)*ey-a)\b; - endfor -endfunction
deleted file mode 100644 --- a/scripts/control/outlist.m +++ /dev/null @@ -1,80 +0,0 @@ -## 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
deleted file mode 100644 --- a/scripts/control/place.m +++ /dev/null @@ -1,126 +0,0 @@ -## Copyright (C) 1997 Jose Daniel Munoz Frias -## -## 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{K} =} place (@var{sys}, @var{P}) -## Computes the matrix K such that if the state -## is feedback with gain K, then the eigenvalues of the closed loop -## system (i.e. A-BK) are those specified in the vector P. -## -## Version: Beta (May-1997): If you have any comments, please let me know. -## (see the file place.m for my address) -## @end deftypefn - -## Author: Jose Daniel Munoz Frias - -## Universidad Pontificia Comillas -## ICAIdea -## Alberto Aguilera, 23 -## 28015 Madrid, Spain -## -## E-Mail: daniel@dea.icai.upco.es -## -## Phone: 34-1-5422800 Fax: 34-1-5596569 -## -## Algorithm taken from "The Control Handbook", IEEE press pp. 209-212 -## -## code adaped by A.S.Hodel (a.s.hodel@eng.auburn.edu) for use in controls -## toolbox - -function K = place (sys, P) - - sav_val = empty_list_elements_ok; - empty_list_elements_ok = 1; - - ## check arguments - - if(!is_struct(sys)) - error("sys must be in system data structure format (see ss2sys)"); - endif - sys = sysupdate(sys,"ss"); # make sure it has state space form up to date - if(!is_controllable(sys)) - error("sys is not controllable."); - elseif( min(size(P)) != 1) - error("P must be a vector") - else - P = reshape(P,length(P),1); # make P a column vector - endif - ## system must be purely continuous or discrete - is_digital(sys); - [n,nz,m,p] = sysdimensions(sys); - nx = n+nz; # already checked that it's not a mixed system. - if(m != 1) - error(["sys has ", num2str(m)," inputs; need only 1"]); - endif - - ## takes the A and B matrix from the system representation - [A,B]=sys2ss(sys); - sp = length(P); - if(nx == 0) - error("place: A matrix is empty (0x0)"); - elseif(nx != length(P)) - error(["A=(",num2str(nx),"x",num2str(nx),", P has ", num2str(length(P)), ... - "entries."]) - endif - - ## arguments appear to be compatible; let's give it a try! - ## The second step is the calculation of the characteristic polynomial ofA - PC=poly(A); - - ## Third step: Calculate the transformation matrix T that transforms the state - ## equation in the controllable canonical form. - - ## first we must calculate the controllability matrix M: - M=B; - AA=A; - for n = 2:nx - M(:,n)=AA*B; - AA=AA*A; - endfor - - ## second, construct the matrix W - PCO=PC(nx:-1:1); - PC1=PCO; # Matrix to shift and create W row by row - - for n = 1:nx - W(n,:) = PC1; - PC1=[PCO(n+1:nx),zeros(1,n)]; - endfor - - T=M*W; - - ## finaly the matrix K is calculated - PD = poly(P); # The desired characteristic polynomial - PD = PD(nx+1:-1:2); - PC = PC(nx+1:-1:2); - - K = (PD-PC)/T; - - ## Check if the eigenvalues of (A-BK) are the same specified in P - Pcalc = eig(A-B*K); - - Pcalc = sortcom(Pcalc); - P = sortcom(P); - - if(max( (abs(Pcalc)-abs(P))./abs(P) ) > 0.1) - disp("Place: Pole placed at more than 10% relative error from specified"); - endif - - empty_list_elements_ok = sav_val; -endfunction -
deleted file mode 100644 --- a/scripts/control/pzmap.m +++ /dev/null @@ -1,86 +0,0 @@ -## 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{zer}, @var{pol}]=} pzmap (@var{sys}) -## Plots the zeros and poles of a system in the complex plane. -## @strong{Inputs} -## @var{sys} system data structure -## -## @strong{Outputs} -## if omitted, the poles and zeros are plotted on the screen. -## otherwise, pol, zer are returned as the system poles and zeros. -## (see sys2zp for a preferable function call) -## @end deftypefn - -function [zer, pol]=pzmap (sys) - - save_emp = empty_list_elements_ok; - - empty_list_elements_ok = 1; - - if(nargin != 1) - usage("pzmap(sys) or [zer,pol] = pzmap(sys)"); - elseif (!is_struct(sys)); - error("sys must be in system format"); - endif - - [zer,pol] = sys2zp(sys); - - ## force to column vectors, split into real, imaginary parts - zerdata = poldata = []; - if(length(zer)) - zer = reshape(zer,length(zer),1); - zerdata = [real(zer(:,1)), imag(zer(:,1))]; - endif - if(length(pol)) - pol = reshape(pol,length(pol),1); - poldata = [real(pol(:,1)), imag(pol(:,1))]; - endif - - ## determine continuous or discrete plane - vars = "sz"; - varstr = vars(is_digital(sys) + 1); - - ## Plot the data - gset nologscale xy; - if(is_siso(sys)) - title(sprintf("Pole-zero map from %s to %s", ... - sysgetsignals(sys,"in",1,1), sysgetsignals(sys,"out",1,1) )); - endif - xlabel(["Re(",varstr,")"]); - ylabel(["Im(",varstr,")"]); - grid; - - ## compute axis limits - axis(axis2dlim([zerdata;poldata])); - grid - ## finally, plot the data - if(length(zer) == 0) - plot(poldata(:,1), poldata(:,2),"@12 ;poles (no zeros);"); - elseif(length(pol) == 0) - plot(zerdata(:,1), zerdata(:,2),"@31 ;zeros (no poles);"); - else - plot(zerdata(:,1), zerdata(:,2),"@31 ;zeros;", ... - poldata(:,1), poldata(:,2),"@12 ;poles;"); - endif - replot - - empty_list_elements_ok = save_emp; - -endfunction
deleted file mode 100644 --- a/scripts/control/starp.m +++ /dev/null @@ -1,126 +0,0 @@ -## 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{outputs} =} starp (@var{inputs}) -## @format -## -## sys = starp(P, K, ny, nu) -## -## Redheffer star product or upper/lower LFT, respectively. -## -## -## +-------+ -## --------->| |---------> -## | P | -## +--->| |---+ ny -## | +-------+ | -## +-------------------+ -## | | -## +----------------+ | -## | | -## | +-------+ | -## +--->| |------+ nu -## | K | -## --------->| |---------> -## +-------+ -## -## If ny and nu "consume" all inputs and outputs of K then the result -## is a lower fractional transformation. If ny and nu "consume" all -## inputs and outputs of P then the result is an upper fractional -## transformation. -## -## ny and/or nu may be negative (= negative feedback) -## @end format -## @end deftypefn - -## Author: Kai P. Mueller <mueller@ifr.ing.tu-bs.de> -## Created: May 1998 - -function sys = starp (P, K, ny, nu); - - if((nargin != 2) && (nargin != 4)) - usage("sys = starp(P, K, ny, nu)"); - endif - if (!is_struct(P)) - error("---> P must be in system data structure"); - endif - if (!is_struct(K)) - error("---> K must be in system data structure"); - endif - - P = sysupdate(P, "ss"); - [n, nz, mp, pp] = sysdimensions(P); - np = n + nz; - K = sysupdate(K, "ss"); - [n, nz, mk, pk] = sysdimensions(K); - nk = n + nz; - ny_sign = 1; - nu_sign = 1; - if (nargin == 2) - ## perform a LFT of P and K (upper or lower) - ny = min([pp, mk]); - nu = min([pk, mp]); - else - if (ny < 0) - ny = -ny; - ny_sign = -1; - endif - if (nu < 0) - nu = -nu; - nu_sign = -1; - endif - endif - if (ny > pp) - error("---> P has not enough outputs."); - endif - if (nu > mp) - error("---> P has not enough inputs."); - endif - if (ny > mk) - error("---> K has not enough inputs."); - endif - if (nu > pk) - error("---> K has not enough outputs."); - endif - nwp = mp - nu; - nzp = pp - ny; - nwk = mk - ny; - nzk = pk - nu; - if ((nwp + nwk) < 1) - error("---> no inputs left for star product."); - endif - if ((nzp + nzk) < 1) - error("---> no outputs left for star product."); - endif - - ## checks done, form sys - if (nzp) Olst = [1:nzp]; endif - if (nzk) Olst = [Olst, pp+nu+1:pp+pk]; endif - if (nwp) Ilst = [1:nwp]; endif - if (nwk) Ilst = [Ilst, mp+ny+1:mp+mk]; endif - Clst = zeros(ny+nu,2); - for ii = 1:nu - Clst(ii,:) = [nwp+ii, nu_sign*(pp+ii)]; - endfor - for ii = 1:ny - Clst(nu+ii,:) = [mp+ii, ny_sign*(nzp+ii)]; - endfor - sys = buildssic(Clst,[],Olst,Ilst,P,K); - -endfunction
new file mode 100644 --- /dev/null +++ b/scripts/control/system/jet707.m @@ -0,0 +1,57 @@ +## Copyright (C) 1997 Kai P. 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, 59 Temple Place, Suite 330, Boston, MA 02111 USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{outsys} =} jet707 () +## Creates linearized state space model of a Boeing 707-321 aircraft +## at v=80m/s. (M = 0.26, Ga0 = -3 deg, alpha0 = 4 deg, kappa = 50 deg) +## System inputs: (1) thrust and (2) elevator angle +## System outputs: (1) airspeed and (2) pitch angle +## Ref: R. Brockhaus: Flugregelung (Flight Control), Springer, 1994 +## @end deftypefn +## @seealso{ord2} + +## Author: Kai P. Mueller <mueller@ifr.ing.tu-bs.de> +## Created: September 28, 1997 + +function outsys = jet707 () + + if (nargin != 0) + usage("outsys = jet707()") + endif + if (nargin > 1) + usage("outsys = jet707()") + endif + + a = [ -0.46E-01, 0.10681415316, 0.0, -0.17121680433; + -0.1675901504661613, -0.515, 1.0, 0.6420630320636088E-02; + 0.1543104215347786, -0.547945, -0.906, -0.1521689385990753E-02; + 0.0, 0.0, 1.0, 0.0 ]; + b = [ 0.1602300107479095, 0.2111848453E-02; + 0.8196877780963616E-02, -0.3025E-01; + 0.9173594317692437E-01, -0.75283075; + 0.0, 0.0 ]; + c = [ 1.0, 0.0, 0.0, 0.0; + 0.0, 0.0, 0.0, 1.0 ]; + d=zeros(2,2); + inam = ["thrust"; "rudder"]; + onam = ["speed"; "pitch"]; + snam = ["x1"; "x2"; "x3"; "x4"]; + outsys = ss2sys(a, b, c, d, 0.0, 4, 0, snam, inam, onam); + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/control/system/ugain.m @@ -0,0 +1,39 @@ +## Copyright (C) 1997 Kai P. 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, 59 Temple Place, Suite 330, Boston, MA 02111 USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{outsys} =} ugain (@var{n}) +## Creates a system with unity gain, no states. +## This trivial system is sometimes needed to create arbitrary +## complex systems from simple systems with buildssic. +## Watch out if you are forming sampled systems since "ugain" +## does not contain a sampling period. +## @end deftypefn +## @seealso{hinfdemo and jet707} + +## Author: Kai P. Mueller <mueller@ifr.ing.tu-bs.de> +## Created: April 1998 + +function outsys = ugain (n) + + if((nargin != 1) || (nargout > 1)) + usage("outsys = ugain(n)") + endif + outsys = ss2sys([],[],[],eye(n)); + +endfunction
deleted file mode 100644 --- a/scripts/control/tzero.m +++ /dev/null @@ -1,125 +0,0 @@ -## 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} {} tzero (@var{a}, @var{b}, @var{c}, @var{d}@{, @var{opt}@}) -## @deftypefnx {Function File} {} tzero (@var{sys}@{,@var{opt}@}) -## Compute transmission zeros of a continuous -## @example -## . -## x = Ax + Bu -## y = Cx + Du -## @end example -## or discrete -## @example -## x(k+1) = A x(k) + B u(k) -## y(k) = C x(k) + D u(k) -## @end example -## system. -## @strong{Outputs} -## @table @var -## @item zer -## transmission zeros of the system -## @item gain -## leading coefficient (pole-zero form) of SISO transfer function -## returns gain=0 if system is multivariable -## @end table -## @strong{References} -## @enumerate -## @item Emami-Naeini and Van Dooren, Automatica, 1982. -## @item Hodel, "Computation of Zeros with Balancing," 1992 Lin. Alg. Appl. -## @end enumerate -## @end deftypefn - -## Author: R. Bruce Tenison <btenison@eng.auburn.edu> -## Created: July 4, 1994 -## A. S. Hodel Aug 1995: allow for MIMO and system data structures - -function [zer, gain] = tzero (A, B, C, D) - - ## get A,B,C,D and Asys variables, regardless of initial form - if(nargin == 4) - Asys = ss2sys(A,B,C,D); - elseif( (nargin == 1) && (! is_struct(A))) - usage("[zer,gain] = tzero(A,B,C,D) or zer = tzero(Asys)"); - elseif(nargin != 1) - usage("[zer,gain] = tzero(A,B,C,D) or zer = tzero(Asys)"); - else - Asys = A; - [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 - - ## see if it's a gain block - if(isempty(A)) - zer = []; - gain = D; - return; - endif - - ## 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 - meps = 2*eps*norm ([A, B; C, D], "fro"); - Asys = zgreduce(Asys,meps); [A, B, C, D] = sys2ss(Asys); # ENVD algorithm - if(!isempty(A)) - ## repeat with dual system - Asys = ss2sys(A', C', B', D'); Asys = zgreduce(Asys,meps); - - ## transform back - [A,B,C,D] = sys2ss(Asys); Asys = ss2sys(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) - ## We can now solve the generalized eigenvalue problem. - [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); - endif - endif - - 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) - gain = D; - elseif ( mz < n ) - gain = C*(A^(n-1-mz))*B; - endif - else - gain = []; - endif -endfunction -
deleted file mode 100644 --- a/scripts/control/tzero2.m +++ /dev/null @@ -1,66 +0,0 @@ -## Copyright (C) 1993 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{zr} =} tzero2 (@var{a}, @var{b}, @var{c}, @var{d}, @var{bal}) -## Compute the transmission zeros of a, b, c, d. -## -## bal = balancing option (see balance); default is "B". -## -## Needs to incorporate @code{mvzero} algorithm to isolate finite zeros; use -## @code{tzero} instead. -## @end deftypefn - -## Author: A. S. Hodel <a.s.hodel@eng.auburn.edu> -## Created: August 1993 - -function zr = tzero2 (a, b, c, d, bal) - - if (nargin == 4) - bal = "B"; - elseif (nargin != 5) - error ("tzero: invalid number of arguments"); - endif - - [n, m, p] = abcddim (a, b, c, d); - - if (n > 0 && m > 0 && p > 0) - if (m != p) - fprintf (stderr, "tzero: number of inputs,outputs differ. squaring up"); - if (p > m) - fprintf (stderr, " by padding b and d with zeros."); - b = [b, (zeros (n, p-m))]; - d = [d, (zeros (p, p-m))]; - m = p; - else - fprintf (stderr, " by padding c and d with zeros."); - c = [c; (zeros (m-p, n))]; - d = [d; (zeros (m-p, m))]; - p = m; - endif - fprintf (stderr, "This is a kludge. Try again with SISO system."); - endif - ab = [-a, -b; c, d]; - bb = [(eye (n)), (zeros (n, m)); (zeros (p, n)), (zeros (p, m))]; - [ab,bb] = balance (ab, bb); - zr = -qz (ab, bb); - else - error ("tzero: a, b, c, d not compatible. exiting"); - endif - -endfunction
deleted file mode 100644 --- a/scripts/control/ugain.m +++ /dev/null @@ -1,39 +0,0 @@ -## Copyright (C) 1997 Kai P. 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, 59 Temple Place, Suite 330, Boston, MA 02111 USA. - -## -*- texinfo -*- -## @deftypefn {Function File} {@var{outsys} =} ugain (@var{n}) -## Creates a system with unity gain, no states. -## This trivial system is sometimes needed to create arbitrary -## complex systems from simple systems with buildssic. -## Watch out if you are forming sampled systems since "ugain" -## does not contain a sampling period. -## @end deftypefn -## @seealso{hinfdemo and jet707} - -## Author: Kai P. Mueller <mueller@ifr.ing.tu-bs.de> -## Created: April 1998 - -function outsys = ugain (n) - - if((nargin != 1) || (nargout > 1)) - usage("outsys = ugain(n)") - endif - outsys = ss2sys([],[],[],eye(n)); - -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
deleted file mode 100644 --- a/scripts/control/wgt1o.m +++ /dev/null @@ -1,59 +0,0 @@ -## Copyright (C) 1998 Kai P. 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, 59 Temple Place, Suite 330, Boston, MA 02111 USA. - -## -*- texinfo -*- -## @deftypefn {Function File} {@var{wsys} =} wgt1o (@var{vl}, @var{vh}, @var{fc}) -## State space description of a first order weighting function. -## -## Weighting function are needed by the H2/H_infinity design procedure. -## These function are part of thye augmented plant P (see hinfdemo -## for an applicattion example). -## -## vl = Gain at low frequencies -## -## vh = Gain at high frequencies -## -## fc = Corner frequency (in Hz, *not* in rad/sec) -## @end deftypefn - -## Author: Kai P. Mueller <mueller@ifr.ing.tu-bs.de> -## Created: September 30, 1997 - -function wsys = wgt1o (vl, vh, fc) - - if (nargin != 3) - usage("wsys = wgt1o(vl, vh, fc)"); - endif - - if(nargout > 1) - usage("wsys = wgt1o(vl, vh, fc)"); - endif - - if (vl == vh) - a = []; - b = []; - c = []; - else - a = [-2*pi*fc]; - b = [-2*pi*fc]; - c = [vh-vl]; - endif - d=[vh]; - - wsys = ss2sys(a,b,c,d); -endfunction