TSi ProPac: Controls

Harry G. Kwatny

Department of Mechanical Engineering & Mechanics
Drexel University

Gilmer L. Blankenship

Department of Electrical & Computer Engineering
University of Maryland

Getting Started

\(Off[General::spell]; 
Off[General::spell1]; \)

This notebook illustrates the tools provided in Controls including a few new functions contained in Version 3.0. Controls is contained in the directory ProPac. ProPac can be placed anywhwere you like, but a search path must be defined for it. It is most convenient to place ProPac in the AddOns\Applications directiory, i.e, C:\Program Files\Wolfram Research\Mathematica\3.0\AddOns\Applications. Once ProPac is loaded all of the functions in Dynamics as well as other packages included in TSi ProPac will be available for use and loaded as when needed.

In Mathematica 3 ProPac can be loaded with the command

<<ProPac`

With Mathematica 2.2 use the command ProPac`Master`.

Basics of The Linear Package

What's in the Package

We will demonstrate some of the functions in the linear control systems package ControL. First, load ProPac

<< ProPac`

To obtain a list of functions available in the linear package ControlL:

LinearControl

?LinearControl

Linearization

Define a nonlinear system:

f = {(x1^2) Sin[x1] - x2 + u1 u2, -x2 + (x2^3) + u1 + 
          3 u2}

g = {x1 + x2, x2 + (u2^2)}

x = {x1, x2}

u = {u1, u2}

Now, linearize it about the equilibrium point x=(0,0), u=(0,0)

{A1, B1, C1, D1} = 
      TaylorLinearize[f, g, x, {0, 0}, u, {0, 0}]

Poles & Zeros

The functions Poles and Zeros compute, respectively, the system poles and (transmission) zeros from its transfer matrix.  Poles computes the roos of the least common denominator of minors of all orders and Zeros computes the roots of the greatest common multiple of minors of order k=Rank[G[s]]. Hence these are the poles and zeros of a minimal realization of G[s].

G[s_] := {{(s - 1)*(s + 2), 
            0, (s - 1)^2}, {-(s + 1)*(s + 2), (s - 
                  1)*(s + 1), (s - 1)*(s + 1)}}/((s + 
                1)*(s + 2)*(s - 1))
H[s_] := {{s - 1, s}, {-6, 
            s - 2}}/(1.25*(s + 1)*(s + 2))

G[s] // MatrixForm

H[s] // MatrixForm

Poles[G, s]

Poles[H, s]

G[s_] := {{(s - 1)*(s + 2), 
            0, (s - 1)^2}, {-(s + 1)*(s + 2), (s - 
                  1)*(s + 1), (s - 1)*(s + 1)}}/((s + 
                1)*(s + 2)*(s - 1))
H[s_] := {{s - 1, s}, {-6, 
            s - 2}}/(1.25*(s + 1)*(s + 2))
K[s_] := {{s - 1, 4}, {4.5, 2*(s - 1)}}/(s + 2)

G[s] // MatrixForm

H[s] // MatrixForm

K[s] // MatrixForm

Zeros[G, s]

Zeros[H, s]

Zeros[K, s]

State Space to Transfer Function and Back

Convert the state space description to a transfer function:

G1[s_] := StateSpaceToTransferFunction[A1, B1, C1, D1, s]
G1[s]

Now, convert back to state space by finding a controllable realization:

{A2, B2, C2} = ControllableRealization[G1, s]

Is the controllable realization observable?

ObservablePair[A2, C2]

Let's construct a transfer function for this realization

G2[s_] := StateSpaceToTransferFunction[A2, B2, C2, s]
G2[s]

Minimal Realization

The function MinimalRealization can be used to construct a minimal (or irreducible) realization. It always retuns the system parameters in Jordan form.

H[s_] := 
        Together[
          Outer[Times, {1, 1, 1}, {1, 0, 0}]/(s + p)^3 + 
            Outer[Times, {0, 1, 0}, {0, 1, 0}]/(s + p)^3 + 
            Outer[Times, {1, 0, 1}, {1, 0, 0}]/(s + p)^2 + 
            Outer[Times, {2, 5, 2}, {0, 0, 1}]/(s + p)^2 + 
            Outer[Times, {2, 2, 3}, {-1, -1, 0}]/(s + 
                  p)];

