Go to the first, previous, next, last section, table of contents.
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
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
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.
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 section System Construction and Interface Functions).
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.
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 section Data structure access functions).
- 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).
- num
-
numerator coefficients (vector)
- den
-
denominator coefficients (vector)
- zer
-
system zeros (vector)
- pol
-
system poles (vector)
- k
-
leading coefficient (scalar)
- 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)
Construction and manipulations of the OCST system data structure
(See section System Data Structure) 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 section System display functions.
- Function File : sys = fir2sys ( num{, tsam, inname, outname } )
-
construct a system data structure from FIR description
Inputs:
- num
-
vector of coefficients @math{[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 section Finite impulse response system interface functions for
parameter descriptions.
- 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 (sampeled-data)
Inputs
- a, b, c, d
-
usual state space matrices.
default: d = zero matrix
- tsam
-
sampling rate. Default: @math{tsam = 0} (continuous system)
- n, nz
-
number of continuous, discrete states in the system
default:
- tsam = 0
-
@math{n =
rows
(a)}, @math{nz = 0}
- tsam > 0
-
@math{ n = 0}, @math{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
-
@math{outlist = []}
- tsam > 0
-
@math{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:
@math{ d/dt xc(t) = acc*xc(t) + acd*xd(k*tsam) + bc*u(t)}
@math{ xd((k+1)*tsam) = adc*xc(k*tsam) + add*xd(k*tsam) + bd*u(k*tsam)}
@math{ yc(t) = ccc*xc(t) + ccd*xd(k*tsam) + dc*u(t)}
@math{ 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 @math{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 section 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) 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
- 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 : 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 section Zero-pole system interface functions 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
- 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.
@math{yd(ii) = 0} if output ii is continous
- 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 : 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.
- Function File : syschnamesl
-
used internally in syschnames
item olist: index list
old_names: original list names
inames: new names
listname: name of index list
combines the two string lists old_names and inames
- 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 : 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"
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"
- 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)
@math{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.
See section System Analysis-Time Domain
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
See also: jet707 (MIMO example, Boeing 707-321 aircraft model)
- 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 = 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; @math{y(out_idx(ii))}
is connected to @math{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 : sys = sysgroup ( Asys, Bsys)
-
Combines two systems into a single system
Inputs
Asys, Bsys: system data structures
Outputs
@math{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)
-
names = sysgroupn(names)
Locate and mark duplicate names
inputs:
names: list of signal names
kind: kind of signal name (used for diagnostic message purposes only)
outputs:
returns names with unique suffixes added; diagnostic warning
message is printed to inform the user of the new signal name
used internally in sysgroup and elsewhere.
- Function File : sys = sysmult( Asys, Bsys)
-
Compute @math{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, {var{list})
-
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 @math{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 : 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)
- 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 @math{b:=b*b'}.
- c
-
nxn matrix or pxm matrix; in the latter case
c is replaced by @math{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 @math{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: 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:
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
@math{ dx/dt = a x + b u}.
m satisfies @math{ 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.
Solves by using the Bartels-Stewart algorithm (1972).
- Loadable Function: pinv (x, tol)
-
Return the pseudoinverse of x. Singular values less than
tol are ignored.
If the second argument is omitted, it is assumed that
tol = max (size (x)) * sigma_max (x) * eps,
where sigma_max (x)
is the maximal singular value of x.
- 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 : 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 : x = 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 : retsys = zgpbal(Asys)
-
used internally in
tzero
; minimal argument checking performed
implementation of zero computation generalized eigenvalue problem
balancing method (Hodel and Tiller, Allerton Conference, 1991)
Based on Ward's balancing algorithm (SIAM J. Sci Stat. Comput., 1981)
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
- Function File : retsys = zgreduce(Asys,meps)
-
Implementation of procedure REDUCE in (Emami-Naeini and Van Dooren,
Automatica, # 1982).
- Function File : [nonz, zer] = zgrownorm (mat, meps)
-
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
References:
ZGEP: Hodel, "Computation of Zeros with Balancing," 1992, submitted to LAA
Generalized CG: Golub and Van Loan, "Matrix Computations, 2nd ed" 1989
- 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
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
- 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 : 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.
- 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.
Default: tol = 10*norm(a,'fro')*eps
- 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 @math{B = [Bw Bu] }
- Cz
-
(nz x n) Qz-transformed error output matrix
- Cy
-
(ny x n) Ry-transformed measured output matrix
Note @math{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,U] = is_observable (a, c{,tol})
-
- Function File : [retval,U] = is_observable (sys{, tol})
-
Logical check for system observability.
Default: 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.
- 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).
Test for stabilizability is performed via an ordered Schur decomposition
that reveals the unstable subspace of the system A matrix.
Returns
retval
= 1 if the system, a
, is stabilizable, if the pair
(a
, b
) is stabilizable, or 0 if not.
U
= orthogonal basis of controllable subspace.
Controllable subspace is determined by applying Arnoldi iteration with
complete re-orthogonalization to obtain an orthogonal basis of the
Krylov subspace.
span ([b,a*b,...,a^ b]).
tol is a roundoff paramter, set to 200*eps if omitted.
- 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)
- 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
@math{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 @math{[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 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, t: impulse response
- Function File : [y, t] = impulse (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,
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
When invoked with the output paramter y the plot is not displayed.
- Function File : [y, t] = stepimp(sitype,sys[, inp, tstop, n])
-
Impulse or step response for a linear system.
The system can be discrete or multivariable (or both).
This m-file contains the "common code" of step and impulse.
Produces a plot or the response data for system sys.
Limited argument checking; "do not attempt to do this at home".
Used internally in
impulse
, step
. Use step
or impulse
instead.
Demonstration/tutorial script
- Function File : frdemo ( )
-
Octave Controls 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 @math{G(jw)} where
@math{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
- @math{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)
- via routine bodquist, isolate all poles and zeros away from
w=0 (jw=0 or @math{
exp
(jwT)}=1) and select the frequency
range based on the breakpoint locations of the frequencies.
- if sys is discrete time, the frequency range is limited
to @math{jwT} in
[0,2 pi /T]
- 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
@math{G(jw)} or @math{G(
exp
(jwT))} at the selected frequency values.
- w
-
the vector of frequency values used
Notes
- 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 = []
).
- If the requested plot is for an MIMO system, mag is set to
@math{||G(jw)||} or @math{||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; @math{f(ii) = f(omega(ii))}
Note bodquist could easily be incorporated into a Nichols
plot function; this is in a "to do" list.
Both bode and nyquist share the same introduction, so the common parts are
in bodquist. It contains the part that finds the number of arguments,
determines whether or not the system is SISO, and computes the frequency
response. Only the way the response is plotted is different between the
two functions.
- 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 @math{G(jw)}
if sys is discrete, then bode evaluates @math{G(exp(jwT))}, where
@math{T=sysgettsam(sys)} (the system sampling time)
- default
-
the default frequency range is selected as follows: (These
steps are NOT performed if w is specified)
- via routine bodquist, isolate all poles and zeros away from
w=0 (jw=0 or @math{exp(jwT)=1}) and select the frequency
range based on the breakpoint locations of the frequencies.
- if sys is discrete time, the frequency range is limited
to jwT in
[0,2p*pi]
- 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 indices of the output(s) and input(s) to be used in
the frequency response; see sysprune.
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
@math{G(jw)} or @math{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 : 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.
- Function File : dgkfdemo ( )
-
Octave Controls toolbox demo: H2/Hinfinity 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=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 discrete system is a
stable and second order.
- SISO plant
-
s - 2
G(s) = --------------
(s + 2)(s - 1)
+----+
-------------------->| W1 |---> v1
z | +----+
----|-------------+ || T || => min.
| | vz infty
| +---+ v y +----+
u *--->| G |--->O--*-->| W2 |---> v2
| +---+ | +----+
| |
| +---+ |
-----| K |<-------
+---+
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 (from y (to K)
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(s) = --------------------------
(z - 0.36788)(z - 0.13534)
+----+
-------------------->| W1 |---> v1
z | +----+
----|-------------+ || T || => min.
| | vz infty
| +---+ v +----+
*--->| G |--->O--*-->| W2 |---> v2
| +---+ | +----+
| |
| +---+ |
-----| K |<-------
+---+
W1 and W2 are the robustness and performancs 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] + 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
- Anderson and Moore, Optimal Control: Linear Quadratic Methods,
Prentice-Hall, 1990, pp. 56-58
- 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
Discrete time control per Zhou, Doyle, and Glover, ROBUST AND OPTIMAL
CONTROL, Prentice-Hall, 1996
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
- Doyle, Glover, Khargonekar, Francis, "State Space Solutions
to Standard H2 and Hinf Control Problems," IEEE TAC August 1989
- Maciejowksi, J.M., "Multivariable feedback design,"
Addison-Wesley, 1989, ISBN 0-201-18243-2
- 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: [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
Do not attempt to use this at home; no argument checking performed.
- 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).
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 : 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.
- @deftypefn{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 = 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 )
-
Octave Controls toolbox demo: Root Locus demo
- 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 = 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 = susball ( inputs )
-
- 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)
- k
-
real scalar (leading coefficient)
Outputs
A, B, C, D
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 : [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 should appear in conjugate pairs
- k
-
real scalar (leading coefficient)
[num,den] = zp2tf(zer,pol,k)
forms the transfer function
num/den
from the vectors of poles and zeros.
% DO NOT EDIT! Generated automatically by munge-texi.
Go to the first, previous, next, last section, table of contents.