--- a/scripts/control/base/dgkfdemo.m	Fri Jan 14 09:52:35 2000 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,354 +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} {} dgkfdemo ()
-## Octave Controls toolbox demo: H2/Hinfinity options demos
-## @end deftypefn
-
-## Author: A. S. Hodel <a.s.hodel@eng.auburn.edu>
-## Created: June 1995
-
-function dgkfdemo ()
-
-  save_val = page_screen_output;
-  page_screen_output = 1;
-  while (1)
-    clc
-    sel = 0;
-    while (sel > 10 || sel < 1)
-      sel = menu ("Octave H2/Hinfinity options demo",
-                  "LQ regulator",
-                  "LG state estimator",
-                  "LQG optimal control design",
-                  "H2 gain of a system",
-                  "H2 optimal controller of a system",
-                  "Hinf gain of a system",
-                  "Hinf optimal controller of a SISO system",
-                  "Hinf optimal controller of a MIMO system",
-                  "Discrete-time Hinf optimal control by bilinear transform",
-                  "Return to main demo menu");
-    endwhile
-    if (sel == 1)
-      disp("Linear/Quadratic regulator design:")
-      disp("Compute optimal state feedback via the lqr command...")
-      help lqr
-      disp(" ")
-      disp("Example:")
-      A = [0, 1; -2, -1]
-      B = [0; 1]
-      Q = [1, 0; 0, 0]
-      R = 1
-      disp("Q = state penalty matrix; R = input penalty matrix")
-      prompt
-      disp("Compute state feedback gain k, ARE solution P, and closed-loop")
-      disp("poles as follows:");
-      cmd = "[k, p, e] = lqr(A,B,Q,R)";
-      run_cmd
-      prompt
-      disp("A similar approach can be used for LTI discrete-time systems")
-      disp("by using the dlqr command in place of lqr (see LQG example).")
-    elseif (sel == 2)
-      disp("Linear/Gaussian estimator design:")
-      disp("Compute optimal state estimator via the lqe command...")
-      help lqe
-      disp(" ")
-      disp("Example:")
-      A = [0, 1; -2, -1]
-      disp("disturbance entry matrix G")
-      G = eye(2)
-      disp("Output measurement matrix C")
-      C = [0, 1]
-      SigW = [1, 0; 0, 1]
-      SigV = 1
-      disp("SigW = input disturbance intensity matrix;")
-      disp("SigV = measurement noise intensity matrix")
-      prompt
-      disp("Compute estimator feedback gain k, ARE solution P, and estimator")
-      disp("poles via the command: ")
-      cmd = "[k, p, e] = lqe(A,G,C,SigW,SigV)";
-      run_cmd
-      disp("A similar approach can be used for LTI discrete-time systems")
-      disp("by using the dlqe command in place of lqe (see LQG example).")
-    elseif (sel == 3)
-      disp("LQG optimal controller of a system:")
-      disp("Input accepted as either A,B,C matrices or in system data structure form")
-      disp("in both discrete and continuous time.")
-      disp("Example 1: continuous time design:")
-      prompt
-      help lqg
-      disp("Example system")
-      A = [0, 1; .5, .5];
-      B = [0; 2];
-      G = eye(2)
-      C = [1, 1];
-      sys = ss2sys(A, [B, G], C);
-      sys = syssetsignals(sys,"in", ...
-                       ["control input"; "disturbance 1"; "disturbance 2"]);
-      sysout(sys)
-      prompt
-      disp("Filtering/estimator parameters:")
-      SigW = eye(2)
-      SigV = 1
-      prompt
-      disp("State space (LQR) parameters Q and R are:")
-      Q = eye(2)
-      R = 1
-      cmd = "[K,Q1,P1,Ee,Er] = lqg(sys,SigW,SigV,Q,R,1);";
-      run_cmd
-      disp("Check: closed loop system A-matrix is")
-      disp(" [A,      B*Cc]")
-      disp(" [Bc*C,   Ac  ]")
-      cmd = "[Ac, Bc, Cc] = sys2ss(K);";
-      run_cmd
-      cmd = "Acl = [A, B*Cc; Bc*C, Ac]";
-      run_cmd
-      disp("Check: poles of Acl:")
-      Acl_poles = sortcom(eig(Acl))
-      disp("Predicted poles from design = union(Er,Ee)")
-      cmd = "pred_poles = sortcom([Er; Ee])";
-      run_cmd
-      disp("Example 2: discrete-time example")
-      cmd1 = "Dsys = ss2sys(A, [G, B], C, [0, 0, 0], 1);";
-      cmd2 = "[K,Q1,P1,Ee,Er] = lqg(Dsys,SigW, SigV,Q,R);";
-      disp("Run commands:")
-      cmd = cmd1;
-      run_cmd
-      cmd = cmd2;
-      run_cmd
-      prompt
-      disp("Check: closed loop system A-matrix is")
-      disp(" [A,      B*Cc]")
-      disp(" [Bc*C,   Ac  ]")
-      [Ac,Bc,Cc] = sys2ss(K);
-      Acl = [A, B*Cc; Bc*C, Ac]
-      prompt
-      disp("Check: poles of Acl:")
-      Acl_poles = sortcom(eig(Acl))
-      disp("Predicted poles from design = union(Er,Ee)")
-      pred_poles = sortcom([Er;Ee])
-    elseif (sel == 4)
-      disp("H2 gain of a system: (Energy in impulse response)")
-      disp("Example 1: Stable plant:")
-      cmd = "A = [0, 1; -2, -1]; B = [0; 1]; C = [1, 0]; sys_poles = eig(A)";
-      run_cmd
-      disp("Put into Packed system form:")
-      cmd = "Asys = ss2sys(A,B,C);";
-      run_cmd
-      disp("Evaluate system 2-norm (impulse response energy):");
-      cmd = "AsysH2 = h2norm(Asys)";
-      run_cmd
-      disp("Compare with a plot of the system impulse response:")
-      tt = 0:0.1:20;
-      for ii=1:length(tt)
-        ht(ii) = C*expm(A*tt(ii))*B;
-      endfor
-      plot(tt,ht)
-      title("impulse response of example plant")
-      prompt
-      disp("Example 2: unstable plant")
-      cmd = "A = [0, 1; 2, 1]";
-      eval(cmd);
-      cmd = "B = [0; 1]";
-      eval(cmd);
-      cmd = "C = [1, 0]";
-      eval(cmd);
-      cmd = "sys_poles = eig(A)";
-      run_cmd
-      prompt
-      disp("Put into system data structure form:")
-      cmd="Bsys = ss2sys(A,B,C);";
-      run_cmd
-      disp("Evaluate 2-norm:")
-      cmd = "BsysH2 = h2norm(Bsys)";
-      run_cmd
-      disp(" ")
-      prompt("NOTICE: program returns a value without an error signal.")
-      disp("")
-
-    elseif (sel == 5)
-      disp("H2 optimal controller of a system: command = h2syn:")
-      prompt
-      help h2syn
-      prompt
-      disp("Example system: double integrator with output noise and")
-      disp("input disturbance:")
-      disp(" ");
-      disp("       -------------------->y2");
-      disp("       |   _________");
-      disp("u(t)-->o-->| 1/s^2 |-->o-> y1");
-      disp("       ^   ---------   ^");
-      disp("       |               |");
-      disp("      w1(t)           w2(t)");
-      disp(" ")
-      disp("w enters the system through B1, u through B2")
-      disp("z = [y1; y2] is obtained through C1, y=y1 through C2");
-      disp(" ")
-      cmd = "A = [0, 1; 0, 0];  B1 = [0, 0; 1, 0]; B2 = [0; 1];";
-      disp(cmd)
-      eval(cmd);
-      cmd = "C1 = [1, 0; 0, 0]; C2 = [1, 0];    D11 = zeros(2);";
-      disp(cmd)
-      eval(cmd);
-      cmd = "D12 = [0; 1];  D21 = [0, 1];  D22 = 0; D = [D11, D12; D21, D22];";
-      disp(cmd)
-      eval(cmd);
-      disp("Design objective: compute U(s)=K(s)Y1(s) to minimize the closed")
-      disp("loop impulse response from w(t) =[w1; w2] to z(t) = [y1; y2]");
-      prompt
-      disp("First: pack system:")
-      cmd="Asys = ss2sys(A, [B1, B2], [C1; C2], D);";
-      run_cmd
-      disp("Open loop multivariable Bode plot: (will take a moment)")
-      cmd="bode(Asys);";
-      run_cmd
-      prompt("Press a key to close plot and continue");
-      closeplot
-      disp("Controller design command: (only need 1st two output arguments)")
-      cmd="[K,gain, Kc, Kf, Pc,  Pf] = h2syn(Asys,1,1);";
-      run_cmd
-      disp("Controller is:")
-      cmd = "sysout(K)";
-      run_cmd
-      disp(["returned gain value is: ",num2str(gain)]);
-      disp("Check: close the loop and then compute h2norm:")
-      prompt
-      cmd="K_loop = sysgroup(Asys,K);";
-      run_cmd
-      cmd = "Kcl = sysconnect(K_loop,[3,4],[4,3]);";
-      run_cmd
-      cmd = "Kcl = sysprune(Kcl,[1,2],[1,2]);";
-      run_cmd
-      cmd="gain_Kcl = h2norm(Kcl)";
-      run_cmd
-      cmd="gain_err = gain_Kcl - gain";
-      run_cmd
-      disp("Check: multivarible bode plot:")
-      cmd="bode(Kcl);";
-      run_cmd
-      prompt
-      disp("Related functions: is_dgkf, is_controllable, is_stabilizable,")
-      disp("                is_observable, is_detectable")
-    elseif (sel == 6)
-      disp("Hinfinity gain of a system: (max gain over all j-omega)")
-      disp("Example 1: Stable plant:")
-      cmd = "A = [0, 1; -2, -1]; B = [0; 1]; C = [1, 0]; sys_poles = eig(A)";
-      run_cmd
-      disp("Pack into system format:")
-      cmd = "Asys = ss2sys(A,B,C);";
-      run_cmd
-      disp("The infinity norm must be computed iteratively by")
-      disp("binary search.  For this example, we select tolerance tol = 0.01, ")
-      disp("min gain gmin = 1e-2, max gain gmax=1e4.")
-      disp("Search quits when upper bound <= (1+tol)*lower bound.")
-      cmd = "tol = 0.01; gmin = 1e-2; gmax = 1e+4;";
-      run_cmd
-      cmd = "[AsysHinf,gmin,gmax] = hinfnorm(Asys,tol,gmin,gmax)"
-      run_cmd
-      disp("Check: look at max value of magntude Bode plot of Asys:");
-      [M,P,w] = bode(Asys);
-      xlabel("Omega")
-      ylabel("|Asys(j omega)| ")
-      grid();
-      semilogx(w,M);
-      disp(["Max magnitude is ",num2str(max(M)), ...
-        ", compared with gmin=",num2str(gmin)," and gmax=", ...
-        num2str(gmax),"."])
-      prompt
-      disp("Example 2: unstable plant")
-      cmd = "A = [0, 1; 2, 1]; B = [0; 1]; C = [1, 0]; sys_poles = eig(A)";
-      run_cmd
-      disp("Pack into system format:")
-      cmd = "Bsys = ss2sys(A,B,C);";
-      run_cmd
-      disp("Evaluate with BsysH2 = hinfnorm(Bsys,tol,gmin,gmax)")
-      BsysH2 = hinfnorm(Bsys,tol,gmin,gmax)
-      disp(" ")
-      disp("NOTICE: program returns a value without an error signal.")
-      disp("")
-
-    elseif (sel == 7)
-      disp("Hinfinity optimal controller of a system: command = hinfsyn:")
-      prompt
-      help hinfsyn
-      prompt
-      disp("Example system: double integrator with output noise and")
-      disp("input disturbance:")
-      A = [0, 1; 0, 0]
-      B1 = [0, 0; 1, 0]
-      B2 = [0; 1]
-      C1 = [1, 0; 0, 0]
-      C2 = [1, 0]
-      D11 = zeros(2);
-      D12 = [0; 1];
-      D21 = [0, 1];
-      D22 = 0;
-      D = [D11, D12; D21, D22]
-      prompt
-      disp("First: pack system:")
-      cmd="Asys = ss2sys(A, [B1, B2], [C1; C2], D);";
-      run_cmd
-      prompt
-      disp("Open loop multivariable Bode plot: (will take a moment)")
-      cmd="bode(Asys);";
-      run_cmd
-      prompt
-      disp("Controller design command: (only need 1st two output arguments)")
-      gmax = 1000
-      gmin = 0.1
-      gtol = 0.01
-      cmd="[K,gain] = hinfsyn(Asys,1,1,gmin,gmax,gtol);";
-      run_cmd
-      disp("Check: close the loop and then compute h2norm:")
-      prompt
-      cmd="K_loop = sysgroup(Asys,K);";
-      run_cmd
-      cmd = "Kcl = sysconnect(K_loop,[3,4],[4,3]);";
-      run_cmd
-      cmd = "Kcl = sysprune(Kcl,[1,2],[1,2]);";
-      run_cmd
-      cmd="gain_Kcl = hinfnorm(Kcl)";
-      run_cmd
-      cmd="gain_err = gain_Kcl - gain";
-      run_cmd
-      disp("Check: multivarible bode plot:")
-      cmd="bode(Kcl);";
-      run_cmd
-      prompt
-      disp("Related functions: is_dgkf, is_controllable, is_stabilizable,")
-      disp("                   is_observable, is_detectable, buildssic")
-    elseif (sel == 8)
-      disp("Hinfinity optimal controller of MIMO system: command = hinfsyn:")
-      prompt
-      help hinfsyn
-      prompt
-      disp("Example system: Boeing 707-321 airspeed/pitch angle control")
-      disp(" ")
-      hinfdemo
-    elseif (sel == 9)
-      disp("Discrete time H-infinity control via bilinear transform");
-      prompt
-      dhinfdemo
-    elseif (sel == 10)
-      return
-    endif
-    prompt
-  endwhile
-  page_screen_output = save_val;
-
-endfunction