H[s] // MatrixForm
{AA, BB, CC, DD} = MinimalRealization[H, s];
AA // MatrixForm
BB // MatrixForm
CC // MatrixForm
DD // MatrixForm

H[s_] := 
        N[Together[
            Outer[Times, {1, 1, 1}, {1, 0, 0}]/(s + p)^3 + 
              Outer[Times, {0, 1, 0}, {0, 1, 0}]/(s + p)^3 + 
              Outer[Times, {1, 0, 1}, {1, 0, 0}]/(s + p)^2 + 
              Outer[Times, {2, 5, 2}, {0, 0, 1}]/(s + p)^2 + 
              Outer[Times, {2, 2, 3}, {-1, -1, 0}]/(s + 
                    p)]];

H[s] // MatrixForm
{AA, BB, CC, DD} = MinimalRealization[H, s];
AA // MatrixForm
BB // MatrixForm
CC // MatrixForm
DD // MatrixForm

G[s_] := {{1/((s + 2)^3 (s + 5)), 
          1/(s + 5)}, {1/(s + 2), 0}}

G[s] // MatrixForm

{AA, BB, CC, DD} = MinimalRealization[G, s];
AA // MatrixForm
BB // MatrixForm
CC // MatrixForm
DD // MatrixForm

G[s_] := 
      Together[{{1 + 1/((s + 2)^3*(s + 5)), 
            1/(s + 5)}, {1/(s + 2), 1}}]

G[s] // MatrixForm

{AA, BB, CC, DD} = MinimalRealization[G, s];
AA // MatrixForm
BB // MatrixForm
CC // MatrixForm
DD // MatrixForm

MinimalRealization can also be used with systems defined in state space form.

{Anew, Bnew, Cnew, Dnew} = 
        MinimalRealization[AA, BB, CC, DD];

Anew // MatrixForm
Bnew // MatrixForm
Cnew // MatrixForm
Dnew // MatrixForm

Controllability, Observability and Pole Placement

Check the controllability and observability of the original linearized model.

ControllablePair[A1, B1]

ObservablePair[A1, C1]

Since the system is controllable, we can place the poles, at say (-2,-2), using state feedback. First, lets do it using both controls

PolePlace[A1, B1, {-2, -2}]

Now, do it using only the second control (port 2):

PolePlace[A1, B1, 2, {-2, -2}]

Kalman Reduction & Minimal (Irreducible) Realizations

KalmanDecomposition - transforms a state space description to a canonical form identifying the four parts: controllable & observable, controllable & not observable, not controllable & observable, not controllable & not observable. An option (IrreducibleRealization->True) returns only the minimal (or irreducible) system, i.e., the controllable & observable part.

AA = {{1, 1, 0}, {0, 1, 0}, {0, 1, 1}};
BB = {{0, 1}, {1, 0}, {0, 1}};
CC = {{1, 1, 1}}; DD = {{0, 0}};

{Anew, Bnew, Cnew, Dnew} = 
        KalmanDecomposition[AA, BB, CC, DD];
Anew // MatrixForm
Bnew // MatrixForm
Cnew // MatrixForm
Dnew // MatrixForm

Anew
Bnew
Cnew
Dnew

{Anew, Bnew, Cnew, Dnew} = 
        KalmanDecomposition[AA, BB, CC, DD, 
          IrreducibleRealization -> True];
Anew // MatrixForm
Bnew // MatrixForm
Cnew // MatrixForm
Dnew // MatrixForm

AA = {{0, -1, 1}, {1, -2, 1}, {0, 1, -1}};
BB = {{1, 0}, {1, 1}, {1, 2}};
CC = {{0, 1, 0}}; DD = {{0, 0}};

{Anew, Bnew, Cnew, Dnew} = 
        KalmanDecomposition[AA, BB, CC, DD];
Anew // MatrixForm
Bnew // MatrixForm
Cnew // MatrixForm
Dnew // MatrixForm

{Anew, Bnew, Cnew, Dnew} = 
        KalmanDecomposition[AA, BB, CC, DD, 
          IrreducibleRealization -> True];
Anew // MatrixForm
Bnew // MatrixForm
Cnew // MatrixForm
Dnew // MatrixForm

AA = {{0, -1, 1}, {1, -2, 1}, {0, 1, -1}};
BB = {{1, 0}, {1, 1}, {1, 2}};
CC = {{0, 1, 0}}; DD = {{0, 0}};

