new file mode 100644
--- /dev/null
+++ b/scripts/control/system/sysmin.m
@@ -0,0 +1,177 @@
+## 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.
+
+## [retsys,nc,no] = sysmin(sys{,flg});
+## return a minimal (or reduced order) system
+## inputs:
+##   sys: system data structure
+##   flg: 0 [default] return minimal system; state names lost
+##      : 1           return system with physical states removed that
+##                    are either uncontrollable or unobservable
+##                    (cannot reduce further without discarding physical
+##                    meaning of states)
+## outputs:
+##   retsys: returned system
+##   nc: number of controllable states in the returned system
+##   no: number of observable states in the returned system
+##   cflg: is_controllable(retsys)
+##   oflg: is_observable(retsys)
+
+## Author: A. S. Hodel <a.s.hodel@eng.auburn.edu>
+
+function [retsys, nc, no, cflg, oflg] = sysmin (sys, flg);
+
+  switch(nargin)
+  case(1), flg = 0;
+  case(2), jnk = flg;    # dummy operation
+  otherwise,
+    usage("[retsys,nc,no] = sysmin(sys{,flg})");
+  endswitch
+  dflg = is_digital(sys,2);
+  [n,nz,m,p] = sysdimensions(sys);
+  if(n*nz > 0)
+    # both continuous and discrete states
+    [aa,bb,cc,dd,tsam,n,nz,stnam,innam,outnam,yd] = sys2ss(sys);
+    crng = 1:n;
+    drng = n+(1:nz);
+
+    # get minimal realization of continuous part
+    Ac  = aa(crng,crng);
+    Acd = aa(crng,drng);
+    Adc = aa(drng,crng);
+    Ad  = aa(drng,drng);
+    Bc  = bb(crng,:);
+    Bd  = bb(drng,:);
+    Cc  = cc(:,crng);
+    Cd  = cc(:,drng);
+
+    cstnam = stnam(crng);
+    dstnam = stnam(drng);
+    cinnam = append(innam,stnam(drng));
+    coutnam = append(outnam,stnam(drng));
+    csys = ss2sys(Ac,[Bc,Acd],[Cc;Adc]);
+    csys = syssetsignals(csys,"st",cstnam);
+    csys = syssetsignals(csys,"in",cinnam);
+    csys = syssetsignals(csys,"out",coutnam);
+
+    # reduce continuous system, recombine with discrete part
+    csys = sysmin(csys,flg);
+    cn = sysdimensions(csys);
+
+    if(cn == 0)
+      # continuous states are removed; just reduce the discrete part
+      sys = sysprune(sys,1:p,1:m,drng);
+      retsys = sysmin(sys,flg);
+    else
+      # extract updated parameters from reduced continuous system
+      [caa,cbb,ccc,cdd,ctsam,cn,cnz,cstnam,cinnam,coutnam] = sys2ss(csys);
+      crng = 1:cn;
+      Ac  = caa;
+      Bc  = cbb(:,1:m);
+      Acd = cbb(:,m+(1:nz));
+      Cc  = ccc(1:p,:);
+      Adc = ccc(p + (1:nz),:);
+
+      # recombine to reduce discrete part of the system
+      dinnam = append(innam,cstnam);
+      doutnam = append(outnam,cstnam);
+      dsys = ss2sys(Ad,[Bd,Adc],[Cd;Acd],[],tsam);
+      dsys = syssetsignals(dsys,"st",dstnam);
+      dsys = syssetsignals(dsys,"in",dinnam);
+      dsys = syssetsignals(dsys,"out",doutnam);
+
+      # reduce discrete subsystem
+      dsys = sysmin(dsys);
+      [n1,nz] = sysdimensions(dsys);
+      if(nz == 0)
+        # discrete subsystem is not needed
+        retsys = sysprune(csys,1:p,1:m);
+      else
+        # combine discrete, continuous subsystems
+        [Ad,dbb,dcc] = sys2ss(dsys);
+        dstnam = sysgetsignals(dsys,"st");
+        Bd  = dbb(:,1:m);
+        Adc = dbb(:,m+(1:cn));
+        Cd  = dcc(1:p,:);
+        Acd = dcc(p+(1:cn),:);
+        stnam = append(cstnam,dstnam);
+        aa = [Ac, Acd; Adc, Ad];
+        bb = [Bc; Bd];
+        cc = [Cc, Cd];
+        retsys = ss2sys([Ac, Acd; Adc, Ad], [Bc ; Bd], [Cc, Cd], dd, tsam, ...
+          cn, nz, stnam, innam, outnam, find(yd == 1));
+      end
+    endif
+  else
+    Ts = sysgettsam(sys);
+    switch(flg)
+    case(0),
+      ## reduce to a minimal system
+      [aa,bb,cc,dd] = sys2ss(sys);
+      [cflg,Uc] = is_controllable(aa,bb);
+      if(!cflg)
+        ## reduce to controllable states
+        if(!isempty(Uc))
+          aa = Uc'*aa*Uc;
+          bb = Uc'*bb;
+          cc = cc*Uc;
+        else
+          aa = bb = cc = [];
+        endif
+      endif
+      if(!isempty(aa))
+        [oflg,Uo] = is_observable(aa,cc);
+        if(!oflg)
+          if(!isempty(Uo))
+            aa = Uo'*aa*Uo;
+            bb = Uo'*bb;
+            cc = cc*Uo;
+          else
+            aa = bb = cc = [];
+          endif
+        endif
+      endif
+      switch(dflg)
+      case(0),
+        nc = no = nn = columns(aa);
+        nz = 0;
+      case(1),
+        nc = no = nz = columns(aa);
+        nn = 0;
+      endswitch
+      innam = sysgetsignals(sys,"in");
+      outnam= sysgetsignals(sys,"out");
+      retsys = ss2sys(aa,bb,cc,dd,Ts,nn,nz,[],innam,outnam);
+    case(1),
+      ## reduced model with physical states
+      [cflg,Uc] = is_controllable(sys); xc = find(max(abs(Uc')) != 0);
+      [oflg,Uo] = is_observable(sys);   xo = find(max(abs(Uo')) != 0);
+      xx = intersection(xc,xo);
+      if(isempty(xx)) xx = 0;  endif    # signal no states in reduced model
+      retsys = sysprune(sys,[],[],xx);
+    otherwise,
+      error ("invalid value of flg = %d", flg);
+    endswitch
+    if(sysdimensions(retsys,"st") > 0)
+      [cflg,Uc] = is_controllable(retsys); nc = columns(Uc);
+      [oflg,Uo] = is_observable(retsys);   no = columns(Uo);
+    else
+      nc = no = 0;
+    endif
+  endif
+endfunction