Control System Design using Maple

Maplesoft, a division of Waterloo Maple Inc., 2004

Introduction

This worksheet demonstrates the use of Maple's LinearAlgebra package to model an inverted pendulum on a cart and create a controller for the system.

The following Maple techniques are highlighted:

Modeling a system

Analyzing the model

Developing and simulating a controller

Animating the results

Introduction to Control System Design

Computer-Aided Control Systems Design tools make extensive use of algorithms from linear algebra.  The LinearAlgebra package in Maple, which combines Maple's state-of-the-art symbolics with NAG's fast and efficient numerics, allows for seamless transition through all the phases of system design.

Control Theory - study of dynamical systems with the intention of causing the output to conform to a desired reference value.

Feedback Control  - measures the output and uses this information to influence its value.  

Examples of Feedback Control Systems  

automobile cruise control

airplane auto-pilots

space-telescopes

computer disk-drive read heads

room temperature control

inverted pendulum on a cart - standard introductory control theory example, which demonstrates basic concepts easily.

Initialization

>	restart: interface( warnlevel=0 ): with(LinearAlgebra): with(plots): with(plottools):

1.  Initial Study of the System

Study the system to be controlled.  

Understand the dynamic requirements.

Select sensors and actuators.

Inverted pendulum on a cart

  sensors - measure the angular displacement of the pendulum and the cart displacement .  

  actuator - motor which moves the cart.

  dynamic requirement - the system must return to its original equilibrium position when disturbances are applied.

2.  Model the System

Apply physical laws (i.e. Newton's laws of motion in the case of mechanical systems) to obtain a mathematical model of the system.  This is usually a system of ordinary differential equations (although can be pde's as well).  

= counter-clockwise angular displacement from upright position

= angular velocity

= position of cart

= velocity of cart

= horizontal force applied to the cart

= half-length of pendulum

= vertical (downward) force exerted on the cart by the pendulum

= horizontal force exerted on the pendulum by the cart

Let and to simplify notation.

>	EQ1 := diff(x(t),t) = y(t);

>	EQ2 := diff(theta(t),t) = phi(t);

Let   = position of center of mass of pendulum

>	r := [ x(t)-L*sin(theta(t)), L*cos(theta(t)) ];

The acceleration of the center of mass is then calculated as

>	rpp := eval( map(diff, r, t$2), {EQ1,EQ2} );

Apply F=ma in the horizontal (x-direction) to the pendulum.

>	eq1 := N = expand( m*rpp[1] );

Apply F=ma in the direction perpendicular to the pendulum.

>	eq2 := (P+g)*sin(theta(t))-N*cos(theta(t)) = m*L*diff(phi(t),t)-m*diff(y(t),t);

Apply F=ma in the horizontal direction to the cart gives:

>	eq3 := u(t) - N = M * diff(y(t),t);

Apply to the pendulum, where is the sum of all moments about the pendulum's

center of mass, is the pendulumn's moment of inertia, and is its angular acceleration.

>	eq4 := P*L*sin(theta(t))-N*L*cos(theta(t)) = inert*diff(phi(t),t);

The moment of inertia is given by:

>	inert := 1/3*m*L^2;

>	eqns := solve({eq1,eq2,eq3,eq4},{diff(y(t),t),diff(phi(t),t),P,N});

>	EQ3 := diff(y(t),t) = simplify( eval(diff(y(t),t),eqns) );

>	EQ4 := diff(phi(t),t) = simplify( eval(diff(phi(t),t),eqns) );

The final nonlinear model is

>	NLEQS := Vector([EQ1,EQ2,EQ3,EQ4]);

3.  Simplify the Model

Create a linearized version of the model.  Although this is only valid in a small interval about an equilibrium point, this range of validity is extended by the use of feedback controllers.

In our example, the  linearized model is only a good approximation of the system when the pendulum is close to the upright position.  However, with a good feedback controller in place the pendulum can go well beyond this range and still be adequately controlled.

Time Domain Representations of Linear Systems

The State-Space Description  of a Linear Time-Invariant (LTI) System  is given by

where

is a vector (of length ) representing the input  signals to the system.

is a vector (of length ) representing the output signals of the system.

is a vector (of length ) representing the internal states of the system.

is an x matrix

is an x matrix

is an x matrix

is an x matrix

Note: if   we have a Single-Input, Single-Output (SISO) system.

If neither is 1, we have a Multi-Input, Multi-Output (MIMO) system.

The solution to the state-space equations is given by

where is a linear operator.

Note: this type of integral is called a convolution integral.

We can't solve the integral because we don't know what is (these functions are disturbances to the system which we can be quite arbitrary).  Instead, we design a controller which for work well for a class of (bounded) functions , and then simulate the controller for specific values u(t), (i.e. a step or impulse function).

Our nonlinear model is of the form

>	f := eval(map(rhs,NLEQS),      {y(t)=y,theta(t)=theta,phi(t)=phi,u(t)=u});

Verify that the pendulum/cart at rest in the inverted position is an equilibrium point.

>	eval( f, {theta=0,y=0,phi=0,u=0} );

Set up the state-vector

>	z := < x,theta,y,phi >;