{Anew, Bnew, Cnew, Dnew} = 
        MinimalRealization[AA, BB, CC, DD];
Anew // MatrixForm
Bnew // MatrixForm
Cnew // MatrixForm
Dnew // MatrixForm

Another Example: lateral aircraft dynamics

The following model is representaive of the linearized lateral dynamics of an aircraft:

AA = {{-10, 0, -10, 0}, {0, - .7, 9, 0}, {0, -1, 
          - .7, 0}, {1, 0, 0, 0}}

BB = {{20, 2.8}, {0, -3.13}, {0, 0}, {0, 0}}

CC = {{0, 0, 1, 0}, {0, 0, 0, 1}}

Lets test controllability and observability:

ControllablePair[AA, BB]

ObservablePair[AA, CC]

Now, we can obtain the transfer function.

K[s_] := StateSpaceToTransferFunction[AA, BB, CC, s]

K[s]

We can obtain the McMillan degree and generate a Laurent series.

McMillanDegree[K, s]

LaurentSeries[K, s, 3]

We can construct controllable or observable realizations from this transfer function.

ControllableRealization[K, s]

ObservableRealization[K, s]

Now, let's look at a SISO system whose input is the first control and output is the second output:

K1[s_] := 
      StateSpaceToTransferFunction[AA, 
        BB[[{1, 2, 3, 4}, {1}]], 
        CC[[{2}, {1, 2, 3, 4}]], s]

K1[s]

McMillanDegree[K1, s]

ObservableRealization[K1, s]

ControllablePair[AA, 
      BB[[{1, 2, 3, 4}, {1}]]]

ControllabilityMatrix[AA, 
      BB[[{1, 2, 3, 4}, {1}]]]

ObservablePair[AA, 
      CC[[{2}, {1, 2, 3, 4}]]]

ObservabilityMatrix[AA, 
      CC[[{2}, {1, 2, 3, 4}]]]

Graphics Functions

Introduction

In the following paragraphs we will illustrate the use of graphics functions that produce Bode, Nyquist, and root locus plots, and a surface plot illustrating the singularities of a rational transefer function.

Bode Plot

Now, define a transfer function

G[s_] := (((s + 10) (s + 4))/((s + 5) ((s^2) + 2 s + 2) ((s^2) + 1 s + 1)))

The Bode plot of G over the frquency range .01 to 10 radians per second is obtained by invoking the function:

Bode[G, s,  .01, 10]

Bode can also be used with the PlotPoints option when finer resolution is required.

\(	G[
        s_] := \(\((s + 10)\)\ \((s + 4)\)\)\/\(s\ \((s\^2 +  .1\ s + 2)\)\ \((s\^2 +  .1\ s + 1)\)\)\)
Bode[G, s,  .2, 10, PlotPoints -> 1000]

Nyquist Plot

The Nyquist plot of G over the frequency range -20 to +20 radians per second is obtained from:

G[s_] := (((s + 10) (s + 4))/((s + 5) ((s^2) + 2 s + 2) ((s^2) + 1 s + 1)))

Nyquist[G, s, -20, +20]

Similarly, the inverse Nyquist plot can be obtained:

InverseNyquist[G, s, -20, +20]

We can also get a color version of the same plot. For variety, we choose a frequency range of -4 to +4 radians per second and specify a frequency increment of 0.01.

ND = ColorNyquist[G, s, -4. , +4. ,  .01]

ColorNyquist has several options. PlotRange can be used when the transfer function has poles on the imaginary axis in order to obtain a visible plot. SensitivityPeak, Bandwidth and Margins place circles on the plot that help analyze performanc and stability margins.

G5[s_] :=  .2 (s + 1) (s^2 + Sqrt[9]  + 9)/( 
            s (s^2  +  2 s  +  2) (s^2 + 2  s + 
                  4)) 
ColorNyquist[G5, s, -3, 3,  .02, 
      PlotRange -> {{-1.25,  .25}, {-1.5, 1.5}}, 
      SensitivityPeak -> 1.5, Bandwidth -> True, Margins -> True]

Root Locus

To compute the root locus plot for values of the gain k ranging from 0 to 40 in steps of nominally unit length (40 steps). Actually, the gain is changed in a nonlinear fashion to achieve more uniform spacing of the computed points. The gain can be roughly estimated from the color scale.

RootLocus[G, s, k, 0, +40. , 1. ]

