Control Theory
**************

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.

OCST demo 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
           O C T A V E    C O N T R O L   S Y S T E M S   T O O L B O X
          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.

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.

System representation demo program
==================================

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

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 sysupdate).

     sys(0)
          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)

     sys(1:3)
          boolean flags to 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).

tf format variables
-------------------

num
     numerator coefficients   (vector)

den
     denominator coefficients   (vector)

zp format variables
-------------------

zer
     system zeros   (vector)

pol
     system poles    (vector)

k
     leading coefficient   (scalar)

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)

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.

Finite impulse response system interface functions
--------------------------------------------------

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

     Inputs:

     num
          vector of coefficients [c_0 c_1 ... c_n] of the SISO FIR transfer function

          C(z) = c0 + c1*z^{-1} + c2*z^{-2} + ... + znz^{-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.

     Outputs 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.

State space system interface functions
--------------------------------------

 -- Function File :  SYS = ss2sys (A,B,C{,D, TSAM, N, NZ, STNAME, INNAME, OUTNAME, OUTLIST})
     Create system structure from state-space data.   May be continous, discrete, or mixed (sampeld-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

          default:

          tsam = 0
               n = rows(A), nz = 0

          tsam > 0
                n = 0,       nz = rows(A)

               see below for system partitioning

     stname
          list of strings of state signal names

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

     inname
          list of strings of input signal names

          default (INNAME = [] on input): u_n

     outname
          list of strings of input signal names

          default (OUTNAME = [] on input): y_n

     outlist

          list of indices of outputs y that are sampled

          default:

          tsam = 0
               outlist = []

          tsam > 0
               outlist = 1:rows(C)

          Unlike states, discrete/continous outputs may appear in any order.

          Note 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 = ss2sys(a,b,c,[],0,3,0,list("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.

     Inputs SYS system data structure (see sysstruct)

     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) 0.

     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

Transfer function system interface functions
--------------------------------------------

 -- Function File :  SYS =  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 list with a single string entry.

     Outputs 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 tf2sys for parameter descriptions.

     Example

          octave:1> sys=ss2sys([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

Zero-pole system interface functions
------------------------------------

 -- Function File :  SYS = 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)

     Outputs 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 zp2sys for parameter descriptions.

     Example

          octave:1> sys=ss2sys([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

Data structure access functions
-------------------------------

 -- Function File : RETSYS = syschnames (SYS, OPT, LIST, NAMES)
     Superseded by syssetsignals

 -- Function File :  retsys = 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 continous

 -- 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 of signal (or indices of signals if signum is a vector)

     strflg
          flag to return a string instead of a list;  Values:

          0
               (default) return a list (even if signum is a scalar)

          1
               return a string.  Exits with an error if signum is not a scalar.

     Outputs

     If SIGID is not specified

          stname, inname, outname
               signal names (lists 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, then

          sigid="in"
               SIGLIST is set to the list of input names

          sigid="out"
               SIGLIST is set to the list of output names

          sigid="st"
               SIGLIST is set to the list 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, then SIGNAME is
          a list 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=ss2sys(rand(4),rand(4,2),rand(3,4));
          octave> [Ast,Ain,Aout,Ayd] = sysgetsignals(sys) i  # get all signal names
          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> Ain = sysgetsignals(sys,"in")   # get only input signal names
          Ain =
          (
            [1] = u_1
            [2] = u_2
          )
          octave> Aout = sysgetsignals(sys,"out",2)   # get name of output 2 (in list)
          Aout =
          (
            [1] = y_2
          )
          octave> Aout = sysgetsignals(sys,"out",2,1)  # get name of output 2 (as string)
          Aout = y_2

 -- Function File :  Tsam = sysgettsam ( sys )
     return the sampling time of the system

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

     Inputs SYS: system data structure

     Outputs SYSTYPE: string indicating how the structure was initially constructed: values: "ss", "zp", or "tf"

     Note FIR initialized systems return systype="tf".

 -- Function File : RETSYS = 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", or "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.

     list
          vector of indices of outputs, yd, inputs, or states whose respective names should be changed.

          Default: replace entire list of names/entire yd vector.

     Outputs RETSYS=SYS with appropriate signal names changed (or yd values, where appropriate)

     Example

          octave:1> sys=ss2sys([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 :  SYS = 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.

Data structure internal functions
---------------------------------

 -- Function File :   syschnamesl
     used internally in syschnames

 -- Function File :  IONAME = sysdefioname (N,STR {,M})
     return default input or output names given N, STR, M. N is the final value, STR is the string prefix, and M is start value

     used internally, minimal argument checking

     Example ioname = sysdefioname(5,"u",3) returns the list:

          ioname =
          (
            [1] = u_3
            [2] = u_4
            [3] = u_5
          )

 -- Function File :  STNAME = sysdefstname (N, NZ)
     return default state names given N, NZ

     used internally, minimal argument checking

 -- Function File :  VEC =  tf2sysl (VEC)
     used internally in tf2sys. strip leading zero coefficients to get the true polynomial length

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); see sysgettype.

          "ss"
               state space form

          "tf"
               transfer function form

          "zp"
               zero-pole form

          "all"
               all of the above

 -- Function File :  Y = polyout ( C{, X})
     write formatted polynomial

             c(x) = c(1) * x^n + ... + c(n) x + c(n+1)

     to string Y or to the screen (if Y is omitted) X defaults to the string "s"

 -- Function File :   tfout (NUM, DENOM{, X})
     print formatted transfer function n(s)/d(s)  to the screen X defaults to the string "s"

 -- Function File :   zpout (ZER, POL, K{, X})
     print formatted zero-pole form to the screen. X defaults to the string "s"

 -- Function File :   outlist (LMAT{, TABCHAR, YD, ILIST })
     Prints an enumerated list of strings. internal use only; minimal argument checking performed

     Inputs

     lmat
          list of strings

     tabchar
          tab character (default: none)

     yd
          indices of strings to append with the string "(discrete)" (used by SYSOUT; minimal checking of this argument) yd = []  indicates all outputs are continuous

     ilist
          index numbers to print with names.

          default: 1:rows(lmat)

     Outputs prints the list to the screen, numbering each string in order.

Block Diagram Manipulations
===========================

See systime

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

 -- Function File :  outputs = bddemo ( inputs )
     Octave Controls toolbox demo: Block Diagram Manipulations demo

 -- Function File : [SYS] = buildssic(CLST, ULST, OLST, ILST, S1, S2, S3, S4, S5, S6, S7, S8)

     Contributed by Kai Mueller.

     Form an arbitrary complex (open or closed loop) system in state-space form from several systems. "buildssic" can easily (despite it's 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 H2 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] ==> new input 1 is old inpout 1 + output 2 - output 1, 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, specifiy the outputs of the resulting systems. Elements are output numbers of s1-s8. The numbers are alowed to be negative and may appear in any order. An empty matrix means all outputs.

     Ilst
          input list, specifiy the inputs of the resulting systems. Elements are input numbers of s1-s8. The numbers are alowed 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 optained by

          GW = buildssic([1 2; 2 -1], [2], [1 2 3], [2], G, K);

     Clst
          (1. row) connect input 1 (G) with output 2 (K). (2. row) connect input 2 (K) with neg. 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   |                    +----+
          ----|-------------+                   || GW   ||     => min.
              |             |                        vz   infty
              |    +---+    v      +----+
              *--->| G |--->O--*-->| W2 |---> v2
              |    +---+       |   +----+
              |                |
              |                v
             u                  y

     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 (auxillary) function with order 0. (e.g. One = ugain(1);)

 -- Function File :  OUTSYS  = jet707 ( )
     Creates linearized state space model of a Boeing 707-321 aircraft at v=80m/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 Ref: R. Brockhaus: Flugregelung (Flight Control), Springer, 1994

     see also: ord2

     Contributed by Kai Mueller

 -- Function File :  OUTSYS = 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 ommitted.

     Outputs 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 :  SYSP = parallel(ASYS, BSYS)
     Forms the parallel connection of two systems.

                   ____________________
                   |      ________    |
          u  ----->|----> | Asys |--->|----> y1
              |    |      --------    |
              |    |      ________    |
              |--->|----> | Bsys |--->|----> y2
                   |      --------    |
                   --------------------
                        Ksys

 -- Function File :  SYS = 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 : RETSYS = sysappend (SYS,B{, C, D, OUTNAME, INNAME, YD})
     appends new inputs and/or outputs to a system

     Inputs

     sys
          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 := [sys.b , b]
             sys.c := [sys.c  ]
                      [ c     ]
             sys.d := [sys.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 = continuous or mixed YD = zeros(1,rows(c))

       * SYS = discrete YD = ones(1,rows(c))

 -- Function File : RETSYS = 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
          list of connections indices; y(out_idx(ii)) is connected to u(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: 200EPS

     Outputs SYS: 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.

     Inputs SYS is a 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.

          returns CSYS empty if no continuous/continous path exists

 -- Function File :  [N_TOT, ST_C, ST_D, Y_C, Y_D] = syscont_disc(SYS)
     Used internally in syscont and sysdisc.

     Inputs SYS is a system data structure.

     Outputs

     n_tot
          total number of states

     st_c
          vector of continuous state indices (empty if none)

     st_d
          vector of discrete state indices (empty if none)

     y_c
          vector of continuous output indices

     y_d
          vector of discrete output indices

 -- Function File :  [DSYS, ADC, CDC] = sysdisc (SYS)

     Inputs 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 (see ss2sys)

     out_idx,in_idx
          list of connections indices; duplicates are made of Y(OUT_IDX(II)) and U(IN_IDX(II)).

     Outputs 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

     Outputs 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 :  NAMES = sysgroupn (NAMES)
     Locate and mark duplicate names.  Used internally in sysgroup (and elsewhere).

 -- 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 of Bsys 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
          list of connections indices; the new system has outputs y(out_idx(ii)) and inputs u(in_idx(ii)). May select as [] (empty matrix) to specify all outputs/inputs.

     Outputs RETSYS: resulting system

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

 -- Function File :  PV = sysreorder( VLEN, {varlist)

     Inputs VLEN=vector length, LIST= a subset of [1:vlen],

     Outputs 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 : SYS = sysscale (SYS, OUTSCALE, INSCALE{, OUTNAME, INNAME})
     scale inputs/outputs of a system.

     Inputs sys: structured system outscale, inscale: constant matrices of appropriate dimension

     Outputs SYS: 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)
     returns 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 :  OUTSYS = 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 (MIMO H_infinty example, Boeing 707-321 aircraft model)

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

     Weighting function are needed by the H2/H_infinity design procedure. These function are part of thye augmented plant P (see hinfdemo for an applicattion example).

     vl = Gain @ low frequencies

     vh = Gain @ high frequencies

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

Numerical Functions
===================

 -- Function File: 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
          NxN matrix.

     b
          NxN matrix or NxM matrix; in the latter case B is replaced by b:=b*b'.

     c
          NxN matrix or PxM 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.

     Outputs X: solution of the ARE.

     Method Laub's Schur method (IEEE Transactions on Automatic Control, 1979) is applied to the appropriate Hamiltonian matrix.

 -- Function File: dare (A, B, C, 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 + c = 0

     Inputs

     a
          N by N.

     b
          N by M.

     c
          N by N, symmetric positive semidefinite, or P by N. In the latter case c:=c'*c 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.

     Outputs 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.

 -- Function File :  M = dgram ( A, B)
     Return controllability grammian of discrete time system

            x(k+1) = a x(k) + b u(k)

     Inputs

     a
          N by N matrix

     b
          N by M matrix

     Outputs M (N by N) satisfies

           a m a' - m + b*b' = 0

 -- Function File: X =  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.

     Outputs 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 = FXF' + 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 :  M = gram (A, B)
     Return controllability grammian M of the continuous time system dx/dt = a x + b u.

     M satisfies  a m + m a' + 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.

 -- Function File :   pinv ( X{,TOL} )
     Returns the pseudoinverse of X; singular values less than tol are ignored.

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

          (A - lambda B).

     A and B must be real matrices.

     Note qzval is obsolete; use qz instead.

 -- Function File :   zgfmul
 -- Function File :   zgfslv
 -- Function File :   zginit
 -- Function File : RETSYS = zgpbal (ASYS)
 -- Function File : zgreduce
 -- Function File :  [NONZ, ZER] = zgrownorm (MAT, MEPS)
 -- Function File :  x = zgscal (F, Z, N, M, P)
 -- Function File :   zgsgiv ( )
 -- Function File :  X = zgshsr( Y)
     Used internally by tzero. Minimal argument checking performed.

     Details involving major subroutines:

     zgpbal
          Implementation of zero computation generalized eigenvalue problem balancing method.  zgpbal computes a state/input/output weighting that attempts to reduced the range of the magnitudes of the nonzero elements of [a,b,c,d] The weighting uses scalar multiplication by powers of 2, so no roundoff will occur.

          zgpbal should be followed by zgpred

     zgreduce
          Implementation of procedure REDUCE in (Emami-Naeini and Van Dooren, Automatica, 1982).

     zgrownorm
          Returns NONZ = number of rows of MAT whose two norm exceeds MEPS, ZER = number of rows of mat whose two norm is less than meps

     zgscal
          Generalized conjugate gradient iteration to solve zero-computation generalized eigenvalue problem balancing equation fx=z; called by zgepbal

     zgsgiv
          apply givens rotation c,s to column vector a,b

     zgshsr
          Apply Householder vector based on e^(m) (all ones) 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

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 : [Y, MY, NY] = abcddims (X)

     Used internally in abcddim.  If X is a zero-size matrix, both dimensions are set to 0 in Y. MY and NY are the row and column dimensions of the result.

 -- Function File : QS = ctrb(SYS {, B})
 -- Function File : QS = 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).

     Note ctrb forms the controllability matrix. The numerical properties of is_controllable are much better for controllability tests.

 -- Function File : RETVAL = h2norm(SYS)
     Computes the H2 norm of a system data structure (continuous time only)

     Reference: Doyle, Glover, Khargonekar, Francis, ``State Space Solutions to Standard H2 and Hinf 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: 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 H2 and Hinf control problems", IEEE TAC August 1989 Iglesias and Glover, "State-Space approach to discrete-time Hinf control," Int. J. Control, vol 54, #5, 1991 Zhou, Doyle, Glover, "Robust and Optimal Control," Prentice-Hall, 1996 $Revision: 1.9 $

 -- Function File :  QB = obsv (SYS{, C})
     Build observability matrix

               | C        |
               | CA       |
          Qb = | CA^2     |
               | ...      |
               | CA^(n-1) |

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

     Note: obsv() forms the observability matrix.

     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. Inputs SYS system data structure

     Outputs if omitted, the poles and zeros are plotted on the screen. otherwise, pol, zer are returned as the system poles and zeros. (see sys2zp for a preferable function call)

Function File outputs =  synKnames (inputs) Return controller signal names based in plant signal names and dimensions

 -- 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.

 -- 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 paramter.  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 is_observable, is_stabilizable, is_detectable, krylov, krylovb

 -- Function File :  [RETVAL, U] = is_detectable (A, C{, TOL})
 -- Function File :  [RETVAL, U] = is_detectable (SYS{, TOL})
     Test for detactability (observability of unstable modes) of (A,C).

     Returns 1 if the system A or the pair (A,C)is detectable, 0 if not.

     See is_stabilizable for detailed description of arguments and computational method.

 -- 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
          threshhold for 0.  Default: 200EPS

     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;

               Note B = [Bw Bu]

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

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

               Note C = [Cz; Cy]

          Dzu, Dyw
               off-diagonal blocks of transformed 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 H2 and Hinf Control Problems," IEEE TAC August 1989

     [2]
          Maciejowksi, J.M.: "Multivariable feedback design,"

 -- Function File :  RETVAL = is_digital ( SYS)
     Return nonzero if system is digital; Exits with an error of 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. 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.

 -- Function File :  RETVAL = is_sample (TS)
     return true if TS is a legal sampling time (real,scalar, > 0)

 -- Function File :  RETVAL = is_siso (SYS)
     return nonzero if the system data structure SYS is single-input, single-output.

 -- Function File : [RETVAL, U] = is_stabilizable (SYS{, TOL})
 -- Function File : [RETVAL, U] = is_stabilizable (A{, B ,TOL})
     Logical check for system stabilizability (i.e., all unstable modes are controllable).

     See is_controllable for description of inputs, outputs. Test for stabilizability is performed via an ordered Schur decomposition that reveals the unstable subspace of the system A matrix.

 -- Function File : FLG = is_signal_list (MYLIST)
     returns true if mylist is a list of individual strings (legal for input to SYSSETSIGNALS).

 -- Function File :  RETVAL = is_stable (A{,TOL,DFLG})
 -- Function File :  RETVAL = is_stable (SYS{,TOL})
     Returns retval = 1 if the matrix A or the system SYS is stable, or 0 if not.

     Inputs

     tol
          is a roundoff paramter, set to 200*EPS if omitted.

     dflg
          Digital system flag (not required for system data structure):

          DFLG != 0
               stable if eig(a) in unit circle

          DFLG == 0
               stable if eig(a) in open LHP (default)

System Analysis-Time Domain
===========================

 -- Function File :  DSYS = c2d (SYS{, OPT, T})
 -- Function File :  DSYS = c2d (SYS{, T})

     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.)

     T
          sampling time; required if sys is purely continuous.

          Note If the 2nd argument is not a string, c2d assumes that the 2nd argument is T and performs appropriate argument checks.

     Outputs DSYS discrete time equivalent via zero-order hold, sample each T sec.

     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

     Note This function adds the suffix  _d to the names of the new discrete states.

 -- Function File : CSYS = d2c (SYS{,TOL})
 -- Function File : CSYS = d2c (SYS, OPT)
     Convert discrete (sub)system to a purely continuous system.  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), 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

     Outputs 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
          list of states with sampling time sysgettsam(SYS) (may be empty)

     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 discrete time sys 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 : GM = 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 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 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, T: impulse response

 -- Function File : [Y, T] = step (SYS{, INP,TSTOP, N})
     Step response of a linear system; calling protocol is identical to impulse.

 -- Function File :   stepimp ( )
     Used internally in impulse, step.

System Analysis-Frequency Domain
================================

Demonstration/tutorial script
 -- Function File :   frdemo ( )

 -- 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=sysgettsam(SYS) (the system sampling time) and

            * 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 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

     Notes

      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 :  [F, W] = bodquist (SYS, W, OUT_IDX, IN_IDX)
     used internally by bode, nyquist; compute system frequency response.

     Inputs

     sys
          input system structure

     w
          range of frequencies; empty if user wants default

     out_idx
          list of outputs; empty if user wants all

     in_idx
          list of inputs; empty if user wants all

     rname
          name of routine that called bodquist ("bode" or "nyquist")

     Outputs

     w
          list of frequencies

     f
          frequency response of sys; f(ii) = f(omega(ii))

     Note bodquist could easily be incorporated into a Nichols plot function; this is in a "to do" list.

 -- Function File :  RETVAL = freqchkw ( W )
     Used by freqresp to check that input frequency vector W is legal. Returns boolean value.

 -- Function File :  OUT = freqresp (SYS,USEW{,W});
     Frequency response function - used internally by bode, nyquist. minimal argument checking; "do not attempt to do this at home"

     Inputs

     sys
          system data structure

     USEW
          returned by freqchkw

     optional
          must be present if USEW is true (nonzero)

     Outputs

     OUT
          vector of finite G(j*w) entries (or ||G(j*w)|| for MIMO)

     w
          vector of corresponding frequencies

 -- 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

     Outputs OUT

                                     -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.

     Arguments are identical to bode with exceptions as noted below:

     Inputs (pass as empty to get default values)

     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: tzero (A, B, C, D{, OPT})
 -- Function File: tzero (SYS{,OPT})
     Compute transmission zeros of a continuous

          .
          x = Ax + Bu
          y = Cx + Du

     or discrete

          x(k+1) = A x(k) + B u(k)
          y(k)   = C x(k) + D u(k)

     system. 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.

Controller Design
=================

 -- Function File :   dgkfdemo ( )
     Octave Controls toolbox demo: H2/Hinfinity options demos

 -- Function File :   hinfdemo ( )
     Non-trivial H_infinity design demo.

     H_infinity optimal controller for the jet707 plant; Linearized state space model of a Boeing 707-321 aircraft at v=80m/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 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 (from                 y (to K)
                     controller
                     K)

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

 -- 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] + w[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+1] = A z[k] + B u[k] + k (y[k] - C z[k] - D 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).

     References

      1. Anderson and Moore, Optimal Control: Linear Quadratic Methods, Prentice-Hall, 1990, pp. 56-58

      2. Kuo, Digital Control Systems, Harcourt Brace Jovanovich, 1992, section 11-5-2.

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

     Inputs input system is passed as either

     Asys
          system data structure (see ss2sys, 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
          threshhold 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_inf optimal controller.

     Inputs

     dgs
          data structure returned by is_dgkf

     F, H
          feedback and filter gain (not partitioned)

     g
          final gamma value

     Outputs controller K (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 ss2sys, 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 threshhold.  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
          threshhold for 0.  Default: 200*eps

          GMAX, MIN, TOL, and TOL must all be postive 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

      1. Doyle, Glover, Khargonekar, Francis, "State Space Solutions to Standard H2 and Hinf 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 and h-infinity-norm bound and relations to risk sensitivity," Systems & Control Letters 11, Oct. 1988, pp 167-172.

 -- Function File : [XX, ERR] = hinfsyn_c (NN, PTOL, S1{, S2})
     used internally in hinfsyn to evaluate hamiltonian/symplectic eigenvalue problems.

     WARNING Argument checking not performed.

     Inputs

     s1 (alone)
          hamiltonian matrix

     (s1,s2)  as a pair
          symplectic matrix pencil

     Outputs

     xx: positive (semi-)definite solution of DARE (set to 0 if err <=2)

     code:  error:

          0
               no error

          1
               (s1): eigenvalues on imaginary axis

               (s1,s2): gen. eigenvalues on unit circle

          2
               unequal number of stable/antistable (generalized) eigenvalues

          3
               (s1): infinite entries in solution x

               (s1,s2): infinite entires in solution x or (I + R X) singular

          4
               x is not symmetric

          5
               x has negative eigenvalues

     Solution method: Either Laub's schur method or Symplectic GEP approach; uses Van Dooren's code to re-order qz decompostion (www.netlib.com - toms/590)

     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.

 -- 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 H2 and Hinf 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-inf ARE Pf: solution of "filter" H-inf ARE

 -- 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 + b u
          dt

          y = c x + d u

     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, Q, P, 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
          indices of controlled inputs

          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)

     P
          Solution of control (state feedback) algebraic Riccati equation

     Q
          Solution of estimation algebraic Riccati equation

     Ee
          estimator poles

     Es
          controller poles

 -- 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.

     p
          The stabilizing solution of appropriate algebraic Riccati equation.

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

 -- Function File :   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 column in u corresponds to a different time step.  Each row 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 with output parameters: [y,x] = lsim(sys,u,t,[x0]) a plot is not displayed, however, the data is returned in y = system output and x = system states.

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

     Written by: Jose Daniel Munoz Frias.

Miscellaneous Functions (Not yet properly filed/documented)
===========================================================

Function File   AXVEC = axis2dlim (AXDATA) determine axis limits for 2-d data(column vectors); leaves a 10% margin around the plots. puts in margins of +/- 0.1 if data is one dimensional (or a single point)

Inputs AXDATA nx2 matrix of data [x,y]

Outputs AXVEC vector of axis limits appropriate for call to axis() function

 -- Function File :  outputs = mb ( inputs )

           $Revision: 1.9 $

 -- Function File :  outputs = moddemo ( inputs )

           Octave Controls toolbox demo: Model Manipulations demo
           Written by David Clem August 15, 1994

 -- Function File :  outputs = prompt ( inputs )

           function prompt([str])
           Prompt user to continue
           str: input string. Default value: "\n ---- Press a key to continue ---"
           Written by David Clem August 15, 1994
           Modified A. S. Hodel June 1995

 -- Function File :  outputs = rldemo ( inputs )

 -- Function File :  outputs = rlocus ( inputs )

           [rldata, k] = rlocus(sys[,increment,min_k,max_k])
           Displays root locus plot of the specified SISO system.

                 -----   ---     --------
             --->| + |---|k|---->| SISO |----------->
                 -----   ---     --------        |
                 - ^                             |
                   |_____________________________|

          inputs: sys = system data structure
                  min_k, max_k,increment: minimum, maximum values of k and
                         the increment used in computing gain values
           Outputs: plots the root locus to the screen.
             rldata: Data points plotted column 1: real values, column 2: imaginary
                     values)
             k: gains for real axis break points.

 -- Function File :  outputs = sortcom ( inputs )

           [yy,idx] = sortcom(xx[,opt]): sort a complex vector
           xx: complex vector
           opt: sorting option:
          	"re": real part (default)
          	"mag": by magnitude
          	"im": by imaginary part

            if opt != "im" then complex conjugate pairs are grouped together,
               a - jb followed by a + jb.
           yy: sorted values
           idx: permutation vector: yy = xx(idx)

 -- Function File :  outputs = ss2tf ( inputs )

           [num,den] = ss2tf(a,b,c,d)
           Conversion from tranfer 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 :  outputs = ss2zp ( inputs )

           Converts a state space representation to a set of poles and zeros.

           [pol,zer,k] = ss2zp(a,b,c,d) returns the poles and zeros of the state space
           system (a,b,c,d).  K is a gain associated with the zeros.

           used internally in system data structure format manipulations

 -- Function File :  outputs = starp ( inputs )

           [sys] = 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 (= negative feedback)

 -- Function File :  outputs = susball ( inputs )

 -- Function File :  outputs = swap ( inputs )

           [a1,b1] = swap(a,b)
           interchange a and b

 -- Function File :  outputs = swapcols ( inputs )

           function B = swapcols(A)
           permute columns of A into reverse order

 -- Function File :  outputs = swaprows ( inputs )

           function B = swaprows(A)
           permute rows of A into reverse order

 -- Function File :  outputs = tf2ss ( inputs )

           Conversion from tranfer 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)

           via the function call [a,b,c,d] = tf2ss(num,den).
           The vector 'den' must contain only one row, whereas the vector 'num'
           may contain as many rows as there are outputs of the system 'y'.
           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 :  outputs = tf2zp ( inputs )

           Converts transfer functions to poles / zeros.

           [zer,pol,k] = tf2zp(num,den) returns the zeros and poles of the SISO system
           defined by num/den.  K is a gain associated with the system zeros.

 -- Function File :   zp2ss
     Conversion from zero / pole to state space. The state space system

          .
          x = Ax + Bu
          y = Cx + Du

     is obtained from a vector of zeros and a vector of poles via the function call [a,b,c,d] = zp2ss(zer,pol,k). The vectors zer and pol may either be row or column vectors.  Each zero and pole that has an imaginary part must have a conjugate in the list. The number of zeros must not exceed the number of poles. k is zp-form leading coefficient.

 -- Function File :  [POLY, RVALS] = zp2ssg2 (RVALS)
     Used internally in zp2ss Extract 2 values from RVALS (if possible) and construct a polynomial with those roots.

 -- Function File :   zp2tf
     Converts zeros / poles to a transfer function.

     [num,den] = zp2tf(zer,pol,k) forms the transfer function num/den from the vectors of poles and zeros.