Octave Control Systems Toolbox (OCST)
*************************************

1 Introduction
**************

The Octave Control Systems Toolbox (OCST) was initially developed by
Dr. A. Scottedward Hodel <a.s.hodel@eng.auburn.edu> with the assistance
of his students
   * R. Bruce Tenison <btenison@dibbs.net>,

   * David C. Clem,

   * John E. Ingram <John.Ingram@sea.siemans.com>, and

   * Kristi McGowan.
   This development was supported in part by NASA's Marshall Space
Flight Center as part of an in-house CACSD environment.  Additional
important contributions were made by Dr. Kai Mueller
<mueller@ifr.ing.tu-bs.de> and Jose Daniel Munoz Frias (`place.m').

   An on-line menu-driven tutorial is available via `DEMOcontrol';
beginning OCST users should start with this program.

 -- Function File:  DEMOcontrol
     Octave Control Systems Toolbox demo/tutorial program.  The demo
     allows the user to select among several categories of OCST
     function:
          octave:1> DEMOcontrol
          Octave Controls System Toolbox Demo

          [ 1] System representation
          [ 2] Block diagram manipulations
          [ 3] Frequency response functions
          [ 4] State space analysis functions
          [ 5] Root locus functions
          [ 6] LQG/H2/Hinfinity functions
          [ 7] End
     Command examples are interactively run for users to observe the use
     of OCST functions See also: bddemo, frdemo, analdemo, moddmeo,
     rldemo

2 System Data Structure
***********************

   The OCST stores all dynamic systems in a single data structure
format that can represent continuous systems, discrete-systems, and
mixed (hybrid) systems in state-space form, and can also represent
purely continuous/discrete systems in either transfer function or
pole-zero form. In order to provide more flexibility in treatment of
discrete/hybrid systems, the OCST also keeps a record of which system
outputs are sampled.

   Octave structures are accessed with a syntax much like that used by
the C programming language.  For consistency in use of the data
structure used in the OCST, it is recommended that the system structure
access m-files be used (see sysinterface).  Some elements of the
data structure are absent depending on the internal system
representation(s) used.  More than one system representation can be
used for SISO systems; the OCST m-files ensure that all representations
used are consistent with one another.

 -- Function File:  sysrepdemo
     Tutorial for the use of the system data structure functions

2.1 Variables common to all OCST system formats
===============================================

The data structure elements (and variable types) common to all  system
representations are listed below; examples of the initialization and
use of the system data structures are given in subsequent sections and
in the online demo `DEMOcontrol'.
N
NZ
     The respective number of continuous and discrete states in the
     system (scalar)

INNAME
OUTNAME
     list of name(s) of the system input, output signal(s). (list of
     strings)

SYS
     System status vector.  (vector)

     This vector indicates both what representation was used to
     initialize the system data structure (called the primary system
     type) and which other representations are currently up-to-date
     with the primary system type (see structaccess).

     The value of the first element of the vector indicates the primary
     system type.

    0
          for tf form (initialized with `tf2sys' or `fir2sys')

    1
          for zp form (initialized with `zp2sys')

    2
          for ss form (initialized with `ss2sys')

     The next three elements are boolean flags that indicate whether
     tf, zp, or ss, respectively, are "up to date" (whether it is safe
     to use the variables associated with these representations).
     These flags are changed when calls are made to the `sysupdate'
     command.

TSAM
     Discrete time sampling period  (nonnegative scalar).   TSAM is set
     to 0 for continuous time systems.

YD
     Discrete-time output list (vector)

     indicates which outputs are discrete time (i.e.,     produced by
     D/A converters) and which are continuous time.      yd(ii) = 0 if
     output ii is continuous, = 1 if discrete.

   The remaining variables of the  system data structure are only
present if the corresponding entry of the `sys' vector is true (=1).

2.2 `tf' format variables
=========================

NUM
     numerator coefficients   (vector)

DEN
     denominator coefficients   (vector)

2.3 `zp' format variables
=========================

ZER
     system zeros   (vector)

POL
     system poles    (vector)

K
     leading coefficient   (scalar)

2.4 `ss' format variables
=========================

A
B
C
D
     The usual state-space matrices. If a system has both
     continuous and discrete states, they are sorted so that
     continuous states come first, then discrete states

     *Note* some functions (e.g., `bode', `hinfsyn') will not accept
     systems with both discrete and continuous states/outputs

STNAME
     names of system states   (list of strings)

3 System Construction and Interface Functions
*********************************************

Construction and manipulations of the OCST system data structure (see
sysstruct) requires attention to many details in order to ensure that
data structure contents remain consistent.  Users are strongly
encouraged to use the system interface functions in this section.
Functions for the formatted display in of system data structures are
given in sysdisp.

3.1 Finite impulse response system interface functions
======================================================

 -- Function File:  fir2sys (NUM, TSAM, INNAME, OUTNAME)
     construct a system data structure from FIR description

     *Inputs*
    NUM
          vector of coefficients [c0, c1, ..., cn] of the SISO FIR
          transfer function C(z) = c0 + c1*z^(-1) + c2*z^(-2) + ... +
          cn*z^(-n)

    TSAM
          sampling time (default: 1)

    INNAME
          name of input signal;  may be a string or a list with a
          single entry

    OUTNAME
          name of output signal; may be a string or a list with a
          single entry

     *Output*
    SYS
          system data structure

     *Example*
          octave:1> sys = fir2sys([1 -1 2 4],0.342,\
          > "A/D input","filter output");
          octave:2> sysout(sys)
          Input(s)
          1: A/D input

          Output(s):
          1: filter output (discrete)

          Sampling interval: 0.342
          transfer function form:
          1*z^3 - 1*z^2 + 2*z^1 + 4
          -------------------------
          1*z^3 + 0*z^2 + 0*z^1 + 0

 -- Function File: [C, TSAM, INPUT, OUTPUT] = sys2fir (SYS)
     Extract FIR data from system data structure; see `fir2sys' for
     parameter descriptions See also: fir2sys

3.2 State space system interface functions
==========================================

 -- Function File: OUTSYS = ss (A, B, C, D, TSAM, N, NZ, STNAME,
          INNAME, OUTNAME, OUTLIST)
     Create system structure from state-space data.   May be continuous,
     discrete, or mixed (sampled data)

     *Inputs*
    A
    B
    C
    D
          usual state space matrices

          default: D = zero matrix

    TSAM
          sampling rate.  Default: tsam = 0 (continuous system)

    N
    NZ
          number of continuous, discrete states in the system

          If TSAM is 0, n = `rows'(A), nz = 0

          If TSAM is greater than zero, n = 0, nz = `rows'(A)

          see below for system partitioning

    STNAME
          cell array of strings of state signal names

          default (STNAME=[] on input): `x_n' for continuous states,
          `xd_n' for discrete states

    INNAME
          cell array of strings of input signal names

          default (INNAME = [] on input): `u_n'

    OUTNAME
          cell array of strings of output signal names

          default (OUTNAME = [] on input): `y_n'

    OUTLIST
          list of indices of outputs y that are sampled

          If TSAM is 0, outlist = []

          If TSAM is greater than 0, outlist = 1:`rows'(C)

     Unlike states, discrete/continuous outputs may appear in any order

     `sys2ss' returns a vector YD where YD(OUTLIST) = 1; all other
     entries of YD are 0

     *Output*
    OUTSYS
          system data structure

     *System partitioning*

     Suppose for simplicity that outlist specified that the first
     several outputs were continuous and the remaining outputs were
     discrete.  Then the system is partitioned as
          x = [ xc ]  (n x 1)
          [ xd ]  (nz x 1 discrete states)
          a = [ acc acd ]  b = [ bc ]
          [ adc add ]      [ bd ]
          c = [ ccc ccd ]  d = [ dc ]
          [ cdc cdd ]      [ dd ]

          (cdc = c(outlist,1:n), etc.)
     with dynamic equations: d/dt xc(t)     = acc*xc(t)      +
     acd*xd(k*tsam) + bc*u(t)

     xd((k+1)*tsam) = adc*xc(k*tsam) + add*xd(k*tsam) + bd*u(k*tsam)

     yc(t)      = ccc*xc(t)      + ccd*xd(k*tsam) + dc*u(t)

     yd(k*tsam) = cdc*xc(k*tsam) + cdd*xd(k*tsam) + dd*u(k*tsam)

     *Signal partitions*
          | continuous      | discrete               |
          ----------------------------------------------------
          states  | stname(1:n,:)   | stname((n+1):(n+nz),:) |
          ----------------------------------------------------
          outputs | outname(cout,:) | outname(outlist,:)     |
          ----------------------------------------------------
     where cout is the list of in 1:`rows'(P) that are not contained in
     outlist. (Discrete/continuous outputs may be entered in any order
     desired by the user.)

     *Example*
          octave:1> a = [1 2 3; 4 5 6; 7 8 10];
          octave:2> b = [0 0 ; 0 1 ; 1 0];
          octave:3> c = eye (3);
          octave:4> sys = ss (a, b, c, [], 0, 3, 0, ..
          >                   {"volts", "amps", "joules"});
          octave:5> sysout(sys);
          Input(s)
          1: u_1
          2: u_2

          Output(s):
          1: y_1
          2: y_2
          3: y_3

          state-space form:
          3 continuous states, 0 discrete states
          State(s):
          1: volts
          2: amps
          3: joules

          A matrix: 3 x 3
          1   2   3
          4   5   6
          7   8  10
          B matrix: 3 x 2
          0  0
          0  1
          1  0
          C matrix: 3 x 3
          1  0  0
          0  1  0
          0  0  1
          D matrix: 3 x 3
          0  0
          0  0
          0  0
     Notice that the D matrix is constructed  by default to the correct
     dimensions.  Default input and output signals names were assigned
     since none were given

 -- Function File:  ss2sys (A, B, C, D, TSAM, N, NZ, STNAME, INNAME,
          OUTNAME, OUTLIST)
     Create system structure from state-space data.   May be continuous,
     discrete, or mixed (sampled data)

     *Inputs*
    A
    B
    C
    D
          usual state space matrices

          default: D = zero matrix

    TSAM
          sampling rate.  Default: tsam = 0 (continuous system)

    N
    NZ
          number of continuous, discrete states in the system

          If TSAM is 0, n = `rows'(A), nz = 0

          If TSAM is greater than zero, n = 0, nz = `rows'(A)

          see below for system partitioning

    STNAME
          cell array of strings of state signal names

          default (STNAME=[] on input): `x_n' for continuous states,
          `xd_n' for discrete states

    INNAME
          cell array of strings of input signal names

          default (INNAME = [] on input): `u_n'

    OUTNAME
          cell array of strings of input signal names

          default (OUTNAME = [] on input): `y_n'

    OUTLIST
          list of indices of outputs y that are sampled

          If TSAM is 0, outlist = []

          If TSAM is greater than 0, outlist = 1:`rows'(C)

     Unlike states, discrete/continuous outputs may appear in any order

     `sys2ss' returns a vector YD where YD(OUTLIST) = 1; all other
     entries of YD are 0

     *Outputs* OUTSYS = system data structure

     *System partitioning*

     Suppose for simplicity that outlist specified that the first
     several outputs were continuous and the remaining outputs were
     discrete.  Then the system is partitioned as
          x = [ xc ]  (n x 1)
          [ xd ]  (nz x 1 discrete states)
          a = [ acc acd ]  b = [ bc ]
          [ adc add ]      [ bd ]
          c = [ ccc ccd ]  d = [ dc ]
          [ cdc cdd ]      [ dd ]

          (cdc = c(outlist,1:n), etc.)
     with dynamic equations: d/dt xc(t)     = acc*xc(t)      +
     acd*xd(k*tsam) + bc*u(t)

     xd((k+1)*tsam) = adc*xc(k*tsam) + add*xd(k*tsam) + bd*u(k*tsam)

     yc(t)      = ccc*xc(t)      + ccd*xd(k*tsam) + dc*u(t)

     yd(k*tsam) = cdc*xc(k*tsam) + cdd*xd(k*tsam) + dd*u(k*tsam)

     *Signal partitions*
          | continuous      | discrete               |
          ----------------------------------------------------
          states  | stname(1:n,:)   | stname((n+1):(n+nz),:) |
          ----------------------------------------------------
          outputs | outname(cout,:) | outname(outlist,:)     |
          ----------------------------------------------------
     where cout is the list of in 1:`rows'(P) that are not contained in
     outlist. (Discrete/continuous outputs may be entered in any order
     desired by the user.)

     *Example*
          octave:1> a = [1 2 3; 4 5 6; 7 8 10];
          octave:2> b = [0 0 ; 0 1 ; 1 0];
          octave:3> c = eye (3);
          octave:4> sys = ss (a, b, c, [], 0, 3, 0,
          >                   {"volts", "amps", "joules"});
          octave:5> sysout(sys);
          Input(s)
          1: u_1
          2: u_2

          Output(s):
          1: y_1
          2: y_2
          3: y_3

          state-space form:
          3 continuous states, 0 discrete states
          State(s):
          1: volts
          2: amps
          3: joules

          A matrix: 3 x 3
          1   2   3
          4   5   6
          7   8  10
          B matrix: 3 x 2
          0  0
          0  1
          1  0
          C matrix: 3 x 3
          1  0  0
          0  1  0
          0  0  1
          D matrix: 3 x 3
          0  0
          0  0
          0  0
     Notice that the D matrix is constructed  by default to the correct
     dimensions.  Default input and output signals names were assigned
     since none were given

 -- Function File: [A, B, C, D, TSAM, N, NZ, STNAME, INNAME, OUTNAME,
YD] = sys2ss (SYS)
     Extract state space representation from system data structure

     *Input*
    SYS
          System data structure

     *Outputs*
    A
    B
    C
    D
          State space matrices for SYS

    TSAM
          Sampling time of SYS (0 if continuous)

    N
    NZ
          Number of continuous, discrete states (discrete states come
          last in state vector X)

    STNAME
    INNAME
    OUTNAME
          Signal names (lists of strings);  names of states, inputs,
          and outputs, respectively

    YD
          Binary vector; YD(II) is 1 if output Y(II) is discrete
          (sampled); otherwise  YD(II) is 0

     A warning massage is printed if the system is a mixed continuous
     and discrete system

     *Example*
          octave:1> sys=tf2sys([1 2],[3 4 5]);
          octave:2> [a,b,c,d] = sys2ss(sys)
          a =
          0.00000   1.00000
          -1.66667  -1.33333
          b =
          0
          1
          c = 0.66667  0.33333
          d = 0

3.3 Transfer function system interface functions
================================================

 -- Function File:  tf (NUM, DEN, TSAM, INNAME, OUTNAME)
     build system data structure from transfer function format data

     *Inputs*
    NUM
    DEN
          coefficients of numerator/denominator polynomials

    TSAM
          sampling interval. default: 0 (continuous time)

    INNAME
    OUTNAME
          input/output signal names; may be a string or cell array with
          a single string entry

     *Outputs* SYS = system data structure

     *Example*
          octave:1> sys=tf([2 1],[1 2 1],0.1);
          octave:2> sysout(sys)
          Input(s)
          1: u_1
          Output(s):
          1: y_1 (discrete)
          Sampling interval: 0.1
          transfer function form:
          2*z^1 + 1
          -----------------
          1*z^2 + 2*z^1 + 1

 -- Function File:  tf2sys (NUM, DEN, TSAM, INNAME, OUTNAME)
     Build system data structure from transfer function format data

     *Inputs*
    NUM
    DEN
          Coefficients of numerator/denominator polynomials

    TSAM
          Sampling interval; default: 0 (continuous time)

    INNAME
    OUTNAME
          Input/output signal names; may be a string or cell array with
          a single string entry

     *Output*
    SYS
          System data structure

     *Example*
          octave:1> sys=tf2sys([2 1],[1 2 1],0.1);
          octave:2> sysout(sys)
          Input(s)
          1: u_1
          Output(s):
          1: y_1 (discrete)
          Sampling interval: 0.1
          transfer function form:
          2*z^1 + 1
          -----------------
          1*z^2 + 2*z^1 + 1

 -- Function File: [NUM, DEN, TSAM, INNAME, OUTNAME] = sys2tf (SYS)
     Extract transfer function data from a system data structure

     See `tf' for parameter descriptions

     *Example*
          octave:1> sys=ss([1 -2; -1.1,-2.1],[0;1],[1 1]);
          octave:2> [num,den] = sys2tf(sys)
          num = 1.0000  -3.0000
          den = 1.0000   1.1000  -4.3000