G[s_] := (((s + 10) (s + 4))/(s (s + 5) (s + 2) ((s^2) +  .2 s + 2) ((s^2) + 1 s + 1)))
RootLocus[G, s, k, 0, +40. , 1. ]

The root locus can also be constructed for a range of gains that includes negative values:

RootLocus[G, s, k, -10, +20. , 1. ]

Pole-Zero Plot

A view of the singlarities of a transfer function is obtained from the pole-zero plot. Try the transfer function

G[s_] := 
      (10/((s + 5) ((s^2) + 5 s + 2) ((s^2) + 1 s + 3)))

Now, the plot

PoleZeroPlot[G, s]

The Linear Quadratic Regulator Problem

Introduction

We will illustrate three functions used for solving the linear quadratic regulator design problem. The main tool is a function which solves the algebraic Riccati equation which is then used to obtain the optimal state feedback controller and the optimal stochastic state estimator.

The Lyapunov Equation

?LyapunovEquation

A = {{0, 1}, {-1, -2}}; Q = IdentityMatrix[2];

Eigenvalues[A]

P = LyapunovEquation[A, Q]

Since A has eignenvalues with negative real parts and Q is positive definite, we expact the solution P to be positive definite.

Eigenvalues[P]

LyapunovEquation[A, IdentityMatrix[3]]

The Algebraic Riccati Equation

It is not our intent to duplicate the many efficient numerical solvers for the algebraic Riccati equation. The method of solution provided here gives infinite precision results when the data is provided in infinite precision (assuming the solution exists). Results of this type are interesting only for rather simple systems. And they are time consuming to obtain. Numerical results are produced - and much more quickly - when numerical data is provided. Consider the solution of the equation
	PA + A'P - PUP + Q = 0
with the following data

A = {{0, 0}, {1, -(Sqrt[2])}}; U = {{1, 0}, {0, 0}}; 
    Q = {{0, 0}, {0, 1}};

The result is obtained with

{P, Eigs} = AlgebraicRiccatiEquation[A, Q, U]

