new file mode 100644
--- /dev/null
+++ b/scripts/control/sysmin.m
@@ -0,0 +1,80 @@
+function [retsys,nc,no,cflg,oflg] = sysmin(sys,flg);
+  # [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)
+  
+  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,flg);
+  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
+    inname = sysgetsignals(sys,"in");
+    outname= sysgetsignals(sys,"out");
+    retsys = ss2sys(aa,bb,cc,dd,Ts,nn,nz,[],inname,outname);
+  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("illegal 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
+endfunction