3.4 Zero-pole system interface functions
========================================

 -- Function File:  zp (ZER, POL, K, TSAM, INNAME, OUTNAME)
     Create system data structure from zero-pole data

     *Inputs*
    ZER
          vector of system zeros

    POL
          vector of system poles

    K
          scalar leading coefficient

    TSAM
          sampling period. default: 0 (continuous system)

    INNAME
    OUTNAME
          input/output signal names (lists of strings)

     *Outputs* sys: system data structure

     *Example*
          octave:1> sys=zp([1 -1],[-2 -2 0],1);
          octave:2> sysout(sys)
          Input(s)
          1: u_1
          Output(s):
          1: y_1
          zero-pole form:
          1 (s - 1) (s + 1)
          -----------------
          s (s + 2) (s + 2)

 -- Function File:  zp2sys (ZER, POL, K, TSAM, INNAME, OUTNAME)
     Create system data structure from zero-pole data

     *Inputs*
    ZER
          Vector of system zeros

    POL
          Vector of system poles

    K
          Scalar leading coefficient

    TSAM
          Sampling period; default: 0 (continuous system)

    INNAME
    OUTNAME
          Input/output signal names (lists of strings)

     *Output*
    SYS
          System data structure

     *Example*
          octave:1> sys=zp2sys([1 -1],[-2 -2 0],1);
          octave:2> sysout(sys)
          Input(s)
          1: u_1
          Output(s):
          1: y_1
          zero-pole form:
          1 (s - 1) (s + 1)
          -----------------
          s (s + 2) (s + 2)

 -- Function File: [ZER, POL, K, TSAM, INNAME, OUTNAME] = sys2zp (SYS)
     Extract zero/pole/leading coefficient information from a system
     data structure

     See `zp' for parameter descriptions

     *Example*
          octave:1> sys=ss([1 -2; -1.1,-2.1],[0;1],[1 1]);
          octave:2> [zer,pol,k] = sys2zp(sys)
          zer = 3.0000
          pol =
          -2.6953
          1.5953
          k = 1

3.5 Data structure access functions
===================================

 -- Function File:  syschnames (SYS, OPT, LIST, NAMES)
     Superseded by `syssetsignals'

 -- Function File:  syschtsam (SYS, TSAM)
     This function changes the sampling time (tsam) of the system.
     Exits with an error if sys is purely continuous time

 -- Function File: [N, NZ, M, P, YD] = sysdimensions (SYS, OPT)
     return the number of states, inputs, and/or outputs in the system
     SYS

     *Inputs*
    SYS
          system data structure

    OPT
          String indicating which dimensions are desired.  Values:
         `"all"'
               (default) return all parameters as specified under
               Outputs below

         `"cst"'
               return N= number of continuous states

         `"dst"'
               return N= number of discrete states

         `"in"'
               return N= number of inputs

         `"out"'
               return N= number of outputs

     *Outputs*
    N
          number of continuous states (or individual requested
          dimension as specified by OPT)

    NZ
          number of discrete states

    M
          number of system inputs

    P
          number of system outputs

    YD
          binary vector; YD(II) is nonzero if output II is discrete
          yd(ii) = 0 if output II is continuous
     See also: sysgetsignals, sysgettsam

 -- Function File: [STNAME, INNAME, OUTNAME, YD] = sysgetsignals (SYS)
 -- Function File: SIGLIST = sysgetsignals (SYS, SIGID)
 -- Function File: SIGNAME = sysgetsignals (SYS, SIGID, SIGNUM, STRFLG)
     Get signal names from a system

     *Inputs*
    SYS
          system data structure for the state space system

    SIGID
          signal id.  String.  Must be one of
         `"in"'
               input signals

         `"out"'
               output signals

         `"st"'
               stage signals

         `"yd"'
               value of logical vector YD

    SIGNUM
          index(indices) or name(s) or signals; see `sysidx'

    STRFLG
          flag to return a string instead of a cell array;  Values:
         `0'
               (default) return a cell array (even if signum specifies
               an individual signal)

         `1'
               return a string.  Exits with an error if signum does not
               specify an individual signal

     *Outputs*
    *If SIGID is not specified:

         STNAME
         INNAME
         OUTNAME
               signal names (cell array of strings);  names of states,
               inputs, and outputs, respectively

         YD
               binary vector; YD(II) is nonzero if output II is discrete

    *If SIGID is specified but SIGNUM is not specified:

         `sigid="in"'
               SIGLIST is set to the cell array of input names

         `sigid="out"'
               SIGLIST is set to the cell array of output names

         `sigid="st"'
               SIGLIST is set to the cell array of state names

               stage signals

         `sigid="yd"'
               SIGLIST is set to logical vector indicating discrete
               outputs; SIGLIST(II) = 0 indicates that output II is
               continuous (unsampled), otherwise it is discrete

    *If the first three input arguments are specified:
          SIGNAME is a cell array of the specified signal names (SIGID
          is `"in"', `"out"', or `"st"'), or else the logical flag
          indicating whether output(s) SIGNUM is(are) discrete
          (SIGVAL=1) or continuous (SIGVAL=0)

     *Examples* (From `sysrepdemo')
          octave> sys=ss(rand(4),rand(4,2),rand(3,4));
          octave># get all signal names
          octave> [Ast,Ain,Aout,Ayd] = sysgetsignals(sys)
          Ast =
          (
          [1] = x_1
          [2] = x_2
          [3] = x_3
          [4] = x_4
          )
          Ain =
          (
          [1] = u_1
          [2] = u_2
          )
          Aout =
          (
          [1] = y_1
          [2] = y_2
          [3] = y_3
          )
          Ayd =

          0  0  0
          octave> # get only input signal names:
          octave> Ain = sysgetsignals(sys,"in")
          Ain =
          (
          [1] = u_1
          [2] = u_2
          )
          octave> # get name of output 2 (in cell array):
          octave> Aout = sysgetsignals(sys,"out",2)
          Aout =
          (
          [1] = y_2
          )
          octave> # get name of output 2 (as string):
          octave> Aout = sysgetsignals(sys,"out",2,1)
          Aout = y_2

 -- Function File:  sysgettype (SYS)
     return the initial system type of the system

     *Input*
    SYS
          System data structure

     *Output*
    SYSTYPE
          String indicating how the structure was initially
          constructed. Values: `"ss"', `"zp"', or `"tf"'

     FIR initialized systems return `systype="tf"'

 -- Function File:  syssetsignals (SYS, OPT, NAMES, SIG_IDX)
     change the names of selected inputs, outputs and states

     *Inputs*
    SYS
          System data structure

    OPT
          Change default name (output)

         `"out"'
               Change selected output names

         `"in"'
               Change selected input names

         `"st"'
               Change selected state names

         `"yd"'
               Change selected outputs from discrete to continuous or
               from continuous to discrete

    NAMES

         `opt = "out", "in", "st"'
               string or string array containing desired signal names
               or values

         `opt = "yd"'
               To desired output continuous/discrete flag Set name to 0
               for continuous, or 1 for discrete

    SIG_IDX
          indices or names of outputs, yd, inputs, or states whose
          respective names/values should be changed

          Default: replace entire cell array of names/entire yd vector

     *Outputs*
    RETSYS
          SYS with appropriate signal names changed (or YD values,
          where appropriate)

     *Example*
          octave:1> sys=ss ([1 2; 3 4],[5;6],[7 8]);
          octave:2> sys = syssetsignals (sys, "st",
          >                              str2mat("Posx","Velx"));
          octave:3> sysout(sys)
          Input(s)
          1: u_1
          Output(s):
          1: y_1
          state-space form:
          2 continuous states, 0 discrete states
          State(s):
          1: Posx
          2: Velx
          A matrix: 2 x 2
          1  2
          3  4
          B matrix: 2 x 1
          5
          6
          C matrix: 1 x 2
          7  8
          D matrix: 1 x 1
          0

 -- Function File:  sysupdate (SYS, OPT)
     Update the internal representation of a system

     *Inputs*
    SYS:
          system data structure

    OPT
          string:
         `"tf"'
               update transfer function form

         `"zp"'
               update zero-pole form

         `"ss"'
               update state space form

         `"all"'
               all of the above

     *Outputs*
    RETSYS
          Contains union of data in sys and requested data If requested
          data in SYS is already up to date then RETSYS=SYS

     Conversion to `tf' or `zp' exits with an error if the system is
     mixed continuous/digital See also: tf, ss, zp, sysout, sys2ss,
     sys2tf, sys2zp

 -- Function File: [SYSTYPE, NOUT, NIN, NCSTATES, NDSTATES] = minfo
          (INMAT)
     Determines the type of system matrix.  INMAT can be a varying, a
     system, a constant, and an empty matrix

     *Outputs*
    SYSTYPE
          Can be one of: varying, system, constant, and empty

    NOUT
          The number of outputs of the system

    NIN
          The number of inputs of the system

    NCSTATES
          The number of continuous states of the system

    NDSTATES
          The number of discrete states of the system

 -- Function File:  sysgettsam (SYS)
     Return the sampling time of the system SYS

4 System display functions
**************************

 -- Function File:  sysout (SYS, OPT)
     print out a system data structure in desired format
    SYS
          system data structure

    OPT
          Display option
         `[]'
               primary system form (default)

         `"ss"'
               state space form

         `"tf"'
               transfer function form

         `"zp"'
               zero-pole form

         `"all"'
               all of the above

 -- Function File:  tfout (NUM, DENOM, X)
     Print formatted transfer function n(s)/d(s) to the screen X
     defaults to the string `"s"' See also: polyval, polyvalm, poly,
     roots, conv, deconv, residue, filter, polyderiv, polyinteg, polyout

 -- Function File:  zpout (ZER, POL, K, X)
     print formatted zero-pole form to the screen X defaults to the
     string `"s"' See also: polyval, polyvalm, poly, roots, conv,
     deconv, residue, filter, polyderiv, polyinteg, polyout

5 Block Diagram Manipulations
*****************************

See systime.

   Unless otherwise noted, all parameters (input,output) are system
data structures.

 -- Function File:  bddemo (INPUTS)
     Octave Controls toolbox demo: Block Diagram Manipulations demo

 -- Function File:  buildssic (CLST, ULST, OLST, ILST, S1, S2, S3, S4,
          S5, S6, S7, S8)
     Form an arbitrary complex (open or closed loop) system in
     state-space form from several systems. `buildssic' can easily
     (despite its cryptic syntax) integrate transfer functions from a
     complex block diagram into a single system with one call This
     function is especially useful for building open loop
     interconnections for H-infinity and H-2 designs or for closing
     loops with these controllers

     Although this function is general purpose, the use of `sysgroup'
     `sysmult', `sysconnect' and the like is recommended for standard
     operations since they can handle mixed discrete and continuous
     systems and also the names of inputs, outputs, and states

     The parameters consist of 4 lists that describe the connections
     outputs and inputs and up to 8 systems S1-S8 Format of the lists:
    CLST
          connection list, describes the input signal of each system.
          The maximum number of rows of Clst is equal to the sum of all
          inputs of s1-s8

          Example: `[1 2 -1; 2 1 0]' means that:  new input 1 is old
          input 1 + output 2 - output 1, and new input 2 is old input 2
          + output 1. The order of rows is arbitrary

    ULST
          if not empty the old inputs in vector ULST will be appended
          to the outputs. You need this if you want to "pull out" the
          input of a system. Elements are input numbers of S1-S8

    OLST
          output list, specifies the outputs of the resulting systems.
          Elements are output numbers of S1-S8 The numbers are allowed
          to be negative and may appear in any order. An empty matrix
          means all outputs

    ILST
          input list, specifies the inputs of the resulting systems.
          Elements are input numbers of S1-S8 The numbers are allowed
          to be negative and may appear in any order. An empty matrix
          means all inputs

     Example:  Very simple closed loop system
          w        e  +-----+   u  +-----+
          --->o--*-->|  K  |--*-->|  G  |--*---> y
          ^  |   +-----+  |   +-----+  |
          - |  |            |            |
          |  |            +----------------> u
          |  |                         |
          |  +-------------------------|---> e
          |                            |
          +----------------------------+

     The closed loop system GW can be obtained by
          GW = buildssic([1 2; 2 -1], 2, [1 2 3], 2, G, K);

    CLST
          1st row: connect input 1 (G) with output 2 (K)

          2nd row: connect input 2 (K) with negative output 1 (G)

    ULST
          Append input of 2 (K) to the number of outputs

    OLST
          Outputs are output of 1 (G), 2 (K) and appended output 3
          (from ULST)

    ILST
          The only input is 2 (K)

     Here is a real example:
          +----+
          -------------------->| W1 |---> v1
          z   |                    +----+
          ----|-------------+
          |             |
          |    +---+    v      +----+
          *--->| G |--->O--*-->| W2 |---> v2
          |    +---+       |   +----+
          |                |
          |                v
          u                  y

          min || GW   ||
          vz   infty

     The closed loop system GW from [z, u]' to [v1, v2, y]' can be
     obtained by (all SISO systems):
          GW = buildssic([1, 4; 2, 4; 3, 1], 3, [2, 3, 5],
          [3, 4], G, W1, W2, One);
     where "One" is a unity gain (auxiliary) function with order 0
     (e.g. `One = ugain(1);')

 -- Function File: SYS = jet707 ()
     Creates a linearized state-space model of a Boeing 707-321 aircraft
     at V=80 m/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

     *Reference*: R. Brockhaus: `Flugregelung' (Flight Control),
     Springer, 1994 See also: ord2

 -- Function File:  ord2 (NFREQ, DAMP, GAIN)
     Creates a continuous 2nd order system with parameters:

     *Inputs*
    NFREQ
          natural frequency [Hz]. (not in rad/s)

    DAMP
          damping coefficient

    GAIN
          dc-gain This is steady state value only for damp > 0 gain is
          assumed to be 1.0 if omitted

     *Output*
    OUTSYS
          system data structure has representation with w = 2 * pi *
          nfreq:
               /                                        \
               | / -2w*damp -w \  / w \                 |
               G = | |             |, |   |, [ 0  gain ], 0 |
               | \   w       0 /  \ 0 /                 |
               \                                        /
     *See also* `jet707' (MIMO example, Boeing 707-321 aircraft model)

 -- Function File:  sysadd (GSYS, HSYS)
     returns SYS = GSYS + HSYS
        * Exits with an error if GSYS and HSYS are not compatibly
          dimensioned

        * Prints a warning message is system states have identical
          names; duplicate names are given a suffix to make them unique

        * SYS input/output names are taken from GSYS

          ________
          ----|  gsys  |---
          u   |    ----------  +|
          -----                (_)----> y
          |     ________   +|
          ----|  hsys  |---
          --------

 -- Function File: SYS = sysappend (SYST, B, C, D, OUTNAME, INNAME, YD)
     appends new inputs and/or outputs to a system

     *Inputs*
    SYST
          system data structure

    B
          matrix to be appended to sys "B" matrix (empty if none)

    C
          matrix to be appended to sys "C" matrix (empty if none)

    D
          revised sys d matrix (can be passed as [] if the revised d is
          all zeros)

    OUTNAME
          list of names for new outputs

    INNAME
          list of names for new inputs

    YD
          binary vector; yd(ii)=0 indicates a continuous output;
          yd(ii)=1 indicates a discrete output

     *Outputs*
    SYS
               sys.b := [syst.b , b]
               sys.c := [syst.c  ]
               [ c     ]
               sys.d := [syst.d | D12 ]
               [ D21   | D22 ]
          where D12, D21, and D22 are the appropriate dimensioned
          blocks of the input parameter D
             * The leading block D11 of D is ignored

             * If INNAME and OUTNAME are not given as arguments, the
               new inputs and outputs are be assigned default names

             * YD is a binary vector of length rows(c) that indicates
               continuous/sampled outputs.  Default value for YD is:
                  - SYS is continuous or mixed YD = `zeros(1,rows(c))'

                  - SYS is discrete YD = `ones(1,rows(c))'

 -- Function File: CLSYS = sysconnect (SYS, OUT_IDX, IN_IDX, ORDER, TOL)
     Close the loop from specified outputs to respective specified
     inputs

     *Inputs*
    SYS
          System data structure

    OUT_IDX
    IN_IDX
          Names or indices of signals to connect (see `sysidx') The
          output specified by out_idx(ii) is connected to the input
          specified by in_idx(ii)

    ORDER
          logical flag (default = 0)
         `0'
               Leave inputs and outputs in their original order

         `1'
               Permute inputs and outputs to the order shown in the
               diagram below

    TOL
          Tolerance for singularities in algebraic loops, default:
          200`eps'

     *Outputs*
    CLSYS
          Resulting closed loop system

     *Method*

     `sysconnect' internally permutes selected inputs, outputs as shown
     below, closes the loop, and then permutes inputs and outputs back
     to their original order
          --------------------
          u_1       ----->|                  |----> y_1
          |        sys       |
          old u_2 |                  |
          u_2* ---->(+)--->|                  |----->y_2
          (in_idx)   ^     --------------------    | (out_idx)
          |                             |
          -------------------------------
     The input that has the summing junction added to it has an * added
     to the end  of the input name

 -- Function File: [CSYS, ACD, CCD] = syscont (SYS)
     Extract the purely continuous subsystem of an input system

     *Input*
    SYS
          system data structure

     *Outputs*
    CSYS
          is the purely continuous input/output connections of SYS

    ACD
    CCD
          connections from discrete states to continuous states,
          discrete states to continuous outputs, respectively

          If no continuous path exists, CSYS will be empty

 -- Function File: [DSYS, ADC, CDC] = sysdisc (SYS)
     *Input*
    SYS
          System data structure

     *Outputs*
    DSYS
          Purely discrete portion of sys (returned empty if there is no
          purely discrete path from inputs to outputs)

    ADC
    CDC
          Connections from continuous states to discrete states and
          discrete outputs, respectively

 -- Function File: RETSYS = sysdup (ASYS, OUT_IDX, IN_IDX)
     Duplicate specified input/output connections of a system

     *Inputs*
    ASYS
          system data structure

    OUT_IDX
    IN_IDX
          indices or names of desired signals (see `sigidx') duplicates
          are made of `y(out_idx(ii))' and `u(in_idx(ii))'

     *Output*
    RETSYS
          Resulting closed loop system: duplicated i/o names are
          appended with a `"+"' suffix

     *Method*

     `sysdup' creates copies of selected inputs and outputs as shown
     below.  U1, Y1 is the set of original inputs/outputs, and U2, Y2
     is the set of duplicated inputs/outputs in the order specified in
     IN_IDX, OUT_IDX, respectively
          ____________________
          u1  ----->|                  |----> y1
          |       asys       |
          u2 ------>|                  |----->y2
          (in_idx)  -------------------- (out_idx)

 -- Function File: SYS = sysgroup (ASYS, BSYS)
     Combines two systems into a single system

     *Inputs*
    ASYS
    BSYS
          System data structures

     *Output*
    SYS
          sys = block diag(asys,bsys)

          __________________
          |    ________    |
          u1 ----->|--> | asys |--->|----> y1
          |    --------    |
          |    ________    |
          u2 ----->|--> | bsys |--->|----> y2
          |    --------    |
          ------------------
          Ksys
     The function also rearranges the internal state-space realization
     of SYS so that the continuous states come first and the discrete
     states come last If there are duplicate names, the second name has
     a unique suffix appended on to the end of the name

 -- Function File: SYS = sysmult (ASYS, BSYS)
     Compute sys = Asys*Bsys (series connection):
          u   ----------     ----------
          --->|  Bsys  |---->|  Asys  |--->
          ----------     ----------
     A warning occurs if there is direct feed-through from an input or
     a continuous state of BSYS, through a discrete output of BSYS, to
     a continuous state or output in ASYS (system data structure does
     not recognize discrete inputs)

 -- Function File: RETSYS = sysprune (ASYS, OUT_IDX, IN_IDX)
     Extract specified inputs/outputs from a system

     *Inputs*
    ASYS
          system data structure

    OUT_IDX
    IN_IDX
          Indices or signal names of the outputs and inputs to be kept
          in the returned system; remaining connections are "pruned" off
          May select as [] (empty matrix) to specify all outputs/inputs

               retsys = sysprune (Asys, [1:3,4], "u_1");
               retsys = sysprune (Asys, {"tx", "ty", "tz"}, 4);

     *Output*
    RETSYS
          Resulting system

          ____________________
          u1 ------->|                  |----> y1
          (in_idx)  |       Asys       | (out_idx)
          u2 ------->|                  |----| y2
          (deleted)-------------------- (deleted)

 -- Function File: PV = sysreorder (VLEN, LIST)
     *Inputs*
    VLEN
          Vector length

    LIST
          A subset of `[1:vlen]'

     *Output*
    PV
          A permutation vector to order elements of `[1:vlen]' in
          `list' to the end of a vector

     Used internally by `sysconnect' to permute vector elements to their
     desired locations

 -- Function File: RETSYS = sysscale (SYS, OUTSCALE, INSCALE, OUTNAME,
          INNAME)
     scale inputs/outputs of a system

     *Inputs*
    SYS
          Structured system

    OUTSCALE
    INSCALE
          Constant matrices of appropriate dimension

    OUTNAME
    INNAME
          Lists of strings with the names of respectively outputs and
          inputs

     *Output*
    RETSYS
          resulting open loop system:
               -----------    -------    -----------
               u --->| inscale |--->| sys |--->| outscale |---> y
               -----------    -------    -----------
     If the input names and output names (each a list of strings) are
     not given and the scaling matrices are not square, then default
     names will be given to the inputs and/or outputs

     A warning message is printed if outscale attempts to add continuous
     system outputs to discrete system outputs; otherwise YD is set
     appropriately in the returned value of SYS

 -- Function File: SYS = syssub (GSYS, HSYS)
     Return sys = Gsys - Hsys

     *Method*

     GSYS and HSYS are connected in parallel The input vector is
     connected to both systems; the outputs are subtracted.  Returned
     system names are those of GSYS
          +--------+
          +--->|  Gsys  |---+
          |    +--------+   |
          |                +|
          u --+                (_)--> y
          |                -|
          |    +--------+   |
          +--->|  Hsys  |---+
          +--------+

 -- Function File:  ugain (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 See also:
     hinfdemo, jet707

 -- Function File: W = wgt1o (VL, VH, FC)
     State space description of a first order weighting function

     Weighting function are needed by the H-2/H-infinity design
     procedure These functions are part of the augmented plant P (see
     `hinfdemo' for an application example)

     *Inputs*
    VL
          Gain at low frequencies

    VH
          Gain at high frequencies

    FC
          Corner frequency (in Hz, *not* in rad/sec)

     *Output*
    W
          Weighting function, given in form of a system data structure

 -- Function File: KSYS = parallel (ASYS, BSYS)
     Forms the parallel connection of two systems

          --------------------
          |      --------    |
          u  ----->|----> | asys |--->|----> y1
          |    |      --------    |
          |    |      --------    |
          |--->|----> | bsys |--->|----> y2
          |      --------    |
          --------------------
          ksys

 -- Function File: [RETSYS, NC, NO] = sysmin (SYS, FLG)
     Returns a minimal (or reduced order) system

     *Inputs*
    SYS
          System data structure

    FLG
          When equal to 0 (default value), returns minimal system, in
          which state names are lost; when equal to 1, returns 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)'

6 Numerical Functions
*********************

 -- Function File: X = are (A, B, C, OPT)
     Solve the Algebraic Riccati Equation
          a' * x + x * a - x * b * x + c = 0

     *Inputs* for identically dimensioned square matrices
    A
          N by N matrix;

    B
          N by N matrix or N by M matrix; in the latter case B is
          replaced by b:=b*b';

    C
          N by N matrix or P by M matrix; in the latter case C is
          replaced by c:=c'*c;

    OPT
          (optional argument; default = `"B"'): String option passed to
          `balance' prior to ordered Schur decomposition

     *Output*
    X
          solution of the ARE

     *Method* Laub's Schur method (IEEE Transactions on Automatic
     Control, 1979) is applied to the appropriate Hamiltonian matrix
     See also: balance, dare

 -- Function File: X = dare (A, B, Q, R, OPT)
     Return the solution, X of the discrete-time algebraic Riccati
     equation
          a' x a - x + a' x b (r + b' x b)^(-1) b' x a + q = 0
     *Inputs*
    A
          N by N matrix;

    B
          N by M matrix;

    Q
          N by N matrix, symmetric positive semidefinite, or a P by N
          matrix, In the latter case q:=q'*q is used;

    R
          M by M, symmetric positive definite (invertible);

    OPT
          (optional argument; default = `"B"'): String option passed to
          `balance' prior to ordered QZ decomposition

     *Output*
    X
          solution of DARE

     *Method* Generalized eigenvalue approach (Van Dooren; SIAM J Sci.
     Stat. Comput., Vol 2) applied  to the appropriate symplectic pencil

     See also: Ran and Rodman, `Stable Hermitian Solutions of Discrete
     Algebraic Riccati Equations', Mathematics of Control, Signals and
     Systems, Vol 5, no 2 (1992), pp 165-194 See also: balance, are

 -- Function File: [TVALS, PLIST] = dre (SYS, Q, R, QF, T0, TF, PTOL,
          MAXITS)
     Solve the differential Riccati equation
          -d P/dt = A'P + P A - P B inv(R) B' P + Q
          P(tf) = Qf
     for the LTI system sys.  Solution of standard LTI state feedback
     optimization
          min int(t0, tf) ( x' Q x + u' R u ) dt + x(tf)' Qf x(tf)
     optimal input is
          u = - inv(R) B' P(t) x
     *Inputs*
    SYS
          continuous time system data structure

    Q
          state integral penalty

    R
          input integral penalty

    QF
          state terminal penalty

    T0
    TF
          limits on the integral

    PTOL
          tolerance (used to select time samples; see below); default =
          0.1

    MAXITS
          number of refinement iterations (default=10)
     *Outputs*
    TVALS
          time values at which P(T) is computed

    PLIST
          list values of P(T); PLIST { I } is P(TVALS(I))
     TVALS is selected so that:
          || Plist{i} - Plist{i-1} || < Ptol
     for every I between 2 and length(TVALS)

 -- Function File:  dgram (A, B)
     Return controllability gramian of discrete time system
          x(k+1) = a x(k) + b u(k)

     *Inputs*
    A
          N by N matrix

    B
          N by M matrix

     *Output*
    M
          N by N matrix, satisfies
               a m a' - m + b*b' = 0

 -- Function File:  dlyap (A, B)
     Solve the discrete-time Lyapunov equation

     *Inputs*
    A
          N by N matrix;

    B
          Matrix: N by N, N by M, or P by N

     *Output*
    X
          matrix satisfying appropriate discrete time Lyapunov equation

     Options:
        * B is square: solve `a x a' - x + b = 0'

        * B is not square: X satisfies either
               a x a' - x + b b' = 0
          or
               a' x a - x + b' b = 0,
          whichever is appropriate

     *Method* Uses Schur decomposition method as in Kitagawa, `An
     Algorithm for Solving the Matrix Equation X = F X F' + S',
     International Journal of Control, Volume 25, Number 5, pages
     745-753 (1977)

     Column-by-column solution method as suggested in Hammarling,
     `Numerical Solution of the Stable, Non-Negative Definite Lyapunov
     Equation', IMA Journal of Numerical Analysis, Volume 2, pages
     303-323 (1982)

 -- Function File: W = gram (SYS, MODE)
 -- Function File: WC = gram (A, B)
     `gram (SYS, 'c')' returns the controllability gramian of the
     (continuous- or discrete-time) system SYS `gram (SYS, 'o')'
     returns the observability gramian of the (continuous- or
     discrete-time) system SYS `gram (A, B)' returns the
     controllability gramian WC of the continuous-time system dx/dt = a
     x + b u; i.e., WC satisfies a Wc + m Wc' + b b' = 0

 -- Function File:  lyap (A, B, C)
 -- Function File:  lyap (A, B)
     Solve the Lyapunov (or Sylvester) equation via the Bartels-Stewart
     algorithm (Communications of the ACM, 1972)

     If A, B, and C are specified, then `lyap' returns the solution of
     the  Sylvester equation
          a x + x b + c = 0
     If only `(a, b)' are specified, then `lyap' returns the solution
     of the Lyapunov equation
          a' x + x a + b = 0
     If B is not square, then `lyap' returns the solution of either
          a' x + x a + b' b = 0
     or
          a x + x a' + b b' = 0
     whichever is appropriate

     Solves by using the Bartels-Stewart algorithm (1972)

 -- Function File:  qzval (A, B)
     Compute generalized eigenvalues of the matrix pencil
          (A - lambda B)

     A and B must be real matrices

     `qzval' is obsolete; use `qz' instead

 -- Function File: Y = zgfmul (A, B, C, D, X)
     Compute product of ZGEP incidence matrix F with vector X Used by
     `zgepbal' (in `zgscal') as part of generalized conjugate gradient
     iteration

 -- Function File:  zgfslv (N, M, P, B)
     Solve system of equations for dense zgep problem

 -- Function File: ZZ = zginit (A, B, C, D)
     Construct right hand side vector ZZ for the zero-computation
     generalized eigenvalue problem balancing procedure.  Called by
     `zgepbal'

 -- Function File:  zgreduce (SYS, MEPS)
     Implementation of procedure REDUCE in (Emami-Naeini and Van Dooren,
     Automatica, # 1982)

 -- Function File: [NONZ, ZER] = zgrownorm (MAT, MEPS)
     Return NONZ = number of rows of MAT whose two norm exceeds MEPS,
     and ZER = number of rows of mat whose two norm is less than MEPS

 -- Function File: X = zgscal (F, Z, N, M, P)
     Generalized conjugate gradient iteration to solve zero-computation
     generalized eigenvalue problem balancing equation fx=z; called by
     `zgepbal'

 -- Function File: [a, b] = zgsgiv (C, S, A, B)
     Apply givens rotation c,s to row vectors A, B No longer used in
     zero-balancing (__zgpbal__); kept for backward compatibility

 -- Function File: X = zgshsr (Y)
     Apply householder vector based on e^(m) to column vector Y Called
     by `zgfslv'

   *References*
*ZGEP*
     Hodel, `Computation of Zeros with Balancing', 1992, Linear Algebra
     and its Applications

**Generalized CG**
     Golub and Van Loan, `Matrix Computations, 2nd ed' 1989.

7 System Analysis-Properties
****************************

 -- Function File:  analdemo ()
     Octave Controls toolbox demo: State Space analysis demo

 -- Function File: [N, M, P] = abcddim (A, B, C, D)
     Check for compatibility of the dimensions of the matrices defining
     the linear system [A, B, C, D] corresponding to

          dx/dt = a x + b u
          y = c x + d u

     or a similar discrete-time system

     If the matrices are compatibly dimensioned, then `abcddim' returns

    N
          The number of system states

    M
          The number of system inputs

    P
          The number of system outputs

     Otherwise `abcddim' returns N = M = P = -1

     Note: n = 0 (pure gain block) is returned without warning See
     also: is_abcd

 -- Function File:  ctrb (SYS, B)
 -- Function File:  ctrb (A, B)
     Build controllability matrix:
          2       n-1
          Qs = [ B AB A B ... A   B ]

     of a system data structure or the pair (A, B)

     `ctrb' forms the controllability matrix The numerical properties
     of `is_controllable' are much better for controllability tests

 -- Function File:  h2norm (SYS)
     Computes the H-2 norm of a system data structure (continuous time
     only)

     Reference: Doyle, Glover, Khargonekar, Francis, `State-Space
     Solutions to Standard' `H-2 and H-infinity' `Control Problems',
     IEEE TAC August 1989

 -- Function File: [G, GMIN, GMAX] = hinfnorm (SYS, TOL, GMIN, GMAX,
          PTOL)
     Computes the H-infinity norm of a system data structure

     *Inputs*
    SYS
          system data structure

    TOL
          H-infinity norm search tolerance (default: 0.001)

    GMIN
          minimum value for norm search (default: 1e-9)

    GMAX
          maximum value for norm search (default: 1e+9)

    PTOL
          pole tolerance:
             * if sys is continuous, poles with  |real(pole))| <
               ptol*||H||  (H is appropriate Hamiltonian) are
               considered to be on the imaginary axis

             * if sys is discrete, poles with |abs(pole)-1| <
               ptol*||[s1,s2]|| (appropriate symplectic pencil) are
               considered to be on the unit circle

             * Default value: 1e-9

     *Outputs*
    G
          Computed gain, within TOL of actual gain.  G is returned as
          Inf if the system is unstable

    GMIN
    GMAX
          Actual system gain lies in the interval [GMIN, GMAX]

     References: Doyle, Glover, Khargonekar, Francis, `State-space
     solutions to standard' `H-2 and H-infinity' `control problems',
     IEEE TAC August 1989; Iglesias and Glover, `State-Space approach
     to discrete-time' `H-infinity' `control', Int. J. Control, vol 54,
     no. 5, 1991; Zhou, Doyle, Glover, `Robust and Optimal Control',
     Prentice-Hall, 1996

 -- Function File:  obsv (SYS, C)
 -- Function File:  obsv (A, C)
     Build observability matrix:
          | C        |
          | CA       |
          Qb = | CA^2     |
          | ...      |
          | CA^(n-1) |
     of a system data structure or the pair (A, C)

     The numerical properties of `is_observable' are much better for
     observability tests

 -- Function File: [ZER, POL] = pzmap (SYS)
     Plots the zeros and poles of a system in the complex plane

     *Input*
    SYS
          System data structure

     *Outputs*
    POL

    ZER
          if omitted, the poles and zeros are plotted on the screen
          otherwise, POL and ZER are returned as the system poles and
          zeros (see `sys2zp' for a preferable function call)

 -- Function File: RETVAL = is_abcd (A, B, C, D)
     Returns RETVAL = 1 if the dimensions of A, B, C, D are compatible,
     otherwise RETVAL = 0 with an appropriate diagnostic message
     printed to the screen.  The matrices B, C, or D may be omitted See
     also: abcddim

 -- Function File: [RETVAL, U] = is_controllable (SYS, TOL)
 -- Function File: [RETVAL, U] = is_controllable (A, B, TOL)
     Logical check for system controllability

     *Inputs*
    SYS
          system data structure

    A
    B
          N by N, N by M matrices, respectively

    TOL
          optional roundoff parameter.  Default value: `10*eps'

     *Outputs*
    RETVAL
          Logical flag; returns true (1) if the system SYS or the pair
          (A, B) is controllable, whichever was passed as input
          arguments

    U
          U is an orthogonal basis of the controllable subspace

     *Method* Controllability is determined by applying Arnoldi
     iteration with complete re-orthogonalization to obtain an
     orthogonal basis of the Krylov subspace
          span ([b,a*b,...,a^{n-1}*b])
     The Arnoldi iteration is executed with `krylov' if the system has
     a single input; otherwise a block Arnoldi iteration is performed
     with `krylovb' See also: size, rows, columns, length, ismatrix,
     isscalar, isvector, is_observable, is_stabilizable, is_detectable,
     krylov, krylovb

 -- Function File: RETVAL = is_detectable (A, C, TOL, DFLG)
 -- Function File: RETVAL = is_detectable (SYS, TOL)
     Test for detectability (observability of unstable modes) of (A, C)

     Returns 1 if the system A or the pair (A, C) is detectable, 0 if
     not, and -1 if the system has unobservable modes at the imaginary
     axis (unit circle for discrete-time systems)

     *See* `is_stabilizable' for detailed description of arguments and
     computational method See also: is_stabilizable, size, rows,
     columns, length, ismatrix, isscalar, isvector

 -- Function File: [RETVAL, DGKF_STRUCT ] = is_dgkf (ASYS, NU, NY, TOL )
     Determine whether a continuous time state space system meets
     assumptions of DGKF algorithm Partitions system into:
          [dx/dt]   [A  | Bw  Bu  ][w]
          [ z   ] = [Cz | Dzw Dzu ][u]
          [ y   ]   [Cy | Dyw Dyu ]
     or similar discrete-time system If necessary, orthogonal
     transformations QW, QZ and nonsingular transformations RU, RY are
     applied to respective vectors W, Z, U, Y in order to satisfy DGKF
     assumptions Loop shifting is used if DYU block is nonzero

     *Inputs*
    ASYS
          system data structure

    NU
          number of controlled inputs

    NY
          number of measured outputs

    TOL
          threshold for 0; default: 200*`eps'
     *Outputs*
    RETVAL
          true(1) if system passes check, false(0) otherwise

    DGKF_STRUCT
          data structure of `is_dgkf' results.  Entries:
         NW
         NZ
               dimensions of W, Z

         A
               system A matrix

         BW
               (N x NW) QW-transformed disturbance input matrix

         BU
               (N x NU) RU-transformed controlled input matrix;

               B = [Bw Bu]

         CZ
               (NZ x N) Qz-transformed error output matrix

         CY
               (NY x N) RY-transformed measured output matrix

               C = [Cz; Cy]

         DZU

         DYW
               off-diagonal blocks of transformed system D matrix that
               enter Z, Y from U, W respectively

         RU
               controlled input transformation matrix

         RY
               observed output transformation matrix

         DYU_NZ
               nonzero if the DYU block is nonzero

         DYU
               untransformed DYU block

         DFLG
               nonzero if the system is discrete-time
     `is_dgkf' exits with an error if the system is mixed
     discrete/continuous

     *References*
    *[1]*
          Doyle, Glover, Khargonekar, Francis, `State Space Solutions
          to Standard' `H-2 and H-infinity' `Control Problems', IEEE
          TAC August 1989

    *[2]*
          Maciejowksi, J.M., `Multivariable Feedback Design',
          Addison-Wesley, 1989

 -- Function File: DIGITAL = is_digital (SYS, EFLG)
     Return nonzero if system is digital

     *Inputs*
    SYS
          System data structure

    EFLG
          When equal to 0 (default value), exits with an error if the
          system is mixed (continuous and discrete components); when
          equal to 1, print a warning if the system is mixed
          (continuous and discrete); when equal to 2, operate silently

     *Output*
    DIGITAL
          When equal to 0, the system is purely continuous; when equal
          to 1, the system is purely discrete; when equal to -1, the
          system is mixed continuous and discrete
     Exits with an error if SYS is a mixed (continuous and discrete)
     system

 -- Function File: [RETVAL, U] = is_observable (A, C, TOL)
 -- Function File: [RETVAL, U] = is_observable (SYS, TOL)
     Logical check for system observability

     Default: tol = `tol = 10*norm(a,'fro')*eps'

     Returns 1 if the system SYS or the pair (A, C) is observable, 0 if
     not

     See `is_controllable' for detailed description of arguments and
     default values See also: size, rows, columns, length, ismatrix,
     isscalar, isvector

 -- Function File:  is_sample (TS)
     Return true if TS is a valid sampling time (real, scalar, > 0)

 -- Function File:  is_siso (SYS)
     Returns nonzero if the system data structure SYS is single-input,
     single-output

 -- Function File: RETVAL = is_stabilizable (SYS, TOL)
 -- Function File: RETVAL = is_stabilizable (A, B, TOL, DFLG)
     Logical check for system stabilizability (i.e., all unstable modes
     are controllable) Returns 1 if the system is stabilizable, 0 if
     the system is not stabilizable, -1 if the system has non
     stabilizable modes at the imaginary axis (unit circle for
     discrete-time systems

     Test for stabilizability is performed via Hautus Lemma. If DFLG!=0
     assume that discrete-time matrices (a,b) are supplied See also:
     size, rows, columns, length, ismatrix, isscalar, isvector,
     is_observable, is_stabilizable, is_detectable

 -- Function File:  is_signal_list (MYLIST)
     Return true if MYLIST is a list of individual strings

 -- Function File:  is_stable (A, TOL, DFLG)
 -- Function File:  is_stable (SYS, TOL)
     Returns 1 if the matrix A or the system SYS is stable, or 0 if not

     *Inputs*
    TOL
          is a roundoff parameter, set to 200*`eps' if omitted

    DFLG
          Digital system flag (not required for system data structure):
         `DFLG != 0'
               stable if eig(a) is in the unit circle

         `DFLG == 0'
               stable if eig(a) is in the open LHP (default)
     See also: size, rows, columns, length, ismatrix, isscalar,
     isvector, is_observable, is_stabilizable, is_detectable, krylov,
     krylovb

8 System Analysis-Time Domain
*****************************

 -- Function File:  c2d (SYS, OPT, T)
 -- Function File:  c2d (SYS, T)
     Converts the system data structure describing:

          x = Ac x + Bc u
     into a discrete time equivalent model:
          x[n+1] = Ad x[n] + Bd u[n]
     via the matrix exponential or bilinear transform

     *Inputs*
    SYS
          system data structure (may have both continuous time and
          discrete time subsystems)

    OPT
          string argument; conversion option (optional argument; may be
          omitted as shown above)
         `"ex"'
               use the matrix exponential (default)

         `"bi"'
               use the bilinear transformation
                    2(z-1)
                    s = -----
                    T(z+1)
               FIXME: This option exits with an error if SYS is not
               purely continuous. (The `ex' option can handle mixed
               systems.)

         `"matched"'
               Use the matched pole/zero equivalent transformation
               (currently only works for purely continuous SISO systems)

    T
          sampling time; required if SYS is purely continuous

          *Note* that if the second argument is not a string, `c2d()'
          assumes that the second argument is T and performs
          appropriate argument checks

     *Output*
    DSYS
          Discrete time equivalent via zero-order hold, sample each T
          sec

     This function adds the suffix  `_d' to the names of the new
     discrete states

 -- Function File:  d2c (SYS, TOL)
 -- Function File:  d2c (SYS, OPT)
     Convert a discrete (sub)system into a purely continuous one The
     sampling time used is `sysgettsam(SYS)'

     *Inputs*
    SYS
          system data structure with discrete components

    TOL
          Scalar value Tolerance for convergence of default `"log"'
          option (see below)

    OPT
          conversion option.  Choose from:
         `"log"'
               (default) Conversion is performed via a matrix logarithm
               Due to some problems with this computation, it is
               followed by a steepest descent algorithm to identify
               continuous time A, B, to get a better fit to the
               original data

               If called as `d2c (SYS, TOL)', with TOL positive scalar,
               the `"log"' option is used.  The default value for TOL
               is `1e-8'

         `"bi"'
               Conversion is performed via bilinear transform z = (1 +
               s T / 2)/(1 - s T / 2) where T is the system sampling
               time (see `sysgettsam')

               FIXME: bilinear option exits with an error if SYS is not
               purely discrete
     *Output*
    CSYS
          continuous time system (same dimensions and signal names as
          in SYS)

 -- Function File: [DSYS, FIDX] = dmr2d (SYS, IDX, SPREFIX, TS2, CUFLG)
     convert a multirate digital system to a single rate digital system
     states specified by IDX, SPREFIX are sampled at TS2, all others
     are assumed sampled at TS1 = `sysgettsam (SYS)'

     *Inputs*
    SYS
          discrete time system; `dmr2d' exits with an error if SYS is
          not discrete

    IDX
          indices or names of states with sampling time
          `sysgettsam(SYS)' (may be empty); see `cellidx'

    SPREFIX
          list of string prefixes of states with sampling time
          `sysgettsam(SYS)' (may be empty)

    TS2
          sampling time of states not specified by IDX, SPREFIX must be
          an integer multiple of `sysgettsam(SYS)'

    CUFLG
          "constant u flag" if CUFLG is nonzero then the system inputs
          are assumed to be constant over the revised sampling interval
          TS2 Otherwise, since the inputs can change during the interval
          T in [k ts2, (k+1) ts2], an additional set of inputs is
          included in the revised B matrix so that these intersample
          inputs may be included in the single-rate system default
          CUFLG = 1

     *Outputs*
    DSYS
          equivalent discrete time system with sampling time TS2

          The sampling time of sys is updated to TS2

          if CUFLG=0 then a set of additional inputs is added to the
          system with suffixes _d1, ..., _dn to indicate their delay
          from the starting time k TS2, i.e u = [u_1; u_1_d1; ...,
          u_1_dn] where u_1_dk is the input k*ts1 units of time after
          u_1 is sampled. (TS1 is the original sampling time of the
          discrete time system and TS2 = (n+1)*ts1)

    FIDX
          indices of "formerly fast" states specified by IDX and
          SPREFIX; these states are updated to the new (slower)
          sampling interval TS2

     *WARNING* Not thoroughly tested yet; especially when CUFLG == 0

 -- Function File:  damp (P, TSAM)
     Displays eigenvalues, natural frequencies and damping ratios of
     the eigenvalues of a matrix P or the A matrix of a system P,
     respectively If P is a system, TSAM must not be specified If P is
     a matrix and TSAM is specified, eigenvalues of P are assumed to be
     in Z-domain See also: eig

 -- Function File:  dcgain (SYS, TOL)
     Returns dc-gain matrix. If dc-gain is infinite an empty matrix is
     returned The argument TOL is an optional tolerance for the
     condition number of the A Matrix in SYS (default TOL = 1.0e-10)

 -- Function File: [Y, T] = impulse (SYS, INP, TSTOP, N)
     Impulse response for a linear system The system can be discrete or
     multivariable (or both) If no output arguments are specified,
     `impulse' produces a plot or the impulse response data for system
     SYS

     *Inputs*
    SYS
          System data structure

    INP
          Index of input being excited

    TSTOP
          The argument TSTOP (scalar value) denotes the time when the
          simulation should end

    N
          the number of data values

          Both parameters TSTOP and N can be omitted and will be
          computed from the eigenvalues of the A Matrix
     *Outputs*
    Y
          Values of the impulse response

    T
          Times of the impulse response
     See also: step

 -- Function File: [Y, T] = step (SYS, INP, TSTOP, N)
     Step response for a linear system The system can be discrete or
     multivariable (or both) If no output arguments are specified,
     `step' produces a plot or the step response data for system SYS

     *Inputs*
    SYS
          System data structure

    INP
          Index of input being excited

    TSTOP
          The argument TSTOP (scalar value) denotes the time when the
          simulation should end

    N
          the number of data values

          Both parameters TSTOP and N can be omitted and will be
          computed from the eigenvalues of the A Matrix
     *Outputs*
    Y
          Values of the step response

    T
          Times of the step response

     When invoked with the output parameter Y the plot is not displayed
     See also: impulse

9 System Analysis-Frequency Domain
**********************************

*Demonstration/tutorial script*

 -- Function File:  frdemo ()
     Octave Control Toolbox demo: Frequency Response demo

 -- Function File: [MAG, PHASE, W] = bode (SYS, W, OUT_IDX, IN_IDX)
     If no output arguments are given: produce Bode plots of a system;
     otherwise, compute the frequency response of a system data
     structure

     *Inputs*
    SYS
          a system data structure (must be either purely continuous or
          discrete; see is_digital)

    W
          frequency values for evaluation

          if SYS is continuous, then bode evaluates G(jw) where G(s) is
          the system transfer function

          if SYS is discrete, then bode evaluates G(`exp'(jwT)), where
             * T is the system sampling time

             * G(z) is the system transfer function

          *Default* the default frequency range is selected as follows:
          (These steps are *not* performed if W is specified)
            1. via routine __bodquist__, isolate all poles and zeros
               away from W=0 (JW=0 or `exp'(jwT)=1) and select the
               frequency range based on the breakpoint locations of the
               frequencies

            2. if SYS is discrete time, the frequency range is limited
               to jwT in [0,2 pi /T]

            3. A "smoothing" routine is used to ensure that the plot
               phase does not change excessively from point to point
               and that singular points (e.g., crossovers from +/- 180)
               are accurately shown

    OUT_IDX
    IN_IDX
          The names or indices of outputs and inputs to be used in the
          frequency response.  See `sysprune'

          *Example*
               bode(sys,[],"y_3", {"u_1","u_4"});
     *Outputs*
    MAG
    PHASE
          the magnitude and phase of the frequency response G(jw) or
          G(`exp'(jwT)) at the selected frequency values

    W
          the vector of frequency values used

       1. If no output arguments are given, e.g.,
               bode(sys);
          bode plots the results to the screen.  Descriptive labels are
          automatically placed

          Failure to include a concluding semicolon will yield some
          garbage being printed to the screen (`ans = []')

       2. If the requested plot is for an MIMO system, mag is set to
          ||G(jw)|| or ||G(`exp'(jwT))|| and phase information is not
          computed

 -- Function File: [WMIN, WMAX] = bode_bounds (ZER, POL, DFLG, TSAM)
     Get default range of frequencies based on cutoff frequencies of
     system poles and zeros Frequency range is the interval [10^WMIN,
     10^WMAX]

     Used internally in `__freqresp__' (`bode', `nyquist')

 -- Function File:  freqchkw (W)
     Used by `__freqresp__' to check that input frequency vector W is
     valid Returns boolean value

 -- Function File: OUT = ltifr (A, B, W)
 -- Function File: OUT = ltifr (SYS, W)
     Linear time invariant frequency response of single-input systems

     *Inputs*
    A
    B
          coefficient matrices of dx/dt = A x + B u

    SYS
          system data structure

    W
          vector of frequencies
     *Output*
    OUT
          frequency response, that is:

          -1
          G(s) = (jw I-A) B
     for complex frequencies s = jw

 -- Function File: [REALP, IMAGP, W] = nyquist (SYS, W, OUT_IDX,
          IN_IDX, ATOL)
 -- Function File:  nyquist (SYS, W, OUT_IDX, IN_IDX, ATOL)
     Produce Nyquist plots of a system; if no output arguments are
     given, Nyquist plot is printed to the screen

     Compute the frequency response of a system

     *Inputs* (pass as empty to get default values)
    SYS
          system data structure (must be either purely continuous or
          discrete; see `is_digital')

    W
          frequency values for evaluation If sys is continuous, then
          bode evaluates G(JW); if sys is discrete, then bode evaluates
          G(exp(JWT)), where T is the system sampling time

    DEFAULT
          the default frequency range is selected as follows: (These
          steps are *not* performed if W is specified)
            1. via routine `__bodquist__', isolate all poles and zeros
               away from W=0 (JW=0 or exp(JWT)=1) and select the
               frequency range based on the breakpoint locations of the
               frequencies

            2. if SYS is discrete time, the frequency range is limited
               to JWT in [0,2p*pi]

            3. A "smoothing" routine is used to ensure that the plot
               phase does not change excessively from point to point
               and that singular points (e.g., crossovers from +/- 180)
               are accurately shown

    ATOL
          for interactive nyquist plots: atol is a change-in-slope
          tolerance for the of asymptotes (default = 0; 1e-2 is a good
          choice).  This allows the user to "zoom in" on portions of
          the Nyquist plot too small to be seen with large asymptotes
     *Outputs*
    REALP
    IMAGP
          the real and imaginary parts of the frequency response G(jw)
          or G(exp(jwT)) at the selected frequency values

    W
          the vector of frequency values used

     If no output arguments are given, nyquist plots the results to the
     screen If ATOL != 0 and asymptotes are detected then the user is
     asked interactively if they wish to zoom in (remove asymptotes)
     Descriptive labels are automatically placed

     Note: if the requested plot is for an MIMO system, a warning
     message is presented; the returned information is of the magnitude
     ||G(jw)|| or ||G(exp(jwT))|| only; phase information is not
     computed

 -- Function File: [MAG, PHASE, W] = nichols (SYS, W, OUTPUTS, INPUTS)
     Produce Nichols plot of a system

     *Inputs*
    SYS
          System data structure (must be either purely continuous or
          discrete; see `is_digital')

    W
          Frequency values for evaluation
             * if sys is continuous, then nichols evaluates G(jw)

             * if sys is discrete, then nichols evaluates G(exp(jwT)),
               where T=SYS. TSAM is the system sampling time

             * the default frequency range is selected as follows (These
               steps are *not* performed if W is specified):
                 1. via routine `__bodquist__', isolate all poles and
                    zeros away from W=0 (jw=0 or exp(jwT)=1) and select
                    the frequency range based on the breakpoint
                    locations of the frequencies

                 2. if sys is discrete time, the frequency range is
                    limited to jwT in [0,2p*pi]

                 3. A "smoothing" routine is used to ensure that the
                    plot phase does not change excessively from point
                    to point and that singular points (e.g., crossovers
                    from +/- 180) are accurately shown

    OUTPUTS
    INPUTS
          the names or indices of the output(s) and input(s) to be used
          in the frequency response; see `sysprune'
     *Outputs*
    MAG
    PHASE
          The magnitude and phase of the frequency response G(jw) or
          G(exp(jwT)) at the selected frequency values

    W
          The vector of frequency values used
     If no output arguments are given, `nichols' plots the results to
     the screen Descriptive labels are automatically placed. See
     `xlabel', `ylabel', and `title'

     Note: if the requested plot is for an MIMO system, MAG is set to
     ||G(jw)|| or ||G(exp(jwT))|| and phase information is not computed

 -- Function File: [ZER, GAIN] = tzero (A, B, C, D, OPT)
 -- Function File: [ZER, GAIN] = tzero (SYS, OPT)
     Compute transmission zeros of a continuous system:

          x = Ax + Bu
          y = Cx + Du
     or of a discrete one:
          x(k+1) = A x(k) + B u(k)
          y(k)   = C x(k) + D u(k)

     *Outputs*
    ZER
          transmission zeros of the system

    GAIN
          leading coefficient (pole-zero form) of SISO transfer function
          returns gain=0 if system is multivariable
     *References*
       1. Emami-Naeini and Van Dooren, Automatica, 1982

       2. Hodel, `Computation of Zeros with Balancing', 1992 Lin. Alg.
          Appl

 -- Function File: ZR = tzero2 (A, B, C, D, BAL)
     Compute the transmission zeros of A, B, C, D

     BAL = balancing option (see balance); default is `"B"'

     Needs to incorporate `mvzero' algorithm to isolate finite zeros;
     use `tzero' instead

10 Controller Design
********************

 -- Function File:  dgkfdemo ()
     Octave Controls toolbox demo: H-2/H-infinity options demos

 -- Function File:  hinfdemo ()
     H-infinity design demos for continuous SISO and MIMO systems and a
     discrete system.  The SISO system is difficult to control because
     it is non-minimum-phase and unstable. The second design example
     controls the `jet707' plant, the linearized state space model of a
     Boeing 707-321 aircraft at V=80 m/s (M = 0.26, GA0 = -3 deg,
     ALPHA0 = 4 deg, KAPPA = 50 deg) Inputs: (1) thrust and (2)
     elevator angle Outputs: (1) airspeed and (2) pitch angle. The
     discrete system is a stable and second order

    SISO plant:
               s - 2
               G(s) = --------------
               (s + 2)(s - 1)

               +----+
               -------------------->| W1 |---> v1
               z   |                    +----+
               ----|-------------+
               |             |
               |    +---+    v   y  +----+
               u *--->| G |--->O--*-->| W2 |---> v2
               |    +---+       |   +----+
               |                |
               |    +---+       |
               -----| K |<-------
               +---+

               min || T   ||
               vz   infty

          W1 und W2 are the robustness and performance weighting
          functions

    MIMO plant:
          The optimal controller minimizes the H-infinity norm of the
          augmented plant P (mixed-sensitivity problem):
               w
               1 -----------+
               |                   +----+
               +---------------------->| W1 |----> z1
               w         |   |                   +----+
               2 ------------------------+
               |   |            |
               |   v   +----+   v      +----+
               +--*-->o-->| G  |-->o--*-->| W2 |---> z2
               |          +----+      |   +----+
               |                      |
               ^                      v
               u                       y (to K)
               (from controller K)

               +    +           +    +
               | z  |           | w  |
               |  1 |           |  1 |
               | z  | = [ P ] * | w  |
               |  2 |           |  2 |
               | y  |           | u  |
               +    +           +    +

    Discrete system:
          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

    The continuous plant:
               1
               G (s) = --------------
               k      (s + 2)(s + 1)

          is discretised with a ZOH (Sampling period = TS = 1 second):

               0.199788z + 0.073498
               G(z) = --------------------------
               (z - 0.36788)(z - 0.13534)

               +----+
               -------------------->| W1 |---> v1
               z   |                    +----+
               ----|-------------+
               |             |
               |    +---+    v      +----+
               *--->| G |--->O--*-->| W2 |---> v2
               |    +---+       |   +----+
               |                |
               |    +---+       |
               -----| K |<-------
               +---+

               min || T   ||
               vz   infty
          W1 and W2 are the robustness and performance weighting
          functions

 -- Function File: [L, M, P, E] = dlqe (A, G, C, SIGW, SIGV, Z)
     Construct the linear quadratic estimator (Kalman filter) for the
     discrete time system

          x[k+1] = A x[k] + B u[k] + G w[k]
          y[k] = C x[k] + D u[k] + v[k]

     where W, V are zero-mean gaussian noise processes with respective
     intensities `SIGW = cov (W, W)' and `SIGV = cov (V, V)'

     If specified, Z is `cov (W, V)'.  Otherwise `cov (W, V) = 0'

     The observer structure is

          z[k|k] = z[k|k-1] + L (y[k] - C z[k|k-1] - D u[k])
          z[k+1|k] = A z[k|k] + B u[k]

     The following values are returned:

    L
          The observer gain, (A - ALC) is stable

    M
          The Riccati equation solution

    P
          The estimate error covariance after the measurement update

    E
          The closed loop poles of (A - ALC)

 -- Function File: [K, P, E] = dlqr (A, B, Q, R, Z)
     Construct the linear quadratic regulator for the discrete time
     system

          x[k+1] = A x[k] + B u[k]

     to minimize the cost functional

          J = Sum (x' Q x + u' R u)

     Z omitted or

          J = Sum (x' Q x + u' R u + 2 x' Z u)

     Z included

     The following values are returned:

    K
          The state feedback gain, (A - BK) is stable

    P
          The solution of algebraic Riccati equation

    E
          The closed loop poles of (A - BK)

 -- Function File: [LP, LF, P, Z] = dkalman (A, G, C, QW, RV, S)
     Construct the linear quadratic estimator (Kalman predictor) for the
     discrete time system

          x[k+1] = A x[k] + B u[k] + G w[k]
          y[k] = C x[k] + D u[k] + v[k]

     where W, V are zero-mean gaussian noise processes with respective
     intensities `QW = cov (W, W)' and `RV = cov (V, V)'

     If specified, S is `cov (W, V)'.  Otherwise `cov (W, V) = 0'

     The observer structure is

          x[k+1|k] = A x[k|k-1] + B u[k] + LP (y[k] - C x[k|k-1] - D u[k])
          x[k|k] = x[k|k-1] + LF (y[k] - C x[k|k-1] - D u[k])

     The following values are returned:

    LP
          The predictor gain, (A - LP C) is stable

    LF
          The filter gain

    P
          The Riccati solution

          P = E [(x - x[n|n-1])(x - x[n|n-1])']

    Z
          The updated error covariance matrix

          Z = E [(x - x[n|n])(x - x[n|n])']

 -- Function File: [K, GAIN, KC, KF, PC, PF] = h2syn (ASYS, NU, NY, TOL)
     Design H-2 optimal controller per procedure in Doyle, Glover,
     Khargonekar, Francis, `State-Space Solutions to Standard' `H-2 and
     H-infinity' `Control Problems', IEEE TAC August 1989

     Discrete-time control per Zhou, Doyle, and Glover, `Robust and
     optimal control', Prentice-Hall, 1996

     *Inputs*
    ASYS
          system data structure (see ss, sys2ss)
             * controller is implemented for continuous time systems

             * controller is *not* implemented for discrete time systems

    NU
          number of controlled inputs

    NY
          number of measured outputs

    TOL
          threshold for 0.  Default: 200*`eps'

     *Outputs*
    K
          system controller

    GAIN
          optimal closed loop gain

    KC
          full information control (packed)

    KF
          state estimator (packed)

    PC
          ARE solution matrix for regulator subproblem

    PF
          ARE solution matrix for filter subproblem

 -- Function File: K = hinf_ctr (DGS, F, H, Z, G)
     Called by `hinfsyn' to compute the H-infinity optimal controller

     *Inputs*
    DGS
          data structure returned by `is_dgkf'

    F
    H
          feedback and filter gain (not partitioned)

    G
          final gamma value
     *Outputs*
    K
          controller (system data structure)

     Do not attempt to use this at home; no argument checking performed

 -- Function File: [K, G, GW, XINF, YINF] = hinfsyn (ASYS, NU, NY,
          GMIN, GMAX, GTOL, PTOL, TOL)
     *Inputs* input system is passed as either
    ASYS
          system data structure (see `ss', `sys2ss')
             * controller is implemented for continuous time systems

             * controller is *not* implemented for discrete time
               systems  (see bilinear transforms in `c2d', `d2c')

    NU
          number of controlled inputs

    NY
          number of measured outputs

    GMIN
          initial lower bound on H-infinity optimal gain

    GMAX
          initial upper bound on H-infinity Optimal gain

    GTOL
          Gain threshold.  Routine quits when GMAX/GMIN < 1+tol

    PTOL
          poles with `abs(real(pole))' < ptol*||H|| (H is appropriate
          Hamiltonian) are considered to be on the imaginary axis
          Default: 1e-9

    TOL
          threshold for 0.  Default: 200*`eps'

          GMAX, MIN, TOL, and TOL must all be positive scalars
     *Outputs*
    K
          System controller

    G
          Designed gain value

    GW
          Closed loop system

    XINF
          ARE solution matrix for regulator subproblem

    YINF
          ARE solution matrix for filter subproblem

     References:
       1. Doyle, Glover, Khargonekar, Francis, `State-Space Solutions
          to Standard' `H-2 and H-infinity' `Control Problems', IEEE
          TAC August 1989

       2. Maciejowksi, J.M., `Multivariable feedback design',
          Addison-Wesley, 1989, ISBN 0-201-18243-2

       3. Keith Glover and John C. Doyle, `State-space formulae for all
          stabilizing controllers that satisfy an' `H-infinity-norm'
          `bound and relations to risk sensitivity', Systems & Control
          Letters 11, Oct. 1988, pp 167-172

 -- Function File: [RETVAL, PC, PF] = hinfsyn_chk (A, B1, B2, C1, C2,
          D12, D21, G, PTOL)
     Called by `hinfsyn' to see if gain G satisfies conditions in
     Theorem 3 of Doyle, Glover, Khargonekar, Francis, `State Space
     Solutions to Standard' `H-2 and H-infinity' `Control Problems',
     IEEE TAC August 1989

     *Warning:* do not attempt to use this at home; no argument
     checking performed

     *Inputs*

     As returned by `is_dgkf', except for:
    G
          candidate gain level

    PTOL
          as in `hinfsyn'

     *Outputs*
    RETVAL
          1 if g exceeds optimal Hinf closed loop gain, else 0

    PC
          solution of "regulator" H-infinity ARE

    PF
          solution of "filter" H-infinity ARE
     Do not attempt to use this at home; no argument checking performed

 -- Function File: [XINF, X_HA_ERR] = hinfsyn_ric (A, BB, C1, D1DOT, R,
          PTOL)
     Forms
          xx = ([bb; -c1'*d1dot]/r) * [d1dot'*c1 bb'];
          Ha = [a 0*a; -c1'*c1 - a'] - xx;
     and solves associated Riccati equation The error code X_HA_ERR
     indicates one of the following conditions:
    0
          successful

    1
          XINF has imaginary eigenvalues

    2
          HX not Hamiltonian

    3
          XINF has infinite eigenvalues (numerical overflow)

    4
          XINF not symmetric

    5
          XINF not positive definite

    6
          R is singular

 -- Function File: [K, P, E] = lqe (A, G, C, SIGW, SIGV, Z)
     Construct the linear quadratic estimator (Kalman filter) for the
     continuous time system

          dx
          -- = A x + G u
          dt

          y = C x + v

     where W and V are zero-mean gaussian noise processes with
     respective intensities

          sigw = cov (w, w)
          sigv = cov (v, v)

     The optional argument Z is the cross-covariance `cov (W, V)'.  If
     it is omitted, `cov (W, V) = 0' is assumed

     Observer structure is `dz/dt = A z + B u + k (y - C z - D u)'

     The following values are returned:

    K
          The observer gain, (A - KC) is stable

    P
          The solution of algebraic Riccati equation

    E
          The vector of closed loop poles of (A - KC)

 -- Function File: [K, Q1, P1, EE, ER] = lqg (SYS, SIGW, SIGV, Q, R,
          IN_IDX)
     Design a linear-quadratic-gaussian optimal controller for the
     system
          dx/dt = A x + B u + G w       [w]=N(0,[Sigw 0    ])
          y = C x + v               [v]  (    0   Sigv ])
     or
          x(k+1) = A x(k) + B u(k) + G w(k)   [w]=N(0,[Sigw 0    ])
          y(k) = C x(k) + v(k)              [v]  (    0   Sigv ])

     *Inputs*
    SYS
          system data structure

    SIGW
    SIGV
          intensities of independent Gaussian noise processes (as above)

    Q
    R
          state, control weighting respectively.  Control ARE is

    IN_IDX
          names or indices of controlled inputs (see `sysidx',
          `cellidx')

          default: last dim(R) inputs are assumed to be controlled
          inputs, all others are assumed to be noise inputs
     *Outputs*
    K
          system data structure format LQG optimal controller (Obtain
          A, B, C matrices with `sys2ss', `sys2tf', or `sys2zp' as
          appropriate)

    P1
          Solution of control (state feedback) algebraic Riccati
          equation

    Q1
          Solution of estimation algebraic Riccati equation

    EE
          Estimator poles

    ES
          Controller poles
     See also: h2syn, lqe, lqr

 -- Function File: [K, P, E] = lqr (A, B, Q, R, Z)
     construct the linear quadratic regulator for the continuous time
     system

          dx
          -- = A x + B u
          dt

     to minimize the cost functional

          infinity
          /
          J = |  x' Q x + u' R u
          /
          t=0

     Z omitted or

          infinity
          /
          J = |  x' Q x + u' R u + 2 x' Z u
          /
          t=0

     Z included

     The following values are returned:

    K
          The state feedback gain, (A - BK) is stable and minimizes the
          cost functional

    P
          The stabilizing solution of appropriate algebraic Riccati
          equation

    E
          The vector of the closed loop poles of (A - BK)

     *Reference* Anderson and Moore, `Optimal control: linear quadratic
     methods', Prentice-Hall, 1990, pp. 56-58

 -- Function File: [Y, X] = lsim (SYS, U, T, 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 row in U
     corresponds to a different time step.  Each column 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 a plot is not displayed;
     however, the data is returned in Y (system output) and X (system
     states)

 -- Function File: K = place (SYS, P)
 -- Function File: K = place (A, B, 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)

11 Miscellaneous Functions (Not yet properly filed/documented)
**************************************************************

 -- Function File:  axis2dlim (AXDATA)
     Determine axis limits for 2-D data (column vectors); leaves a 10%
     margin around the plots Inserts margins of +/- 0.1 if data is
     one-dimensional (or a single point)

     *Input*
    AXDATA
          N by 2 matrix of data [X, Y]

     *Output*
    AXVEC
          Vector of axis limits appropriate for call to `axis' function

 -- Function File:  moddemo (INPUTS)
     Octave Control toolbox demo: Model Manipulations demo

 -- Function File:  prompt (STR)
     Prompt user to continue

     *Input*
    STR
          Input string. Its default value is:
               \n ---- Press a key to  continue ---

 -- Function File:  rldemo (INPUTS)
     Octave Control toolbox demo: Root Locus demo

 -- Function File: [RLDATA, K] = rlocus (SYS[, INCREMENT, MIN_K, MAX_K])
     Display root locus plot of the specified SISO system
          -----   ---     --------
          --->| + |---|k|---->| SISO |----------->
          -----   ---     --------        |
          - ^                             |
          |_____________________________|

     *Inputs*
    SYS
          system data structure

    MIN_K
          Minimum value of K

    MAX_K
          Maximum value of K

    INCREMENT
          The increment used in computing gain values

     *Outputs*

     Plots the root locus to the screen
    RLDATA
          Data points plotted: in column 1 real values, in column 2 the
          imaginary values

    K
          Gains for real axis break points

 -- Function File: [YY, IDX] = sortcom (XX[, OPT])
     Sort a complex vector

     *Inputs*
    XX
          Complex vector

    OPT
          sorting option:
         `"re"'
               Real part (default);

         `"mag"'
               By magnitude;

         `"im"'
               By imaginary part
          if OPT is not chosen as `"im"', then complex conjugate pairs
          are grouped together, a - jb followed by a + jb

     *Outputs*
    YY
          Sorted values

    IDX
          Permutation vector: `yy = xx(idx)'

 -- Function File: [NUM, DEN] = ss2tf (A, B, C, D)
     Conversion from transfer function to state-space The state space
     system:

          x = Ax + Bu
          y = Cx + Du

     is converted to a transfer function:

          num(s)
          G(s)=-------
          den(s)

     used internally in system data structure format manipulations

 -- Function File: [POL, ZER, K] = ss2zp (A, B, C, D)
     Converts a state space representation to a set of poles and zeros;
     K is a gain associated with the zeros

     Used internally in system data structure format manipulations

 -- Function File:  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 (i.e. negative feedback)

 -- Function File: [A, B, C, D] = tf2ss (NUM, DEN)
     Conversion from transfer function to state-space The state space
     system:

          x = Ax + Bu
          y = Cx + Du
     is obtained from a transfer function:
          num(s)
          G(s)=-------
          den(s)

     The vector DEN must contain only one row, whereas the vector NUM
     may contain as many rows as there are outputs Y of the system. The
     state space system matrices obtained from this function will be in
     controllable canonical form as described in `Modern Control
     Theory', (Brogan, 1991)

 -- Function File: [ZER, POL, K] = tf2zp (NUM, DEN)
     Converts transfer functions to poles-and-zero representations

     Returns the zeros and poles of the SISO system defined by NUM/DEN
     K is a gain associated with the system zeros

 -- Function File: [A, B, C, D] = zp2ss (ZER, POL, K)
     Conversion from zero / pole to state space

     *Inputs*
    ZER
    POL
          Vectors of (possibly) complex poles and zeros of a transfer
          function. Complex values must come in conjugate pairs (i.e.,
          x+jy in ZER means that x-jy is also in ZER) The number of
          zeros must not exceed the number of poles

    K
          Real scalar (leading coefficient)

     *Outputs*
    A
    B
    C
    D
          The state space system, in the form:

               x = Ax + Bu
               y = Cx + Du

 -- Function File: [NUM, DEN] = zp2tf (ZER, POL, K)
     Converts zeros / poles to a transfer function

     *Inputs*
    ZER
    POL
          Vectors of (possibly complex) poles and zeros of a transfer
          function.  Complex values must appear in conjugate pairs

    K
          Real scalar (leading coefficient)