The matrix P is the solution of the Riccati equations and the list Eigs contains the eigenvalues of the associated Hamiltonian matrix that have negative real parts (the "closed loop" eigenvalues. Notice that the results are indeed given in infinite precision. Also, a message is received indicating that the Hamiltonian matrix has an incomplete set of eigenvectors so that a more time consuming method must be used. Consider another set of data

A2 = {{0, 1, 0}, {0, 0, 1}, {0, 0, 0}};
Q2 = {{1, 0, 0}, {0, 0, 0}, {0, 0, 0}};
U2 = {{0, 0, 0}, {0, 0, 0}, {0, 0, 1}};

{P, Eigs} = AlgebraicRiccatiEquation[A2, Q2, U2]

In this case, a complete set of eigenvectors does exist so the preferred method is used and no messages appear. Even in this simple case the increase in speed obtained when using finite precision is quite evident. Try this:

{P, Eigs} = AlgebraicRiccatiEquation[N[A2], Q2, U2]

Optimal State Feedback

Consider the linear control system with (A,B) parameter matrices

A = {{0, 0}, {1, -(Sqrt[2])}}; B = {{1}, {0}};

We wish to design a state feedback controller that minimizes a quadratic performance index with control weighting matrix R and state weighting matrix Q:

R = {{1}}; Q = {{0, 0}, {0, 1}};

This is accomplished with the LQR function:

{K, P, Eigs} = LQR[A, B, Q, R]

LQR returns the state feedback gain F, the Riccati solution P and the closed loop eigenvalues, i.e., the eigenvalues of the matrix (A+BK).

Optimal Stochastic State Estimator

When the states are not directly measured, it is necessary to use a state estimator. The function LQE computes the optimal least squares estimator parameters for a system with (white) driving noise and measurement noise. As an example consider a system with state equations

	             x_dot = Ax + Bu + Gw
	
	output equations
	
		 y = Cx + Du + v
												
	where v,w are white Gaussian noise vectors with zero mean and covariance matrices Q, R, respectively. Then the estimator is of the form:
	
		xest_dot = Axest +Bu + L(Cxest + Du - y)
											
	and the matrix parameter L is obtained from the function LQE. For example, consider the system used above with additional parameters:

CC = {{0, 1}}; G = {{1}, {1}}; Q = {{1}}; R = {{1}};

LQE is used to compute L:

{L, S, Eigs} = LQE[A, G, CC, Q, R]

We can get a better look at the results

N[%]

As does LQR, LQE returns three objects, the observer parameter matrix L, the associated Riccati solution S Ä which is the error covariance matrix, and the estimator eigenvalues Eigs Ä the eigenvalues of (A+LC).

Here again, finite precision computation is observed to be substantially faster.

{L, S, Eigs} = LQE[N[A], G, CC, Q, R]

The Nonlinear Control Package

What's in the Package

We will demonstrate some of the functionality of TSiConN. First load the package and the examples file.

Remove[B];

<< ProPac`

The functions available in the nonlinear package are

NonlinearControl

?NonlinearControl

Controllability & Observability

Controllability

?ControlDistribution

A = {{0, 0, -14}, {0, 0, 0}, {0, 0, -1}};
B = {{1, 2, 4}, {0, 2, 0}, {0, 0, 3}};
x = {x1, x2, x3};
f = A . x;
g = B . x;
ControlDistribution[f, g, x] // MatrixForm

ControlDistribution[f, g, x, IncludeDrift -> False] // 
      MatrixForm

?Controllability

Controllability[f, g, x]

Controllability[f, g, x, LocalControllability -> True]

Observability

Example 65 (Section 7.3) from Vidyasagar (shows how control can be used to establish observability)

?Observability

f = {x2, x3, 0};
g = {0, x1, 0};
h1 = x1;
h2 = x2;
var = {x1, x2, x3};

Observability[f, g, h2, var]

Try it without any input:

Observability[f, {{}}, {h2}, var]

ObservabilityCodistribution[f, Transpose[{g}], {h2}, var]

ObservabilityCodistribution[f, {{}}, {h2}, var]

Local Decomposition

?LocalDecomposition

f = {x1 x3 + x2 Exp[x2], x3, x4 - x2 x3, 
          x3^2 + x2 x4 - x2^2 x3};
G = Transpose[{{x1, 1, 0, x3}}]; u = {u1};
h = {x3};
var = {x1, x2, x3, x4};

LocalDecomposition[f, G, h, var, u, Infinity]

f = {x2, x3, 0};
g = {0, x1, 0}; u = {u1}
h = x2;
var = {x1, x2, x3};

LocalDecomposition[f, Transpose[{g}], {h}, var, u, Infinity]

Exact Feedback Linearization

Some nonlinear systems can be exactly linearized by a combination of a state transformation and nonlinear feedback. TSi Controls implements some of the basic constructions of the 'feedback linearization' method. To determine if a system is feedback linearizable, We can use the test function FeedbackLinearizable

As an example, consider the following:

f30 = {theta (x1^3) + x2, x3, 0};
g30 = {0, 0, 1};
var30 = {x1, x2, x3};
FeedbackLinearizable[f30, g30, var30]

Since this system is linearizable, we can proceed to obtain the exact feedback linearizing state transformation with the function SIExactFBL

Trans = SIExactFBL[f30, g30, var30, True]

To obtain the linearizeable system in normal form coordinates, we need to first invert the transformation using InverseTrans and then use TransSystem to obtain the system representation in the new coordinates.

InvTrans = InverseTransformation[var30, {z1, z2, z3}, Trans]

TransformSystem[f30, g30, var30, {z1, z2, z3}, Trans, 
      InvTrans]

In these coordinates the feedback linearizing controller is obvious.

Another illustration is:

f33 = {1 + x1 + x3, 1 + x2, -x3};
g33 = {{1 + x1 + x3, 0}, {1 + x2, 0}, {0, 
            Exp[x3] + x1 x2}};
var33 = {x1, x2, x3};

FeedbackLinearizable[f33, g33, var33]

The fact that the system is not exactly linearizable is not remarkable, but the way in which linearizability fails is interesting. First, let's test controllability:

Controllability[f33, g33, var33]

Now, we test for involutivity of the requisite distributions. There are only two. Let g1,g2 denote the columns of g33:

{g1, g2} = Transpose[g33];
Involutive[{g1, g2}, var33]

Involutive[{g1, g2, Ad[f33, g1, var33], Ad[f33, g2, var33]}, 
      var33]

Input-Output Linearization, Relative Degree & Decoupling

I-O Linearization & Zero Dynamics for SISO Systems

If a system is not exactly feedback linearizable, it may still be linearizable in an input-output sense. In this event, we can find a state transformation and a nonlinear state feedback control such that the input-output behavior is described by a linear dynamical system. When a system is not exactly feedback linearizable, but is linearizable in the input-output sense, it is possible to use a state transformation to reduce the system to a normal form similar to that illustrated above. However, in the this case there remain residual nonlinear dynamics, called the internal dynamics, which are decoupled from the output. Hence the input-output behavior is linear even though the entire state dynamics are not. In the SISO case, the coordinate transformation can always be chosen such that the internal dynamics are independent of the control.

As an example, consider the following (example 4.1.5 in Isidori):

f29 = {x1 x2 - (x1^3), x1, -x3, (x1^2) + x2};
g29 = {0, 2 + 2 x3, 1, 0}; h29 = {x4};
var29 = {x1, x2, x3, x4}; newvar29 = {z1, z2, z3, z4};
{T1, T2} = SISONormalFormTrans[f29, g29, h29, var29]

Trans = Join[T1, T2];
InvTrans = InverseTransformation[var29, newvar29, Trans]

{Newf, Newg} = 
        TransformSystem[f29, g29, var29, newvar29, Trans, 
          InvTrans];

Now, the output y=x4=z1 is identically equal to zero if and only if z1=0 and z2=0. When this obtains, the internal dynamics reduce to the zero (output) dynamics:

ZeroDyn = {(∂_t){z3[t], 
              z4[t]} == (Newf[[{3, 
                  4}]] /. {z1 -> 0, z2 -> 0, 
                z3 -> z3[t], z4 -> z4[t]})}

As another example, consider ( Example 4.1.4 from Isidori)

f28 = {-x1, x1 x2, x2}; g28 = {Exp[x2], 1, 0}; 
    h28 = {x3};
var28 = {x1, x2, x3}; newvar28 = {z1, z2, z3};

{T1, T2} = SISONormalFormTrans[f28, g28, h28, var28]

Trans = Join[T1, T2]

InvTrans = InverseTransformation[var28, newvar28, Trans]

This illustrates a difficulty in using Mathematica's Solve function. Mathematica finds the 'bad equation' to be inconstant with the equation z2=x2. However, the bad equation is indeed equivalent. In fact, the inverse transformation is obvious:

InvTrans = {1 - Exp[z2] + z3, z2, z1}

Now, we can proceed to obtain the transformed system.

{Newf, Newg} = 
      TransformSystem[f28, g28, var28, newvar28, Trans, InvTrans]

ZeroDyn = {(∂_t) z3[
              t] == (Newf[[3]] /. {z1 -> 0, z2 -> 0, z3 -> z3[t]})}

Relative Degree, Decoupling & Zero Dynamics for MIMO Systems

For some SISO systems and most MIMO systems it is difficult or impossible to find complete global transformations to the linearizable normal form. Fortunately, it is not necessary to do so. For SISO systems the difference between the order of the system and the order of the zero dynamics is called the relative degree or relative order. In the above examples

RelativeOrder[f29, g29, h29, var29]

RelativeOrder[f28, g28, h28, var28]

For MIMO systems relative order is a vector (or list) of integers, in a sense characterizing the relative degree associated with each output. The function relative VectorRelativeOrder returns a vector, even for the SISO case.

VectorRelativeOrder[f28, g28, h28, var28]

Recall that the existance of an I-O feedback linearizing controller depends on the nonsingularity of the 'decoupling matrix' which can be systematically computed by differentiation. Consider the following example (adapted from Isidori).

var2 := {x1, x2, x3}
f2 := {0,  x1  +  x2^2,  x1  -  x2}
g2 := {Exp[x2], Exp[x2], 0}
h2 := {x3}
p2 := {Sin[x1], 0, Exp[x2]}

We can use the functions VectorRelativeDegree and then DecouplingMatrix to obtain:

ro = VectorRelativeOrder[f2, g2, h2, var2]

DecouplingMatrix[f2, g2, h2, var2, ro]

Since the decoupling matrix is nonsingular, we can proceed to compute linearizing controller. The function IOLinearize does this.

{rho, alpha, ro, control} = IOLinearize[f2, g2, h2, var2]

It is also possible to compute the partial normal form transformation:

z = NormalCoordinates[f2, g2, h2, var2, ro]

And it is always possible to obtain a local representation of the zero dynamics.

u0 = control /. {v1 -> 0}

LocalZeroDynamics[f2, g2, h2, var2, u0, z]

Dynamic extension

The following example (Slotine, p.269) does not have well defined relative degree. Thus, direct computation leads to a singular decoupling matrix. however, dynamic extension does produce a decouplable system.

f12 := {0, 0, 0};
g12 := Transpose[{{Cos[x3], Sin[x3], 0}, {0, 0, 1}}];
h12 := {x1, x2};
var12 := {x1, x2, x3};
yyd12 := {E^t, t};

DynamicExtension[f12, g12, h12, var12]

Dynamic Inverse

When a system has a singular decoupling matrix, an (equivalent) alternative to dynamic extension is to construct a dynamic inverse. Dynamic inverses are interesting in their own right and have useful application such as in the design of tracking controllers. ProPac implements the Hirschorn-Singh constructions in the function StructureAlgorithm. The following example is from Neijmeijer and Van der Schaft.

x = {x1, x2, x3, x4};
f = {0, x3, x4, 0}; G = {{1, 0}, {x4, 0}, {0, 0}, {0, 1}};
h = {x1, x2};

{DD, CC, HH, ZZ} = StructureAlgorithm[f, G, h, x, t]

The dynamic inverse is defined by f+Gu where:

u = RightInverse[DD] . (-CC + ZZ)

f + G . u

Dynamic extension can also be used to construct the inverse.

IOLinearize[f, G, h, x]

{fnew, Gnew, hnew, xnew} = DynamicExtension[f, G, h, x]

{rho, alpha, r0, u} = IOLinearize[fnew, Gnew, hnew, xnew]

u = (u /. {v1 -> D[y1[t], {t, r0[[1]]}], 
              v2 -> D[y2[t], {t, r0[[2]]}]});
(fnew + Gnew . u) // MatrixForm

Adaptive Control

Adaptive Regulation

We will illustrate adaptive regulation using the following example (adapted from Isidori).

var32 := {x1, x2, x3};
f32 := {0,  x1  +  x2^2,  theta*x1  -  x2};
g32 := {Exp[x2], Exp[x2], 0};
h32 := {x3};
p32 := {Sin[x1], 0, Exp[x2]};

First, we need to specify the desired closed loop pole locations. It is necessary to do this in groups according to the vector relative degree:

ro = VectorRelativeOrder[f32, g32, h32, var32]

Poles = {{-2, -2}};

Now, compute the adaptive control

AdaptiveRegulator[f32, g32, h32, var32, t, {theta}, {}, 
      Poles]

AdaptiveRegulation via Backstepping

TSi Controls implements tool for designing adaptive controllers for SISO systems via backstepping. Backstepping can be applied to systems that are in parametric-strict-feedback-form (PSFF) or parametric-pure-feedback-form (PPFF). Methods exist to test a given system to determine if it can be transformed to PSFF and to determine the required transformation if it can be. The functions PSFFCond and PSFFSolve provide an implementation of these constructions. The function AdaptBackstepReg computes the adaptive regulator for a system in PSFF.

Once again, consider example 30 from above. The function PSFFSolve tests to see if the system can be placed in PSFF and if, so generates the transformed system.

f30 = {theta (x1^3) + x2, x3, 0};
g30 = {0, 0, 1};

Information["PSFFSolve", LongForm -> False]

PSFFSolve[f30, g30, {x1, x2, x3}, {theta}]

In fact, we see that the transformation is an identity and that the system was in PSFF to begin with. Now, we design an adaptive regulator:

Information["AdaptiveBackstepRegulator", 
      LongForm -> False]

{r1, r2, r3} = 
        AdaptiveBackstepRegulator[f30, g30, 
          var30, {theta}, {{AdGain}}, {c1, c2, c3}];

r1

r2

r3

Adaptive Tracking

Adaptive tracking control can also be illustrated using the following example (from Ghanadan and Blankenship).

\(\(f14 := {x2, 
\ \ \ \ \ \ \ \ 0, 
\ \ \ \ \ \ \ \ 0, 
\ \ \ \ \ \ \ \ \((\(-theta\))\)*Sin[x1]^2 + 
            x3*Cos[x1], 
\ \ \ \ \ \ \ \ \((theta/2)\)\ Sin[2\ x1] + \ 
            x3\ Sin[x1]};\)
\)
g14 := Transpose[{{0, 0, 1, 0, 0}, {0, 1, 0, 0, 0}}];
h14 := {x4, x5};
\(\(var14 := {x1, x2, x3, x4, x5};\)
\)

First define a reference signal:

RefSig = 
      Table[ToExpression["yd" <> ToString[i1] <> "[t]"], {i1, 1, 
          Length[h14]}]

Specify the desired closed loop pole locations:

Poles = {{-1, -2}, {-1, -2, -2}};

Now, compute the adaptive controller.

AdaptiveTracking[f14, g14, h14, var14, t, {theta}, RefSig, {}, 
      Poles]

Creating SIMULINK S-Functions

TSi provides a convenient interconnection to MATLAB/SIMULINK that enables a user to incorporate systems or controllers developed in Mathematica as subsystems in SIMULINK. This in turn allows controllers developed using TSi Controls to be easily implemented using MATLAB's Real Time Workshop. The main tool provided in TSi Controls is the function CreateControllerMEX that can be used to generate S-functions for SIMULINK 1.3 (the function CreateMEXFile can be used with the earlier version, SIMULINK 1.2).
	CreateControllerMEX generates C source code that compiles as a 'mex' function in the structure of a SIMULINK S-function. Hence, it is a standard SIMULINK block with inputs and outputs that can be interconnected with any other SIMULINK block using SIMULINK's block diagramming facility. As an illustration consider example 13, which is a 4 state system with 2 inputs and 2 outputs:

f13 := {x2, x1 (x4^2) - ((x4^2)/2), 
          x4, ((4 x2 x4 - 8 x1 x2 x4)/(4 (x1^2) - 4 x1 + 5))};
h13 := {x1, x3};
g13 := 
        Transpose[{{0, 1, 0, 0}, {0, 0, 0, 
              (4/(4 (x1^2) - 4 x1 + 5))}}];
var13 := {x1, x2, x3, x4};
u = {u1, u2};

StateEqns = f13 + g13 . u

The S-function source code, with filename MXTst.c is obtained using CreateControllerMEX:

CreateControllerMEX[var13, u, h13, 
      StateEqns, {"X0"}, {{4, 1}}, "MXTst.c"]

The system for which the S-function representation is to be created can consist of differential equations and an output map as is the case above. The output map can have direct feedthrough of the inputs. In fact, it is not necessary to have any differential equations - the system can be a memoryless. Here is a 2 input - 2 output map that merely interchanges the order of the inputs:

CreateControllerMEX[{}, {u1, u2}, {u2, 
        u1}, {}, {}, {}, "MXTst2.c"]

Bibliography

Effective use of the tools provided in this package require a basic understanding of modern control system analysis and design concepts. Some general references for the methods underlying this package are given below.

Books on linear control abound. Here are just a few titles, including some classics: 
Nise, N., Control Systems Engineering, Benjamin/Cummings, Redwood City, 1995.
Friedland, B. Control System Design, McGraw-Hill, New York, 1986.
Kailath, T., Linear Systems, Prentice-Hall, Englewood Cliffs, 1980.
Wonham, W. M., Linear Multivariable Control: a Geometric Approach, Springer-Verlag, New York, 1979.
Anderson, B. and J. Moore, Linear Optimal Control, Prentice-Hall, Englewood Cliffs, 1971.
Brockett, R. W., Finite Dimensional Linear Systems, John Wiley, New York, 1970.

Nonlinear and adaptive control theory has made import strides in the past decade. Usefull references for the methods implemented here are:
Isidori, A., Nonlinear Control Systems, Springer-Verlag, New York, 1995.
Niejmeijer, H. and A. J. van der Schaft, Nonlinear Dynamical Control Systems, Springer-Verlag, Berlin, 1990.
Slotine, J-J. E. and W. Li, Applied Nonlinear Control, Prentice-Hall, Englewood Cliffs, 1991.
Friedland, B., Advanced Control System Design, Prentice Hall, Englewood Cliffs, 1996.
Astrom, K. and B. Wittenmark, Adaptive Control, Addison-Wesley, Reading, 1995.
Tao, G. and P. Kokotovic, Adaptive Control of Systems with Actuator and Sensor Nonlinearities, John Wiley, New York, 1996.
Young, K-K. D., Variable Structure Control for Robotics and Aerospace Applications, Elsevier, Amsterdam, 1993.