Using the Polysys Class

This is a quick demonstration of the capabilities of the @polysys class.

Timothy J. Wheeler
Dept. of Mechanical Engineering
University of California, Berkeley

Creating a polysys object.
Simulating the system.
Converting other objects to polysys objects.
Interconnections.
Discrete-time systems.
Other Utilities

Creating a `polysys` object.

Since the polysys class is built on the polynomial class, we first create some polynomial objects to work with:

pvar x1 x2 u1 u2

The equations of the system are of the form

Define the polynomial objects f and g

mu = -1;
f = [ x2; mu*(1-x1^2)*x2 - x1 ];
g = [x1;x2];

The polynomial objects states and inputs specify the ordering of the variables. For example, specifying states(1)=x1 indicates that f(1) is the time derivative of x1.

states = [x1;x2];
inputs = [];

Finally, the polynomial objects are used to create a polysys object:

vdp = polysys(f,g,states,inputs)

Continuous-time polynomial dynamic system.
States: x1,x2
State transition map is x'=f(x,u) where
  f1 = x2
  f2 = x1^2*x2 - x1 - x2
Output response map is y=g(x,u) where
  g1 = x1
  g2 = x2

Simulating the system.

The system is simulated over for a given time interval using the sim command. Note that the syntax is similar to ode45.

T = 10;
x0 = randn(2,1);
[t,x] = sim(vdp,[0,T],x0);

plot(x(:,1),x(:,2),'k-')
xlabel('x_1')
ylabel('x_2')
title('Trajectory for the Van der Pol oscillator')

Converting other objects to `polysys` objects.

The simplest object that can be "promoted" to a polysys is a double.

gainsys = polysys(rand(2,2))

Static polynomial map.
Inputs: u1,u2
Output response map is y=g(x,u) where
  g1 = 0.54722*u1 + 0.14929*u2
  g2 = 0.13862*u1 + 0.25751*u2

LTI objects can also be converted to polysys objects.

linearsys = rss(2,2,2);
linearpolysys = polysys(linearsys)

Continuous-time polynomial dynamic system.
States: x1,x2
Inputs: u1,u2
State transition map is x'=f(x,u) where
  f1 = -1.4751*u1 + 0.11844*u2 - 1.0515*x1 - 0.097639*x2
  f2 = -0.234*u1 + 0.31481*u2 - 0.097639*x1 - 1.9577*x2
Output response map is y=g(x,u) where
  g1 = 1.4435*x1 + 0.62323*x2
  g2 = -0.99209*u1 + 0.79905*x2

Polynomial objects can also be converted into a "static" polysys objects.

p = x1^2 - x1*x2;
staticsys = polysys(p)

Static polynomial map.
Inputs: u1,u2
Output response map is y=g(x,u) where
  g1 = u1^2 - u1*u2

Interconnections.

Polysys supports most of the same interconnections as the LTI class with the same syntax and the same semantics. Here are some examples:

append(linearpolysys,staticsys)

Continuous-time polynomial dynamic system.
States: x1,x2
Inputs: u1,u2,u3,u4
State transition map is x'=f(x,u) where
  f1 = -1.4751*u1 + 0.11844*u2 - 1.0515*x1 - 0.097639*x2
  f2 = -0.234*u1 + 0.31481*u2 - 0.097639*x1 - 1.9577*x2
Output response map is y=g(x,u) where
  g1 = 1.4435*x1 + 0.62323*x2
  g2 = -0.99209*u1 + 0.79905*x2
  g3 = u3^2 - u3*u4

series(linearpolysys,gainsys)

Continuous-time polynomial dynamic system.
States: x1,x2
Inputs: u1,u2
State transition map is x'=f(x,u) where
  f1 = -1.4751*u1 + 0.11844*u2 - 1.0515*x1 - 0.097639*x2
  f2 = -0.234*u1 + 0.31481*u2 - 0.097639*x1 - 1.9577*x2
Output response map is y=g(x,u) where
  g1 = -0.14811*u1 + 0.78991*x1 + 0.46034*x2
  g2 = -0.25547*u1 + 0.20011*x1 + 0.29216*x2

The methods append, feedback, parallel, and series are used to interconnect polysys objects.

Discrete-time systems.

It is also possible to create discrete-time polysys objects, as follows:

a = 1;
b = 1;
fduff = [ x2; -b*x1 + a*x2 - x2^3 ];
gduff = [ x1; x2 ];

xduff = [ x1; x2];
uduff = [];
Tsample = 1;

duff = polysys(fduff,gduff,xduff,uduff,Tsample)

Discrete-time polynomial dynamic system.
Sampling time: 1
States: x1,x2
State transition map is x(k+1)=f(x(k),u(k)) where
  f1 = x2
  f2 = -x2^3 - x1 + x2
Output response map is y(k)=g(x(k),u(k)) where
  g1 = x1
  g2 = x2

Discrete-time systems are simulated using the command dsim. Note that simulation time points are specified as (0:T), rather than [0,T].

T = 100;
x0 = [.1;.1];
[t,x] = dsim(duff,(0:T),x0);

plot(x(:,1),x(:,2),'k-')
xlabel('x_1')
ylabel('x_2')
title('Trajectory for the Duffing map')

Other Utilities

Polysys object can be linearized at a given point. This syntax returns an SS object:

xe = [1;2];
vdplin = linearize(vdp,xe);
class(vdplin)

ans =
ss

This syntax returns the state-space data of the linearization:

[A,B,C,D] = linearize(vdp);

Check if a polysys object is linear.

islinear(linearpolysys)

ans =
     1

islinear(vdp)

ans =
     0

Subsystems are referenced using the same syntax as LTI objects:

linearpolysys(1,1)

Continuous-time polynomial dynamic system.
States: x1,x2
Inputs: u1
State transition map is x'=f(x,u) where
  f1 = -1.4751*u1 - 1.0515*x1 - 0.097639*x2
  f2 = -0.234*u1 - 0.097639*x1 - 1.9577*x2
Output response map is y=g(x,u) where
  g1 = 1.4435*x1 + 0.62323*x2

We can also get function handles to the system's state transition and output response maps. This is mostly used to build simulation routines that require handles to functions with a certain syntax (i.e., ode45).

[F,G] = function_handle(vdp);

xeval = randn(2,1);
ueval = []; % VDP is autonomous
teval = []; % The time input is just for compatibility with ode solvers
xdot = F(teval,xeval,ueval)

xdot =
   -0.7420
    0.9962

We can multiply polysys objects by scalars or matrices.

M = diag([2,3]);
M*vdp

Continuous-time polynomial dynamic system.
States: x1,x2
State transition map is x'=f(x,u) where
  f1 = x2
  f2 = x1^2*x2 - x1 - x2
Output response map is y=g(x,u) where
  g1 = 2*x1
  g2 = 3*x2

12*vdp

Continuous-time polynomial dynamic system.
States: x1,x2
State transition map is x'=f(x,u) where
  f1 = x2
  f2 = x1^2*x2 - x1 - x2
Output response map is y=g(x,u) where
  g1 = 12*x1
  g2 = 12*x2

linearpolysys*M

Continuous-time polynomial dynamic system.
States: x1,x2
Inputs: u1,u2
State transition map is x'=f(x,u) where
  f1 = -2.9503*u1 + 0.35533*u2 - 1.0515*x1 - 0.097639*x2
  f2 = -0.46801*u1 + 0.94443*u2 - 0.097639*x1 - 1.9577*x2
Output response map is y=g(x,u) where
  g1 = 1.4435*x1 + 0.62323*x2
  g2 = -1.9842*u1 + 0.79905*x2