TSi Controls Version 1.1

Getting Started

This notebook illustrates the tools provided in TSi Controls including several new functions contained in Version 1.1 for linear system graphics and functions that implement linear quadratic regulator design. TSi Controls consists of three packages: a package of differential geometric functions GeoTools, the linear control package TSiConL, and the nonlinear control package TSiConN. These are contained in the directory TSiCon_1. This directory also contains an examples file named Examples. TSiCon_1 can be placed anywhwere you like, but a search path must be defined for it so that it will be loaded properly. For example, if you place it in the mathmatica packages directory with full path name "c:\math\packages\tsicon_1", then the following command does this.

$Path=Join[$Path,{"c:/math/packages/tsicon_1/"}];

If you don't want to do this every time you want to use the package, then the above command can be placed in the init.m notebook found in the Packages directory. That way the path will be set up every time Mathematica is initialized.

Basics of The Linear Package

What's in the Package

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

<<TSiConL`

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

?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}=Linearize[f,g,x,{0,0},u,{0,0}]

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, covert back to state space by finding a controllable realization:

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

Is the realization observable?

Observable[A2,C2]

Let's construct a transfer function for this realization

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

Controllability, Observability and Pole Placement

Check the controllability and observability of the original linearized model.

Controllable[A1,B1]

Observable[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}]

Another Example: lateral aircraft dynamics

The following model is representaive of th linearizede 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}}

Controllable[AA,BB]

Observable[AA,CC]

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

K[s]

McMillanDegree[K,s]

LaurentSeries[K,s,3]

ControllableRealization[K,s]

ObservableRealization[K,s]

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]

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

Observable[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 frquebcy range .01 to 10 radians per second is obtained by invoking the function:

BodePlot[G,s,.01,10]

Nyquist Plot

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

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

Similarly, the inverse Nyquist plot can be obtained:

InverseNyquistPlot[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=ColorNyquistPlot[G,s,-4.,+4.,.01]

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.

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

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

RootLocusPlot[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 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}=AlgRiccatiEq[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}=AlgRiccatiEq[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}=AlgRiccatiEq[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 functionality of of TSiConN. First load the package and the examples file.

<<TSiConN`
<<Examples`

The functions available in the nonlinear package TSiConN are

?NonlinearControl

Relative Degree and Decoupling

We will work with example 2 in the examples file.

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

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

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

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

u0=control/.{v1->0}

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

Dynamic extension

Example 12 in the examples file does not have well defined relative degree. Thus, direct computation leads to a singular decoupling matrix. however, dynamic extension does produce a decouplable system.

DynamicExtension[f12,g12,h12,var12]

Adaptive Control

We will illustrate adaptive tracking control for example 14 in the examples file. 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]