Using objects for automatic control

In Sysquake, models of dynamic systems can be represented by low-level data structures (vectors and matrices) or by object, i.e. higher-level structures which hide their detailed implementation. This latter approach has several advantages:

Objects and methods (the functions which act on them) related to control are defined in library lti. They support linear time-invariant systems described by transfer functions or state-space model in the Laplace or z domain.

The next sections are a tutorial which show you how to create and use these objects. More information is available in the LTI reference. You are supposed to know how to make basic computations with Sysquake (the LME tutorial will show you how.)

Note: Commands you can type are written in bold face. The layout used here reproduces what would be obtained in the command window of Sysquake. But you can also type several commands in a command text field for Sysquake Remote, or in a script file.

Invoking LTI

Depending on your environment, you may have to tell Sysquake to load library LTI.

> use lti

Continuous-time transfer functions

To create an object for the transfer function 1/(s^2+2*s+3) and store it in variable G, the quickest way is

> G = tf(1, [1, 2, 3])
G =
continuous-time transfer function
1/(s^2+2s+3)

Arithmetic operators can be used with transfer functions:

> 2 * G + 1
ans =
continuous-time transfer function
(s^2+2s+5)/(s^2+2s+3)

The poles of the transfer function are obtained with function pole:

> pole(G)
ans =
      -1 + 1.4142j
      -1 - 1.4142j

Discrete-time transfer functions

Discrete-time transfer functions in the domain of z, where z is a variable which represents a forward shift of one sample, are obtained by adding the sampling period Ts as a third argument to tf:

> H = tf(0.5, [1,-1.6,0.7], 0.1)
H =
discrete-time transfer function, Ts=0.1
0.5/(z^2-1.6z+0.7)

You can also use arithmetic operators on discrete-time transfer functions.

Another way to obtain a discrete-time transfer function is to sample a continuous-time one. If G is a continuous-time system whose input is driven by a zero-order hold and whose output is sampled with a sampling period Ts=0.1, the whole, including the zero-order hold element and the sampler, can be modeled with the transfer function H2:

> H2 = c2d(G, 0.1)
H2 =
discrete-time transfer function, Ts=0.1
(4.6712e-3z+4.3697e-3)/(z^2-1.7916z+0.8187)

Graphics

Graphics can be obtained in the time domain (impulse and step response), in the frequency domain (Bode (magnitude and phase), Nyquist and Nichols diagrams), and in the complex plane of the Laplace and z variable (zeros/poles and root locus). While low-level commands have a specific version for discrete-time systems which begin with a 'd' (step and dstep for the step response, for instance), commands which accept lti objects recognize the kind of system and act accordingly: dstep is not implemented for LTI objects.

So to plot the step response of G, type

> step(G)

With Sysquake, you have to clear the figure with clf before plotting something new (in Sysquake Remote, each command creates a new figure by default):

> clf

Many plots are associated with scaling options which make their interpretation easier. For instance, the complex plane should be displayed with the same scale along the x and the y axis, so that the unit circle is actually a circle and not an ellipse. In Sysquake, scaling options are given with command scale, before the command which actually produce graphical output (this is so because the graphical commands take advantage of the scaling options to optimize the values where the functions are evaluated).

The most useful scaling options are equal for having the same scale along x and y axis, loglin for a logarithmic scale along x axis and a linear scale along y axis, and logdb for a decibel scale along y axis.

In addition, special grids can be displayed in the background: hgrid (Hall grid) for Nyquist plots, ngrid for Nichols plots, sgrid for the complex plane of the Laplace variable s, and zgrid for the complex plane of the variable z. The default is a minimal version (a circle or a line); for the complete version, add the command plotoption fullgrid. For simple, generic grids, just type plotoption xygrid.

Taking this together, here is how the Nyquist plot of G could be displayed:

> scale equal
> hgrid
> plotoption fullgrid
> nyquist(G)

The following commands can be used: they are shown with the corresponding scale option and grid command which make most sense.

ScaleGridCommandEffect
impulse(G)impulse response
step(G)step response
scale logdbbodemag(G)magnitude of the Bode plot
scale loglinbodephase(G)phase of the Bode plot
scale equalhgridnyquist(G)Nyquist plot
scale lindbngridnichols(G)Nichols plot
scale equalsgrid or zgridpzmap(G)poles/zeros map
scale equalsgrid or zgridrlocus(G)root locus

System connections

Systems can be connected together in several ways. For serial and parallel connections of SISO systems, operators * and + are used, respectively. A feedback connection can be calculated explicitly with simple formulas which depend on which signals we are interested in.

To control G, we want to add a proportional-integral (PI) controller K. The integrator is necessary in order to avoid static error, as we can see by observing that the frequency response of G is finite at null frequency:

> evalfr(G, 0)
G =
  0.3333

The time constant of the integrator Ti can be set to the slowest mode of G:

> pole(G)
ans =
      -1 + 1.4142j
      -1 - 1.4142j
> Ti = -1/real(ans(1))
Ti =
  1

With a unit gain, the PI controller K is 1+1/(Ti*s)

> K = 1 + tf(1, [Ti, 0])
K =
continuous-time transfer function
(s+1)/s

We already know how to display the Nyquist plot of KG:

> scale equal
> hgrid
> nyquist(K*G)

By plotting successively nyquist(2*K*G), nyquist(4*K*G), etc., we quickly find that a gain of 4 gives satisfactory gain and phase margins. We can change K:

> K = 4 * K
K =
continuous-time transfer function
(4s+4)/s

The numeric value of the margins is computed with margin, whose argument is the transfer function of the loop, i.e. K*G:

> (gm, psi)=margin(K * G)
gm =
 -inf  inf
psi =
    1.6704

Margins are defined only for systems stable in closed loop; since margin gives non-empty values, controller K stabilizes G. The gain margin is given as an interval of gains: here the system remains stable for any gain. The phase margin is 1.67 radians, i.e. 95 degrees.

The closed-loop transfer function Gc between the set-point and the output can be computed explicitly:

> Gc = K * G / (1 + K * G)
Gc =
continuous-time transfer function
(s^4+3s^3+5s^2+3s)/(s^6+4s^5+11s^4+15s^3+14s^2+3s)

We notice that a pair of pole/zero at the origin can be canceled. We can expect other such pairs; they come from the way we calculated the closed-loop transfer function, not from hidden dynamics. We can safely use the function minreal to discard them:

> Gc = minreal(Gc)
Gc =
continuous-time transfer function
(s+1)/(s^3+2s^2+4s+1)

Copyright 2002-2014, Calerga.
All rights reserved.