Assign values to the parameters

>	param_values := g=9.8,M=10.,m=2.,L=3.;

Linearize the equations by computing Jacobians.

>	A := Matrix( [ seq( [ seq( eval( diff( f[i], z[j] ), {theta=0, phi=0, param_values} ), j=1..4 ) ], i=1..4 ) ], datatype=float[8] );

>	B := Matrix( eval( map( diff, f, u ), {theta=0, param_values} ), datatype=float[8] );

Since we are interested only in x and theta as the output variables, we have

>	C := Matrix( [ [ 1, 0, 0, 0 ], [ 0, 1, 0, 0 ] ], datatype=float[8] );

Frequency Domain Representations of Linear Systems

Taking Laplace Transforms of the state-space representation yields the Frequency Domain Representation  of the system.

where is a ( matrix) rational function, called the Transfer Function.

Working in the frequency domain is convenient, since the differential equations have been turned into algebraic ones, and the convolution operator turns into simple multiplication.

Also many control systems have input/output signals that are sinusoidal and thus transfer functions can be obtained directly from experimental data (this is known as system identification).

>	G:=map(normal,C.(ScalarMatrix(s,4)-A)^(-1).B,expanded);

4.  Analyze the Model

A system is internally stable if the real parts of the eigenvalues of are in the left half-plane (LHP),

i.e. .  This keeps the matrix exponential bounded.

A system is input-output stable if the poles of the transfer function are in the LHP.  

Roughly speaking, a system is controllable if it can obtain any desired state by an appropriate choice of input signal.  A system is observable if given any input and output signals, you can uniquely determine the state of the system.  Controllability and observability (and the less strict forms stabilizability and detectability) play a central role in control theory.  For example, a controllable and observable state-space representation is a minimal realization (having the smallest number of states) of the transfer function.  In this case, the eigenvalues of A are equivalent to the poles of the transfer function.

Maple Procedures to compute Controllability and Observability Matrices

>

ControllabilityMatrix:=proc(A::Matrix,B::{Matrix,Vector}) # #  Compute the Controllability matrix of (A,B) = [B AB A^2B ... A^(n-1)B] #  where A is an nxn matrix and B is an nxm matrix #   local i,nA,mA,nB,mB,C,G;   G := Matrix(B);   (nA,mA):=op(1,A);   (nB,mB):=op(1,G);   if nA <> mA then      error "first matrix must be square"   elif nB <> nA then      error "input matrices must have same row dimension";   else      G:=LinearAlgebra:-LA_Main:-MatrixMatrixMultiply(                            A,G,'inplace'=false,'outputoptions'=[]);      C := Matrix(nB,nB*mB,'datatype'=rtable_options(G,'datatype'));      C[1..nB,1..mB] := B;      C[1..nB,mB+1..2*mB] := G;      for i from 2 to (nB-1) do         G:=LinearAlgebra:-LA_Main:-MatrixMatrixMultiply(                            A,G,'inplace'=false,'outputoptions'=[]);         C[1..nB,i*mB+1..(i+1)*mB] := G;      od;   fi;   return C; end:

>

ObservabilityMatrix:=proc(A::Matrix,C::{Matrix,Vector}) # #  Compute the Observability matrix of (A,C), # #   [    C     ] #   [    CA    ] #   [   CA^2   ] #   [    .     ] #   [    .     ] #   [    .     ] #   [ CA^(n-1) ] # #  where A is an nxn matrix and C is an mxn matrix #   local i,nA,mA,nC,mC,O,G;   G := Matrix(C);   (nA,mA):=op(1,A);   (nC,mC):=op(1,G);   if nA <> mA then      error "first matrix must be square"   elif mC <> mA then      error "input matrices must have same column dimension";   else      LinearAlgebra:-LA_Main:-MatrixMatrixMultiply(                             G,A,'inplace'=true,'outputoptions'=[]);      O := Matrix(nC*mC,mC,'datatype'=rtable_options(G,'datatype'));      O[1..nC,1..mC] := C;      O[nC+1..2*nC,1..mC] := G;      for i from 2 to (mC-1) do         LinearAlgebra:-LA_Main:-MatrixMatrixMultiply(                             G,A,'inplace'=true,'outputoptions'=[]);         O[i*nC+1..(i+1)*nC,1..mC] := G;      od;   fi;   return O; end:

>	CC := ControllabilityMatrix(A,B);

>	OO := ObservabilityMatrix(A,C);

Since the controllability and observability matrices have full rank, the system is controllable and observable.

>	Rank(CC);  Rank(OO);

The eigenvalues of A are

>	Eigenvalues(A);

The poles of the transfer function are

>	fsolve(denom(G[1,1])=0,s);

The eigenvalues are equal to the poles of the transfer function since the system is controllable and observable.  Since there is an eigenvalue in the left half plane the system is unstable.

5.  Decide on Performance Specifications

We intend to design a controller that will maintain or create stability in the system while at the same time meeting certain performance criteria.  For example we may specify a controller that will track a certain type of input signal, such as a step function, to within a specified error range.  Or perhaps we require an optimal controller that will cover a more broad class of bounded inputs.

