--- a/scripts/control/sysmin.m
+++ b/scripts/control/sysmin.m
@@ -42,60 +42,136 @@
   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;
+  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
-        aa = bb = cc = [];
-      endif
+        # 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
-    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
+    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
-    endif
-    switch(dflg)
-    case(0),
-      nc = no = nn = columns(aa);
-      nz = 0;
+      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),
-      nc = no = nz = columns(aa);
-      nn = 0;
+      ## 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
-    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;
+    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