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
ComputerAided Control Systems Design tools make extensive use of algorithms from linear algebra. The LinearAlgebra package in Maple, which combines Maple's stateoftheart 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 autopilots
spacetelescopes
computer diskdrive 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).
= counterclockwise angular displacement from upright position
= angular velocity
= position of cart
= velocity of cart
= horizontal force applied to the cart
= halflength 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 (xdirection) 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:
> 
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 StateSpace Description of a Linear TimeInvariant (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 SingleInput, SingleOutput (SISO) system.
If neither is 1, we have a MultiInput, MultiOutput (MIMO) system.
The solution to the statespace 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 statevector
> 
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 statespace 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 halfplane (LHP),
i.e.
. This keeps the matrix exponential
bounded.
A system is inputoutput 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 statespace 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^(n1)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 (nB1) 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^(n1) ]
#
# 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 (mC1) 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.
The eigenvalues of A are
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
ProportionalIntegralDerivative (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; 
> 
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 redesigned 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); 
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([x1,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:=xr1*st: b:=y0+r1*ct: c:=xr2*st: d:=y0r2*ct:
pt1:=[a+wc,b+ws]: pt2:=[awc,bws]:
pt3:=[cwc,dws]: 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
Pi:=pendplot(gx(i/2),gtheta(i/2)):
od:
display([seq(Pi,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 noncommercial, nonprofit use only. Contact the author for permission if you wish to use this application in forprofit activities.