6.  Design a Controller

The choice of controller depends on the performance specifications.  Some typical controllers include

Proportional-Integral-Derivative (PID) Controllers

   these are controllers of the form

   (it is not possible to implement an electronic device with transfer function for the derivative term)

Phase Lead Compensators

   these are controllers of the form     where .

Phase Lag Compensators

   these are controllers of the form   where .

Linear Quadratic Regulators (LQR)

   these are optimal controllers which minimize the cost functional

and Controllers

   the goal here is to produce optimal controllers in appropriate normed spaces.

A Maple LQR routine

>

LQR:=proc(A::Matrix(numeric),          B::Matrix(numeric),          Q::Matrix(numeric),          R) local i,n,invR,ham,evals,evects,M1,M2,j,k,K,P,output_type; n:=op([1,1],A);  # RowDimension(A) if type(R,'Matrix(numeric)') then   invR:=LinearAlgebra:-MatrixInverse(R); elif type(R,'numeric') then   invR:=Matrix([[1/R]]); else   error "%-1 parameter must be of type numeric"         " or a matrix with numeric entries",4; fi; output_type := LinearAlgebra:-GetResultDataType(                 rtable_options(A,'datatype'),                 rtable_options(B,'datatype'),                 UseHardwareFloats); # form the Hamiltonian matrix ham := Matrix([[A,-B . invR . LinearAlgebra:-Transpose(B)],               [-Q, - LinearAlgebra:-Transpose(A)]],                'datatype'=output_type); (evals,evects) := LinearAlgebra:-Eigenvectors(ham); M1:=Matrix(n,n); M2:=Matrix(n,n); j:=0; for i from 1 to 2*n do   if Re(evals[i]) < 0 then      j:=j+1;      for k from 1 to n do         M1[k,j]:=evects[k,i];         M2[k,j]:=evects[n+k,i];      od;   fi; od; P:=map(Re, M2 . LinearAlgebra:-MatrixInverse(M1) ); K:=invR . LinearAlgebra:-Transpose(B) . P; end:

Specify the weights on the input variables.

>	Q := DiagonalMatrix([1,0,1,0]); R := 1;

>	K := LQR(A,B,Q,R);

>	controller := u(t) = eval( (-K . z)[1],                   {x=x(t),y=y(t),theta=theta(t),phi=phi(t)} );

7.  Simulate the Controlled System

Compute the performance of the design, using the original nonlinear model.  This step is repeated with different input signals and if the necessary the controller is re-designed until the performance criteria are met.

>	NLEQS2 := eval( NLEQS, {param_values,controller} ); #NLEQS2 := eval( NLEQS, {param_values,u(t)=0} ); # no controller

Specify disturbance to the system as initial conditions to the model.

>	ICS := x(0)=0, theta(0)=evalf(Pi/4), y(0)=0, phi(0)=0;

>	g := dsolve( {seq(NLEQS2[i],i=1..4),ICS},             {x(t),theta(t),y(t),phi(t)},             type=numeric, output=listprocedure);

>	gx := eval(x(t),g); gtheta := eval(theta(t),g);

>	gx(0), gtheta(0);

Plot the angular displacement of the pendulum with respect to time.

>	odeplot(g, [t,theta(t)], 0..20, numpoints=150);

Plot the cart displacement with respect to time.

>	odeplot(g, [t,x(t)], 0..20, numpoints=100);

Note that the pendulum returns to its upright position and the cart returns to its initial position, as was intended by the design of the controller.

8.  Implement the Controller

Once an appropriate controller has been designed, it is implemented on the real system (or sometimes first as a prototype).

Maple animation

>

pendplot:=proc(x,t)   local ct,st,wc,ws,a,b,c,d,pt1,pt2,pt3,pt4,pendulum,         wheel1,wheel2,cart,w,r1,r2,y0;   cart:=rectangle([x-1,1],[x+1,0],color=red):   wheel1:=ellipse([x-.6,0],.25,.25,color=black,filled=true):   wheel2:=ellipse([x+.6,0],.25,.25,color=black,filled=true):   w:=.1: r1:=3.8: r2:=.2: y0:=.7:   ct:=cos(t): st:=sin(t):   wc:=w*ct: ws:=w*st:   a:=x-r1*st: b:=y0+r1*ct: c:=x-r2*st: d:=y0-r2*ct:   pt1:=[a+wc,b+ws]: pt2:=[a-wc,b-ws]:   pt3:=[c-wc,d-ws]: pt4:=[c+wc,d+ws]:   pendulum:=polygon([pt1,pt2,pt3,pt4],color=blue):   display([pendulum,wheel1,wheel2,cart],       scaling=constrained,axes=none): end:

>	for i from 1 to 60 do   P||i:=pendplot(gx(i/2),gtheta(i/2)): od: display([seq(P||i,i=1..60)],insequence=true,scaling=constrained);

Legal Notice: The copyright for this application is owned by the author(s). Neither Maplesoft nor the author are responsible for any errors contained within and are not liable for any damages resulting from the use of this material.. This application is intended for non-commercial, non-profit use only. Contact the author for permission if you wish to use this application in for-profit activities.