--- a/scripts/control/base/is_dgkf.m	Fri Jan 14 09:52:35 2000 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,258 +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{retval}, @var{dgkf_struct} ] =} is_dgkf (@var{Asys}, @var{nu}, @var{ny}, @var{tol} )
-## Determine whether a continuous time state space system meets
-## assumptions of DGKF algorithm.
-## Partitions system into:
-## @example
-## [dx/dt] = [A  | Bw  Bu  ][w]
-## [ z   ]   [Cz | Dzw Dzu ][u]
-## [ y   ]   [Cy | Dyw Dyu ]
-## @end example
-## or similar discrete-time system.
-## If necessary, orthogonal transformations @var{Qw}, @var{Qz} and nonsingular
-## transformations @var{Ru}, @var{Ry} are applied to respective vectors
-## @var{w}, @var{z}, @var{u}, @var{y} in order to satisfy DGKF assumptions.
-## Loop shifting is used if @var{Dyu} block is nonzero.
-##
-## @strong{Inputs}
-## @table @var
-## @item         Asys
-## system data structure
-## @item           nu
-## number of controlled inputs
-## @item        ny
-## number of measured outputs
-## @item        tol
-## threshhold for 0.  Default: 200@var{eps}
-## @end table
-## @strong{Outputs}
-## @table @var
-## @item    retval
-## true(1) if system passes check, false(0) otherwise
-## @item    dgkf_struct
-## data structure of @code{is_dgkf} results.  Entries:
-## @table @var
-## @item      nw
-## @itemx     nz
-## dimensions of @var{w}, @var{z}
-## @item      A
-## system @var{A} matrix
-## @item      Bw
-## (@var{n} x @var{nw}) @var{Qw}-transformed disturbance input matrix
-## @item      Bu
-## (@var{n} x @var{nu}) @var{Ru}-transformed controlled input matrix;
-##
-## @strong{Note} @math{B = [Bw Bu] }
-## @item      Cz
-## (@var{nz} x @var{n}) Qz-transformed error output matrix
-## @item      Cy
-## (@var{ny} x @var{n}) @var{Ry}-transformed measured output matrix
-##
-## @strong{Note} @math{C = [Cz; Cy] }
-## @item      Dzu
-## @item      Dyw
-## off-diagonal blocks of transformed @var{D} matrix that enter
-## @var{z}, @var{y} from @var{u}, @var{w} respectively
-## @item      Ru
-## controlled input transformation matrix
-## @item      Ry
-## observed output transformation matrix
-## @item      Dyu_nz
-## nonzero if the @var{Dyu} block is nonzero.
-## @item      Dyu
-## untransformed @var{Dyu} block
-## @item      dflg
-## nonzero if the system is discrete-time
-## @end table
-## @end table
-## @code{is_dgkf} exits with an error if the system is mixed
-## discrete/continuous
-##
-## @strong{References}
-## @table @strong
-## @item [1]
-## Doyle, Glover, Khargonekar, Francis, "State Space Solutions
-## to Standard H2 and Hinf Control Problems," IEEE TAC August 1989
-## @item [2]
-## Maciejowksi, J.M.: "Multivariable feedback design,"
-## @end table
-## @end deftypefn
-
-## Author: A. S. Hodel <a.s.hodel@eng.auburn.edu>
-## Updated by John Ingram July 1996 to accept structured systems
-
-## Revised by Kai P. Mueller April 1998 to solve the general H_infinity
-## problem using unitary transformations Q (on w and z)
-## and non-singular transformations R (on u and y) such
-## that the Dzu and Dyw matrices of the transformed plant
-##
-##    ~
-##    P  (the variable Asys here)
-##
-## become
-##
-##    ~            -1         T
-##    D  = Q   D   R   = [ 0 I ]  or [ I ],
-##     12   12  12  12
-##
-##    ~            T
-##    D  = R   D   Q   = [ 0 I ] or [ I ].
-##     21   21  21  21
-##
-## This transformation together with the algorithm in [1] solves
-## the general problem (see [2] for example).
-
-function [retval, dgkf_struct] = is_dgkf (Asys, nu, ny, tol)
-
-  if (nargin < 3) | (nargin > 4)
-    usage("[retval,dgkf_struct] = is_dgkf(Asys,nu,ny{,tol})");
-  elseif (! is_scalar(nu) | ! is_scalar(ny) )
-    error("is_dgkf: arguments 2 and 3 must be scalars")
-  elseif (! is_struct(Asys) )
-    error("Argument 1 must be a system data structure");
-  endif
-  if(nargin < 4)
-    tol = 200*eps;
-  elseif( !is_sample(tol) )
-    error("is_dgkf: tol must be a positive scalar")
-  endif
-
-  retval = 1;           # assume passes test
-
-  dflg = is_digital(Asys);
-  [Anc, Anz, nin, nout ] = sysdimensions(Asys);
-
-  if( Anz == 0 & Anc == 0 )
-    error("is_dgkf: no system states");
-  elseif( nu >= nin )
-    error("is_dgkf: insufficient number of disturbance inputs");
-  elseif( ny >= nout )
-    error("is_dgkf: insufficient number of regulated outputs");
-  endif
-
-  nw = nin - nu;           nw1 = nw + 1;
-  nz = nout - ny;          nz1 = nz + 1;
-
-  [A,B,C,D] = sys2ss(Asys);
-  ## scale input/output for numerical reasons
-  if(norm (C, "fro") * norm (B, "fro") == 0)
-    error("||C||*||B|| = 0; no dynamic connnection from inputs to outputs");
-  endif
-  xx = sqrt(norm(B, Inf) / norm(C, Inf));
-  B = B / xx;  C = C * xx;
-
-  ## partition matrices
-                        Bw = B(:,1:nw);         Bu = B(:,nw1:nin);
-  Cz = C(1:nz,:);       Dzw = D(1:nz,1:nw);     Dzu = D(1:nz,nw1:nin);
-  Cy = C(nz1:nout,:);   Dyw = D(nz1:nout,1:nw); Dyu = D(nz1:nout,nw1:nin);
-
-  ## Check for loopo shifting
-  Dyu_nz = (norm(Dyu,Inf) != 0);
-  if (Dyu_nz)
-    warning("is_dgkf: D22 nonzero; performing loop shifting");
-  endif
-
-  ## 12 - rank condition at w = 0
-  xx =[A, Bu; Cz, Dzu];
-  [nr, nc] = size(xx);
-  irank = rank(xx);
-  if (irank != nc)
-    retval = 0;
-    warning(sprintf("rank([A Bu; Cz Dzu]) = %d, need %d; n=%d, nz=%d, nu=%d", ...
-        irank,nc,(Anc+Anz),nz,nu));
-    warning(" *** 12-rank condition violated at w = 0.");
-  endif
-
-  ## 21 - rank condition at w = 0
-  xx =[A, Bw; Cy, Dyw];
-  [nr, nc] = size(xx);
-  irank = rank(xx);
-  if (irank != nr)
-    retval = 0;
-    warning(sprintf("rank([A Bw; Cy Dyw]) = %d, need %d; n=%d, ny=%d, nw=%d", ...
-        irank,nr,(Anc+Anz),ny,nw));
-    warning(" *** 21-rank condition violated at w = 0.");
-  endif
-
-  ## can Dzu be transformed to become [0 I]' or [I]?
-  ## This ensures a normalized weight
-  [Qz, Ru] = qr(Dzu);
-  irank = rank(Ru);
-  if (irank != nu)
-    retval = 0;
-    warning(sprintf("*** rank(Dzu(%d x %d) = %d", nz, nu, irank));
-    warning(" *** Dzu does not have full column rank.");
-  endif
-  if (nu >= nz)
-    Qz = Qz(:,1:nu)';
-  else
-    Qz = [Qz(:,(nu+1):nz), Qz(:,1:nu)]';
-  endif
-  Ru = Ru(1:nu,:);
-
-  ## can Dyw be transformed to become [0 I] or [I]?
-  ## This ensures a normalized weight
-  [Qw, Ry] = qr(Dyw');
-  irank = rank(Ry);
-  if (irank != ny)
-    retval = 0;
-    warning(sprintf("*** rank(Dyw(%d x %d) = %d", ny, nw, irank));
-    warning(" *** Dyw does not have full row rank.");
-  endif
-
-  if (ny >= nw)
-    Qw = Qw(:,1:ny);
-  else
-    Qw = [Qw(:,(ny+1):nw), Qw(:,1:ny)];
-  endif
-  Ry = Ry(1:ny,:)';
-
-  ## transform P by Qz/Ru and Qw/Ry
-  Bw  = Bw*Qw;
-  Bu  = Bu/Ru;
-  B   = [Bw, Bu];
-  Cz  = Qz*Cz;
-  Cy  = Ry\Cy;
-  C   = [Cz; Cy];
-  Dzw = Qz*Dzw*Qw;
-  Dzu = Qz*Dzu/Ru;
-  Dyw = Ry\Dyw*Qw;
-
-  ## pack the return structure
-  dgkf_struct.nw        = nw;
-  dgkf_struct.nu        = nu;
-  dgkf_struct.nz        = nz;
-  dgkf_struct.ny        = ny;
-  dgkf_struct.A         = A;
-  dgkf_struct.Bw        = Bw;
-  dgkf_struct.Bu        = Bu;
-  dgkf_struct.Cz        = Cz;
-  dgkf_struct.Cy        = Cy;
-  dgkf_struct.Dzw       = Dzw;
-  dgkf_struct.Dzu       = Dzu;
-  dgkf_struct.Dyw       = Dyw;
-  dgkf_struct.Dyu       = Dyu;
-  dgkf_struct.Ru        = Ru;
-  dgkf_struct.Ry        = Ry;
-  dgkf_struct.Dyu_nz    = Dyu_nz;
-  dgkf_struct.dflg      = dflg;
-
-endfunction

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/control/hinf/dgkfdemo.m	Fri Jan 14 22:00:20 2000 +0000
@@ -0,0 +1,354 @@
+## 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} {} dgkfdemo ()
+## Octave Controls toolbox demo: H2/Hinfinity options demos
+## @end deftypefn
+
+## Author: A. S. Hodel <a.s.hodel@eng.auburn.edu>
+## Created: June 1995
+
+function dgkfdemo ()
+
+  save_val = page_screen_output;
+  page_screen_output = 1;
+  while (1)
+    clc
+    sel = 0;
+    while (sel > 10 || sel < 1)
+      sel = menu ("Octave H2/Hinfinity options demo",
+                  "LQ regulator",
+                  "LG state estimator",
+                  "LQG optimal control design",
+                  "H2 gain of a system",
+                  "H2 optimal controller of a system",
+                  "Hinf gain of a system",
+                  "Hinf optimal controller of a SISO system",
+                  "Hinf optimal controller of a MIMO system",
+                  "Discrete-time Hinf optimal control by bilinear transform",
+                  "Return to main demo menu");
+    endwhile
+    if (sel == 1)
+      disp("Linear/Quadratic regulator design:")
+      disp("Compute optimal state feedback via the lqr command...")
+      help lqr
+      disp(" ")
+      disp("Example:")
+      A = [0, 1; -2, -1]
+      B = [0; 1]
+      Q = [1, 0; 0, 0]
+      R = 1
+      disp("Q = state penalty matrix; R = input penalty matrix")
+      prompt
+      disp("Compute state feedback gain k, ARE solution P, and closed-loop")
+      disp("poles as follows:");
+      cmd = "[k, p, e] = lqr(A,B,Q,R)";
+      run_cmd
+      prompt
+      disp("A similar approach can be used for LTI discrete-time systems")
+      disp("by using the dlqr command in place of lqr (see LQG example).")
+    elseif (sel == 2)
+      disp("Linear/Gaussian estimator design:")
+      disp("Compute optimal state estimator via the lqe command...")
+      help lqe
+      disp(" ")
+      disp("Example:")
+      A = [0, 1; -2, -1]
+      disp("disturbance entry matrix G")
+      G = eye(2)
+      disp("Output measurement matrix C")
+      C = [0, 1]
+      SigW = [1, 0; 0, 1]
+      SigV = 1
+      disp("SigW = input disturbance intensity matrix;")
+      disp("SigV = measurement noise intensity matrix")
+      prompt
+      disp("Compute estimator feedback gain k, ARE solution P, and estimator")
+      disp("poles via the command: ")
+      cmd = "[k, p, e] = lqe(A,G,C,SigW,SigV)";
+      run_cmd
+      disp("A similar approach can be used for LTI discrete-time systems")
+      disp("by using the dlqe command in place of lqe (see LQG example).")
+    elseif (sel == 3)
+      disp("LQG optimal controller of a system:")
+      disp("Input accepted as either A,B,C matrices or in system data structure form")
+      disp("in both discrete and continuous time.")
+      disp("Example 1: continuous time design:")
+      prompt
+      help lqg
+      disp("Example system")
+      A = [0, 1; .5, .5];
+      B = [0; 2];
+      G = eye(2)
+      C = [1, 1];
+      sys = ss2sys(A, [B, G], C);
+      sys = syssetsignals(sys,"in", ...
+                       ["control input"; "disturbance 1"; "disturbance 2"]);
+      sysout(sys)
+      prompt
+      disp("Filtering/estimator parameters:")
+      SigW = eye(2)
+      SigV = 1
+      prompt
+      disp("State space (LQR) parameters Q and R are:")
+      Q = eye(2)
+      R = 1
+      cmd = "[K,Q1,P1,Ee,Er] = lqg(sys,SigW,SigV,Q,R,1);";
+      run_cmd
+      disp("Check: closed loop system A-matrix is")
+      disp(" [A,      B*Cc]")
+      disp(" [Bc*C,   Ac  ]")
+      cmd = "[Ac, Bc, Cc] = sys2ss(K);";
+      run_cmd
+      cmd = "Acl = [A, B*Cc; Bc*C, Ac]";
+      run_cmd
+      disp("Check: poles of Acl:")
+      Acl_poles = sortcom(eig(Acl))
+      disp("Predicted poles from design = union(Er,Ee)")
+      cmd = "pred_poles = sortcom([Er; Ee])";
+      run_cmd
+      disp("Example 2: discrete-time example")
+      cmd1 = "Dsys = ss2sys(A, [G, B], C, [0, 0, 0], 1);";
+      cmd2 = "[K,Q1,P1,Ee,Er] = lqg(Dsys,SigW, SigV,Q,R);";
+      disp("Run commands:")
+      cmd = cmd1;
+      run_cmd
+      cmd = cmd2;
+      run_cmd
+      prompt
+      disp("Check: closed loop system A-matrix is")
+      disp(" [A,      B*Cc]")
+      disp(" [Bc*C,   Ac  ]")
+      [Ac,Bc,Cc] = sys2ss(K);
+      Acl = [A, B*Cc; Bc*C, Ac]
+      prompt
+      disp("Check: poles of Acl:")
+      Acl_poles = sortcom(eig(Acl))
+      disp("Predicted poles from design = union(Er,Ee)")
+      pred_poles = sortcom([Er;Ee])
+    elseif (sel == 4)
+      disp("H2 gain of a system: (Energy in impulse response)")
+      disp("Example 1: Stable plant:")
+      cmd = "A = [0, 1; -2, -1]; B = [0; 1]; C = [1, 0]; sys_poles = eig(A)";
+      run_cmd
+      disp("Put into Packed system form:")
+      cmd = "Asys = ss2sys(A,B,C);";
+      run_cmd
+      disp("Evaluate system 2-norm (impulse response energy):");
+      cmd = "AsysH2 = h2norm(Asys)";
+      run_cmd
+      disp("Compare with a plot of the system impulse response:")
+      tt = 0:0.1:20;
+      for ii=1:length(tt)
+        ht(ii) = C*expm(A*tt(ii))*B;
+      endfor
+      plot(tt,ht)
+      title("impulse response of example plant")
+      prompt
+      disp("Example 2: unstable plant")
+      cmd = "A = [0, 1; 2, 1]";
+      eval(cmd);
+      cmd = "B = [0; 1]";
+      eval(cmd);
+      cmd = "C = [1, 0]";
+      eval(cmd);
+      cmd = "sys_poles = eig(A)";
+      run_cmd
+      prompt
+      disp("Put into system data structure form:")
+      cmd="Bsys = ss2sys(A,B,C);";
+      run_cmd
+      disp("Evaluate 2-norm:")
+      cmd = "BsysH2 = h2norm(Bsys)";
+      run_cmd
+      disp(" ")
+      prompt("NOTICE: program returns a value without an error signal.")
+      disp("")
+
+    elseif (sel == 5)
+      disp("H2 optimal controller of a system: command = h2syn:")
+      prompt
+      help h2syn
+      prompt
+      disp("Example system: double integrator with output noise and")
+      disp("input disturbance:")
+      disp(" ");
+      disp("       -------------------->y2");
+      disp("       |   _________");
+      disp("u(t)-->o-->| 1/s^2 |-->o-> y1");
+      disp("       ^   ---------   ^");
+      disp("       |               |");
+      disp("      w1(t)           w2(t)");
+      disp(" ")
+      disp("w enters the system through B1, u through B2")
+      disp("z = [y1; y2] is obtained through C1, y=y1 through C2");
+      disp(" ")
+      cmd = "A = [0, 1; 0, 0];  B1 = [0, 0; 1, 0]; B2 = [0; 1];";
+      disp(cmd)
+      eval(cmd);
+      cmd = "C1 = [1, 0; 0, 0]; C2 = [1, 0];    D11 = zeros(2);";
+      disp(cmd)
+      eval(cmd);
+      cmd = "D12 = [0; 1];  D21 = [0, 1];  D22 = 0; D = [D11, D12; D21, D22];";
+      disp(cmd)
+      eval(cmd);
+      disp("Design objective: compute U(s)=K(s)Y1(s) to minimize the closed")
+      disp("loop impulse response from w(t) =[w1; w2] to z(t) = [y1; y2]");
+      prompt
+      disp("First: pack system:")
+      cmd="Asys = ss2sys(A, [B1, B2], [C1; C2], D);";
+      run_cmd
+      disp("Open loop multivariable Bode plot: (will take a moment)")
+      cmd="bode(Asys);";
+      run_cmd
+      prompt("Press a key to close plot and continue");
+      closeplot
+      disp("Controller design command: (only need 1st two output arguments)")
+      cmd="[K,gain, Kc, Kf, Pc,  Pf] = h2syn(Asys,1,1);";
+      run_cmd
+      disp("Controller is:")
+      cmd = "sysout(K)";
+      run_cmd
+      disp(["returned gain value is: ",num2str(gain)]);
+      disp("Check: close the loop and then compute h2norm:")
+      prompt
+      cmd="K_loop = sysgroup(Asys,K);";
+      run_cmd
+      cmd = "Kcl = sysconnect(K_loop,[3,4],[4,3]);";
+      run_cmd
+      cmd = "Kcl = sysprune(Kcl,[1,2],[1,2]);";
+      run_cmd
+      cmd="gain_Kcl = h2norm(Kcl)";
+      run_cmd
+      cmd="gain_err = gain_Kcl - gain";
+      run_cmd
+      disp("Check: multivarible bode plot:")
+      cmd="bode(Kcl);";
+      run_cmd
+      prompt
+      disp("Related functions: is_dgkf, is_controllable, is_stabilizable,")
+      disp("                is_observable, is_detectable")
+    elseif (sel == 6)
+      disp("Hinfinity gain of a system: (max gain over all j-omega)")
+      disp("Example 1: Stable plant:")
+      cmd = "A = [0, 1; -2, -1]; B = [0; 1]; C = [1, 0]; sys_poles = eig(A)";
+      run_cmd
+      disp("Pack into system format:")
+      cmd = "Asys = ss2sys(A,B,C);";
+      run_cmd
+      disp("The infinity norm must be computed iteratively by")
+      disp("binary search.  For this example, we select tolerance tol = 0.01, ")
+      disp("min gain gmin = 1e-2, max gain gmax=1e4.")
+      disp("Search quits when upper bound <= (1+tol)*lower bound.")
+      cmd = "tol = 0.01; gmin = 1e-2; gmax = 1e+4;";
+      run_cmd
+      cmd = "[AsysHinf,gmin,gmax] = hinfnorm(Asys,tol,gmin,gmax)"
+      run_cmd
+      disp("Check: look at max value of magntude Bode plot of Asys:");
+      [M,P,w] = bode(Asys);
+      xlabel("Omega")
+      ylabel("|Asys(j omega)| ")
+      grid();
+      semilogx(w,M);
+      disp(["Max magnitude is ",num2str(max(M)), ...
+        ", compared with gmin=",num2str(gmin)," and gmax=", ...
+        num2str(gmax),"."])
+      prompt
+      disp("Example 2: unstable plant")
+      cmd = "A = [0, 1; 2, 1]; B = [0; 1]; C = [1, 0]; sys_poles = eig(A)";
+      run_cmd
+      disp("Pack into system format:")
+      cmd = "Bsys = ss2sys(A,B,C);";
+      run_cmd
+      disp("Evaluate with BsysH2 = hinfnorm(Bsys,tol,gmin,gmax)")
+      BsysH2 = hinfnorm(Bsys,tol,gmin,gmax)
+      disp(" ")
+      disp("NOTICE: program returns a value without an error signal.")
+      disp("")
+
+    elseif (sel == 7)
+      disp("Hinfinity optimal controller of a system: command = hinfsyn:")
+      prompt
+      help hinfsyn
+      prompt
+      disp("Example system: double integrator with output noise and")
+      disp("input disturbance:")
+      A = [0, 1; 0, 0]
+      B1 = [0, 0; 1, 0]
+      B2 = [0; 1]
+      C1 = [1, 0; 0, 0]
+      C2 = [1, 0]
+      D11 = zeros(2);
+      D12 = [0; 1];
+      D21 = [0, 1];
+      D22 = 0;
+      D = [D11, D12; D21, D22]
+      prompt
+      disp("First: pack system:")
+      cmd="Asys = ss2sys(A, [B1, B2], [C1; C2], D);";
+      run_cmd
+      prompt
+      disp("Open loop multivariable Bode plot: (will take a moment)")
+      cmd="bode(Asys);";
+      run_cmd
+      prompt
+      disp("Controller design command: (only need 1st two output arguments)")
+      gmax = 1000
+      gmin = 0.1
+      gtol = 0.01
+      cmd="[K,gain] = hinfsyn(Asys,1,1,gmin,gmax,gtol);";
+      run_cmd
+      disp("Check: close the loop and then compute h2norm:")
+      prompt
+      cmd="K_loop = sysgroup(Asys,K);";
+      run_cmd
+      cmd = "Kcl = sysconnect(K_loop,[3,4],[4,3]);";
+      run_cmd
+      cmd = "Kcl = sysprune(Kcl,[1,2],[1,2]);";
+      run_cmd
+      cmd="gain_Kcl = hinfnorm(Kcl)";
+      run_cmd
+      cmd="gain_err = gain_Kcl - gain";
+      run_cmd
+      disp("Check: multivarible bode plot:")
+      cmd="bode(Kcl);";
+      run_cmd
+      prompt
+      disp("Related functions: is_dgkf, is_controllable, is_stabilizable,")
+      disp("                   is_observable, is_detectable, buildssic")
+    elseif (sel == 8)
+      disp("Hinfinity optimal controller of MIMO system: command = hinfsyn:")
+      prompt
+      help hinfsyn
+      prompt
+      disp("Example system: Boeing 707-321 airspeed/pitch angle control")
+      disp(" ")
+      hinfdemo
+    elseif (sel == 9)
+      disp("Discrete time H-infinity control via bilinear transform");
+      prompt
+      dhinfdemo
+    elseif (sel == 10)
+      return
+    endif
+    prompt
+  endwhile
+  page_screen_output = save_val;
+
+endfunction

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/control/hinf/is_dgkf.m	Fri Jan 14 22:00:20 2000 +0000
@@ -0,0 +1,258 @@
+## 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{retval}, @var{dgkf_struct} ] =} is_dgkf (@var{Asys}, @var{nu}, @var{ny}, @var{tol} )
+## Determine whether a continuous time state space system meets
+## assumptions of DGKF algorithm.
+## Partitions system into:
+## @example
+## [dx/dt] = [A  | Bw  Bu  ][w]
+## [ z   ]   [Cz | Dzw Dzu ][u]
+## [ y   ]   [Cy | Dyw Dyu ]
+## @end example
+## or similar discrete-time system.
+## If necessary, orthogonal transformations @var{Qw}, @var{Qz} and nonsingular
+## transformations @var{Ru}, @var{Ry} are applied to respective vectors
+## @var{w}, @var{z}, @var{u}, @var{y} in order to satisfy DGKF assumptions.
+## Loop shifting is used if @var{Dyu} block is nonzero.
+##
+## @strong{Inputs}
+## @table @var
+## @item         Asys
+## system data structure
+## @item           nu
+## number of controlled inputs
+## @item        ny
+## number of measured outputs
+## @item        tol
+## threshhold for 0.  Default: 200@var{eps}
+## @end table
+## @strong{Outputs}
+## @table @var
+## @item    retval
+## true(1) if system passes check, false(0) otherwise
+## @item    dgkf_struct
+## data structure of @code{is_dgkf} results.  Entries:
+## @table @var
+## @item      nw
+## @itemx     nz
+## dimensions of @var{w}, @var{z}
+## @item      A
+## system @var{A} matrix
+## @item      Bw
+## (@var{n} x @var{nw}) @var{Qw}-transformed disturbance input matrix
+## @item      Bu
+## (@var{n} x @var{nu}) @var{Ru}-transformed controlled input matrix;
+##
+## @strong{Note} @math{B = [Bw Bu] }
+## @item      Cz
+## (@var{nz} x @var{n}) Qz-transformed error output matrix
+## @item      Cy
+## (@var{ny} x @var{n}) @var{Ry}-transformed measured output matrix
+##
+## @strong{Note} @math{C = [Cz; Cy] }
+## @item      Dzu
+## @item      Dyw
+## off-diagonal blocks of transformed @var{D} matrix that enter
+## @var{z}, @var{y} from @var{u}, @var{w} respectively
+## @item      Ru
+## controlled input transformation matrix
+## @item      Ry
+## observed output transformation matrix
+## @item      Dyu_nz
+## nonzero if the @var{Dyu} block is nonzero.
+## @item      Dyu
+## untransformed @var{Dyu} block
+## @item      dflg
+## nonzero if the system is discrete-time
+## @end table
+## @end table
+## @code{is_dgkf} exits with an error if the system is mixed
+## discrete/continuous
+##
+## @strong{References}
+## @table @strong
+## @item [1]
+## Doyle, Glover, Khargonekar, Francis, "State Space Solutions
+## to Standard H2 and Hinf Control Problems," IEEE TAC August 1989
+## @item [2]
+## Maciejowksi, J.M.: "Multivariable feedback design,"
+## @end table
+## @end deftypefn
+
+## Author: A. S. Hodel <a.s.hodel@eng.auburn.edu>
+## Updated by John Ingram July 1996 to accept structured systems
+
+## Revised by Kai P. Mueller April 1998 to solve the general H_infinity
+## problem using unitary transformations Q (on w and z)
+## and non-singular transformations R (on u and y) such
+## that the Dzu and Dyw matrices of the transformed plant
+##
+##    ~
+##    P  (the variable Asys here)
+##
+## become
+##
+##    ~            -1         T
+##    D  = Q   D   R   = [ 0 I ]  or [ I ],
+##     12   12  12  12
+##
+##    ~            T
+##    D  = R   D   Q   = [ 0 I ] or [ I ].
+##     21   21  21  21
+##
+## This transformation together with the algorithm in [1] solves
+## the general problem (see [2] for example).
+
+function [retval, dgkf_struct] = is_dgkf (Asys, nu, ny, tol)
+
+  if (nargin < 3) | (nargin > 4)
+    usage("[retval,dgkf_struct] = is_dgkf(Asys,nu,ny{,tol})");
+  elseif (! is_scalar(nu) | ! is_scalar(ny) )
+    error("is_dgkf: arguments 2 and 3 must be scalars")
+  elseif (! is_struct(Asys) )
+    error("Argument 1 must be a system data structure");
+  endif
+  if(nargin < 4)
+    tol = 200*eps;
+  elseif( !is_sample(tol) )
+    error("is_dgkf: tol must be a positive scalar")
+  endif
+
+  retval = 1;           # assume passes test
+
+  dflg = is_digital(Asys);
+  [Anc, Anz, nin, nout ] = sysdimensions(Asys);
+
+  if( Anz == 0 & Anc == 0 )
+    error("is_dgkf: no system states");
+  elseif( nu >= nin )
+    error("is_dgkf: insufficient number of disturbance inputs");
+  elseif( ny >= nout )
+    error("is_dgkf: insufficient number of regulated outputs");
+  endif
+
+  nw = nin - nu;           nw1 = nw + 1;
+  nz = nout - ny;          nz1 = nz + 1;
+
+  [A,B,C,D] = sys2ss(Asys);
+  ## scale input/output for numerical reasons
+  if(norm (C, "fro") * norm (B, "fro") == 0)
+    error("||C||*||B|| = 0; no dynamic connnection from inputs to outputs");
+  endif
+  xx = sqrt(norm(B, Inf) / norm(C, Inf));
+  B = B / xx;  C = C * xx;
+
+  ## partition matrices
+                        Bw = B(:,1:nw);         Bu = B(:,nw1:nin);
+  Cz = C(1:nz,:);       Dzw = D(1:nz,1:nw);     Dzu = D(1:nz,nw1:nin);
+  Cy = C(nz1:nout,:);   Dyw = D(nz1:nout,1:nw); Dyu = D(nz1:nout,nw1:nin);
+
+  ## Check for loopo shifting
+  Dyu_nz = (norm(Dyu,Inf) != 0);
+  if (Dyu_nz)
+    warning("is_dgkf: D22 nonzero; performing loop shifting");
+  endif
+
+  ## 12 - rank condition at w = 0
+  xx =[A, Bu; Cz, Dzu];
+  [nr, nc] = size(xx);
+  irank = rank(xx);
+  if (irank != nc)
+    retval = 0;
+    warning(sprintf("rank([A Bu; Cz Dzu]) = %d, need %d; n=%d, nz=%d, nu=%d", ...
+        irank,nc,(Anc+Anz),nz,nu));
+    warning(" *** 12-rank condition violated at w = 0.");
+  endif
+
+  ## 21 - rank condition at w = 0
+  xx =[A, Bw; Cy, Dyw];
+  [nr, nc] = size(xx);
+  irank = rank(xx);
+  if (irank != nr)
+    retval = 0;
+    warning(sprintf("rank([A Bw; Cy Dyw]) = %d, need %d; n=%d, ny=%d, nw=%d", ...
+        irank,nr,(Anc+Anz),ny,nw));
+    warning(" *** 21-rank condition violated at w = 0.");
+  endif
+
+  ## can Dzu be transformed to become [0 I]' or [I]?
+  ## This ensures a normalized weight
+  [Qz, Ru] = qr(Dzu);
+  irank = rank(Ru);
+  if (irank != nu)
+    retval = 0;
+    warning(sprintf("*** rank(Dzu(%d x %d) = %d", nz, nu, irank));
+    warning(" *** Dzu does not have full column rank.");
+  endif
+  if (nu >= nz)
+    Qz = Qz(:,1:nu)';
+  else
+    Qz = [Qz(:,(nu+1):nz), Qz(:,1:nu)]';
+  endif
+  Ru = Ru(1:nu,:);
+
+  ## can Dyw be transformed to become [0 I] or [I]?
+  ## This ensures a normalized weight
+  [Qw, Ry] = qr(Dyw');
+  irank = rank(Ry);
+  if (irank != ny)
+    retval = 0;
+    warning(sprintf("*** rank(Dyw(%d x %d) = %d", ny, nw, irank));
+    warning(" *** Dyw does not have full row rank.");
+  endif
+
+  if (ny >= nw)
+    Qw = Qw(:,1:ny);
+  else
+    Qw = [Qw(:,(ny+1):nw), Qw(:,1:ny)];
+  endif
+  Ry = Ry(1:ny,:)';
+
+  ## transform P by Qz/Ru and Qw/Ry
+  Bw  = Bw*Qw;
+  Bu  = Bu/Ru;
+  B   = [Bw, Bu];
+  Cz  = Qz*Cz;
+  Cy  = Ry\Cy;
+  C   = [Cz; Cy];
+  Dzw = Qz*Dzw*Qw;
+  Dzu = Qz*Dzu/Ru;
+  Dyw = Ry\Dyw*Qw;
+
+  ## pack the return structure
+  dgkf_struct.nw        = nw;
+  dgkf_struct.nu        = nu;
+  dgkf_struct.nz        = nz;
+  dgkf_struct.ny        = ny;
+  dgkf_struct.A         = A;
+  dgkf_struct.Bw        = Bw;
+  dgkf_struct.Bu        = Bu;
+  dgkf_struct.Cz        = Cz;
+  dgkf_struct.Cy        = Cy;
+  dgkf_struct.Dzw       = Dzw;
+  dgkf_struct.Dzu       = Dzu;
+  dgkf_struct.Dyw       = Dyw;
+  dgkf_struct.Dyu       = Dyu;
+  dgkf_struct.Ru        = Ru;
+  dgkf_struct.Ry        = Ry;
+  dgkf_struct.Dyu_nz    = Dyu_nz;
+  dgkf_struct.dflg      = dflg;
+
+endfunction