Analysis and Design of Feedback Control Systems Using Mathematica

Pradipto Ghosh

Wolfram Research, Inc.

Talk Outline

System Model Representations

Operations on Models

Analysis of Model Properties

Design of Control Systems

Model Representations

Transfer function models and state-space models

Tachometer

A very simple system : A tachometer, which is a device used for measuring angular speed in terms of electric voltage

tfTacho = TransferFunctionModel[s,s]

Electric Network

A slightly more complicated system, the state-space model of an electrical network. The SystemsModelLabels option labels the input, output and state variables:

ssNtwrk = StateSpaceModel[{(({{-((Subscript[R, 1])/(Subscript[L, 1])), 0, -(1/(Subscript[L, 1]))}, {0, -((Subscript[R, 2])/(Subscript[L, 2])), (1/(Subscript[L, 2]))}, {(1/C), -(1/C), 0}})),(({{(1/(Subscript[L, 1]))}, {0}, {0}})),(({{0, (Subscript[R, 2]), 0}}))},SystemsModelLabels->{{"e"},{"y"},{"(Subscript[i, 1])","(Subscript[i, 2])","(Subscript[e, c])"}}]

Model Representations

DC Motor

The state-space model of a separately excited DC motor

ssMotor = StateSpaceModel[{(({{0, 1, 0}, {0, -(f/J), ((Subscript[k, 1])/J)}, {0, -((Subscript[k, 2])/L), -(R/L)}})),(({{0}, {0}, {(1/L)}})),(({{1, 0, 0}}))}]

The corresponding transfer-function description of the same system can be obtained by converting the state-space model object to a transfer-function object

tfMotor = TransferFunctionModel[ssMotor]//Simplify

Operations on Models

B-747 Longitudinal Dynamics

Linearized longitudinal equations of motion flying near sea level. Elevator and thrust are control inputs, wind velocity is the disturbance input, pitch rate is
the measured output.

(OverDot[x])(t)=a x(t)+(Subscript[b, u]){{(Subscript[δ, e])(t)},{(Subscript[δ, T])(t)}} +(Subscript[b, w]){{(Subscript[u, w])(t)},{(Subscript[w, w])(t)}} 
y(t)=c x(t)

a=(({{-0.021, 0.122, 0, -0.322, 1}, {-0.209, -0.530, 2.21, 0, -0.044}, {0.017, -0.164, -0.412, 0, 0.544}, {0, 0, 1, 0, 0}, {0, 0, 0, 0, -0.25}}));(Subscript[b, u])=(({{0.01, 0}, {-0.064, 0}, {-0.378, 0}, {0, 0}, {0, 0.25}}));(Subscript[b, w])=(({{0.021, -0.122}, {0.209, 0.530}, {-0.017, 0.164}, {0, 0}, {0, 0}}));
c=(({{0, 0, 1, 0, 0}}));
ssb= StateSpaceModel[{a,ArrayFlatten[{{(Subscript[b, u]),(Subscript[b, w])}}],c}]

Discretization

Discretize the model with a sampling period of τ = 0.1 using method "ZeroOrderHold"

ToDiscreteTimeModel[ssb,0.1,Method-> "ZeroOrderHold"]

Transfer-function representation

Obtain the transer-function matrix of the discrete model

TransferFunctionModel[%]//Simplify

Model Extraction

Extract the subsystem associated with the first two inputs

SystemsModelExtract[%,{1,2}]

Model Connections

Thyroid Hormone Control System

Model built up from simpler TransferFunctionModel subsystems

Quit[]

tfHypPit= TransferFunctionModel[(2/(s+2)),s];
tfThyroid = TransferFunctionModel[(1/(s+a)),s];
tfMetabolism = TransferFunctionModel[10,s];

Interconnection

Equivalent system from the reference and the medication inputs to the output

tfry= SystemsModelFeedbackConnect[SystemsModelSeriesConnect[tfHypPit,SystemsModelSeriesConnect[tfThyroid,tfMetabolism]]]

tfdy=SystemsModelFeedbackConnect[tfMetabolism,SystemsModelSeriesConnect[tfHypPit,tfThyroid]]

Model Connections

A More Complicated Model

Simplify an interconnected system of blocks constituting a motor control system:

Quit[];

SystemsModelSeriesConnect[SystemsModelFeedbackConnect[
SystemsModelSeriesConnect[SystemsModelSeriesConnect[ TransferFunctionModel[(Subscript[k, a]),s],SystemsModelSeriesConnect[ TransferFunctionModel[(Subscript[k, r]),s],SystemsModelFeedbackConnect[SystemsModelSeriesConnect[TransferFunctionModel[((1/(Subscript[R, a]))/((Subscript[T, 1])s+1)),s],SystemsModelSeriesConnect[TransferFunctionModel[(Subscript[k, m]),s],TransferFunctionModel[((1/c)/((Subscript[T, 2])s+1)),s]]],TransferFunctionModel[(Subscript[k, b]),s]]]], TransferFunctionModel[(1/s),s]],
TransferFunctionModel[(Subscript[c, 2]),s]],TransferFunctionModel[(Subscript[c, 1]),s]]

Linearization

Start with non-linear systems, and find linear approximations before proceeding with linear control analysis and design.

Obtain a state-space model by linearizing the two-body equations of motion for a low-thrust spacecraft orbit transfer:

r'=(Subscript[v, r])
θ'=((Subscript[v, t])/r)
(Subscript[v, r])'=(((Subscript[v, t])^2)/r)-(1/(r^2))+a Sin[ϕ]
(Subscript[v, t])'=-(((Subscript[v, r]) (Subscript[v, t]))/r)+ a Cos[ϕ]

StateSpaceModel[{r'[t]==(Subscript[v, r])[t], θ'[t]== (Subscript[v, t])[t]/r[t],(Subscript[v, r])'[t]== (Subscript[v, t])[t]*(Subscript[v, t])[t]/r[t]-1/(r[t])^2 + a Sin[ϕ[t]],(Subscript[v, t])'[t]== -(Subscript[v, r])[t] (Subscript[v, t])[t]/r[t]+ a Cos[ϕ[t]]},{{r[t],0.75},{θ[t],Pi/3},{(Subscript[v, r])[t],0},{(Subscript[v, t])[t],0.05}},{{ϕ[t],Pi/4}},{r[t]},t]//Chop

Analysis

Mathematica provides tools for both state-space and analysis and classical analysis. With a model, verify controllablity, observability, absolute and relative stability...

Functions to compute state-space transformations and decompositions, analyze frequency-domain behavior, visualize system responses...

State-space Analysis

A state-space model

ss= StateSpaceModel[{(({{0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}, {-99.65, -153.7, -74.96, -14.6}})),(({{0}, {0}, {0}, {1}})),(({{26, 36, 11, 1}}))}];

Controllable and Observable?

Through@{ControllableModelQ,ObservableModelQ}@ss

Gramians

Compute the controllability and observability gramians

{wc,wo}=Through@{ControllabilityGramian,ObservabilityGramian}@ss;
MatrixForm/@{wc,wo}//Chop

Balanced realization

The realization is not balanced. Compute a balanced realization:

{s1,ib}=InternallyBalancedDecomposition[ss]

State-space Analysis

Jordan Decomposition

Map[MatrixForm@Normal[JordanModelDecomposition[#][[2]]][[1]]&,{ss,ib}]

InternallyBalancedDecomposition has effected a similarity transformation : Jordan decompositions of the original and transformed models are identical

New gramians

After balancing, the new gramians are also identical

{wC,wO}=Through@{ControllabilityGramian,ObservabilityGramian}@ib;
MatrixForm/@{wC,wO}//Chop

Lower Order Model

A model of reduced order by retaining the two most controllable and observable states

ssred=SystemsModelDelete[ib,None,None,{3,4}]

State-space Analysis

Controllability and Observability

The reduced system is controllable and observable

Through@{ControllableModelQ,ObservableModelQ}@ssred

Bode Plot

The Bode plots of the actual and reduced systems are almost coincident in this scale

BodePlot[{ss,ssred},ImageSize-> 250]

State space Transformation

The original system can be retieved from the balanced system by a state-space transformation

ss==StateSpaceTransform[ib,Inverse@s1]//Chop

State-space Analysis

Non-linear system

Objective: analyze stability of a non-linear system based on local linearization.

(Subscript[x, 1])'[t]== -1.1 (Subscript[x, 1])[t]+ (2/π)ArcTan[0.7 π (Subscript[x, 2])[t]]
(Subscript[x, 2])'[t]== -1.1 (Subscript[x, 2])[t]+ (2/π)ArcTan[0.7 π (Subscript[x, 1])[t]]

Equilibrium points

eqpt=FindInstance[{-1.1 (Subscript[x, 1])+ (2/Pi)ArcTan[0.7 Pi (Subscript[x, 2])]==0, (2/Pi)ArcTan[0.7 Pi (Subscript[x, 1])]-1.1(Subscript[x, 2])==0},{(Subscript[x, 1]),(Subscript[x, 2])},Reals,3]

Linear Model

Obtain linear models about equilibrium points :

MapThread[StateSpaceModel[{(Subscript[x, 1])'[t]==-1.1 (Subscript[x, 1])[t]+ (2/Pi)ArcTan[0.7 Pi (Subscript[x, 2])[t]],(Subscript[x, 2])'[t]== (2/Pi)ArcTan[0.7 Pi (Subscript[x, 1])[t]]-1.1 (Subscript[x, 2])[t]},{{(Subscript[x, 1])[t],#1},{(Subscript[x, 2])[t],#2}},{{u[t],0}},{0},t]&,Transpose[{(Subscript[x, 1]),(Subscript[x, 2])}/.eqpt]]

Lyapunov equation

Verify asymptotic stability of the equilibrium points by checking positive-definiteness of the Lyapunov equation solutions :

(Subscript[a, i]).(Subscript[p, i]) + (Subscript[p, i]).(Subscript[a, i])†==c

Map[PositiveDefiniteMatrixQ[LyapunovSolve[Normal[#][[1]],DiagonalMatrix[{-1,-2}]]]&,%]

Classical Analysis

System Model

tfol =  TransferFunctionModel[((2.7778(s^2+0.192s+1.92))/(s(s+1)^2(s^2+0.384s+2.56))),s]

Open-loop Stability

Open loop system is marginally stable, not BIBO stable

Flatten[TransferFunctionPoles[tfol]]//TableForm

Closed-loop Stability

Closed-loop system stable under unity negative feedback since no net encirclement of critical point

Show[{Graphics[Line[{{-1,0},{-2,2}}]],Graphics[{PointSize[Medium],Point[{-1,0}]}],NyquistPlot[tfol]},ImageSize-> 300,Axes-> True,PlotRange-> {{-2.5,0},{-2,2}}]

Classical Analysis

Step response of closed-loop system is highly oscillatory:

Closed-loop Performance

tfcl = SystemsModelFeedbackConnect[tfol];
OutputResponse[tfcl,1,{t,0,500}];
Plot[%,{t,0,500},PlotRange-> All,Frame-> True]

Verify Stability Margins

Bode plot shows poor stability margins

BodePlot[tfol,StabilityMargins-> True,ImageSize-> 300,StabilityMarginsStyle-> Directive[Red,Thick]]

Numerical values:

{gm,pm}={GainMargins[tfol][[1,2]],Map[{#[[1]],((#[[2]])/Degree)}&,PhaseMargins[tfol]][[1,2]]}

Classical Analysis

Proportional Control

Attempt to improve performance by proportional control.

Conditionally stable system. Gain reduction moves dominant poles into LHP, improving stability margins.

Manipulate[RootLocusPlot[k((2.7778(s^2+0.192s+1.92))/(s(s+1)^2(s^2+0.384s+2.56))),{k,0,5},PlotRange->{{-1.2,0.2},{-1.5,1.5}},AspectRatio->1,PlotLabel-> kValue,PoleZeroMarkers->{Automatic,"ParameterValues"->kValue}],{kValue,0,5,0.001}]

Classical Analysis

Closed-loop Simulation

Performance with gain = 0.27 :

tfnew=TransferFunctionModel[0.27((2.7778(s^2+0.192s+1.92))/(s(s+1)^2(s^2+0.384s+2.56))),s];
OutputResponse[SystemsModelFeedbackConnect[tfnew],1,{t,0,500}];
Plot[%,{t,0,60},PlotRange-> All,Frame-> True]

New Stability Margins

The gain and phase margins have also improved to about 4 and 40 degrees respectively:

{gmnew,pmnew}={GainMargins[tfnew][[1,2]],Map[{#[[1]],((#[[2]])/Degree)}&,PhaseMargins[tfnew]][[1,2]]}

Classical Analysis

Sensitivity

Alternately, search for a gain which gives a maximum sensitivity of about 2, signifying robust stability against plant parameter variations.

Manipulate[Show[{NyquistPlot[k((2.7778(s^2+0.192s+1.92))/(s(s+1)^2(s^2+0.384s+2.56))),{0.01,100},PlotRange-> {{-2,0},{-2,0.2}},PlotLabel-> StringJoin["k=",ToString[k],", |S(jω)(Subscript[|, max])= ",ToString[2]]],Graphics[{Red,Circle[{-1,0},0.5]}]}],{k,0.001,1}]

Design of a 2-DOF DC Servo Motor Controller

Objective

Select state-feedback gain κ and reference gain f for the discretized model of a dc motor speed control system (τ = 0.02 s) to obtain:
i) zero steady-state error due to a unit step
ii) damping ratio ~ 0.7
iii) settling time ~ 1 s

Model

Motor model and its discretized state-space representation

tfmotor = TransferFunctionModel[1/((s+1)(s+10)),s];

ssmotor=StateSpaceModel@ToDiscreteTimeModel[tfmotor,0.02,z,Method-> "ZeroOrderHold"]

Compute Pole Locations

Translate performance specifications to desired pole locations in z-domain. With settling time = 0.9 s,

specsToPoles[ζ_,ts_]:=With[{wn=4.6/(ζ ts)},Map[E^(0.02#)&,{-ζ wn + I wn Sqrt[1-ζ^2],-ζ wn - I wn Sqrt[1-ζ^2]}]];
{z1,z2}=specsToPoles[0.7,0.9]

Design of a 2-DOF DC Servo Motor Controller

Verify Controllability

Pole placement possible since the system is controllable :

ControllableModelQ[ssmotor]

State-feedback Controller

State-feedback gain matrix for pole placement :

κ= StateFeedbackGains[ssmotor,{z1,z2}]

Feedforward Gain

Compute the reference gain matrix to ensure perfect tracking: f=(((c.((((Subscript[I, n])-(a-b. κ)))^(-1))b))^(-1))

sskm=SystemsModelStateFeedbackConnect[ssmotor,κ];
{acl,bcl,ccl}= Normal[sskm][[1;;3]];
f=Inverse[ccl.Inverse[(IdentityMatrix[2]-acl)].bcl]

Construct 2 DOF Control System

sscl = SystemsModelSeriesConnect[sskm,StateSpaceModel[{{{0}},{{0}},{{0}},f},SamplingPeriod-> 0.02]]

Design of a 2-DOF DC Motor Servo Controller

Closed-loop response

Simulate the closed-loop behavior. Tracking is achieved.

OutputResponse[sscl,Table[1,{k,0,120}]];
ListPlot[Map[Thread[{(Range[Length[#]]-1)0.02,#}]&,%],Joined->True,InterpolationOrder->0,PlotRange-> All,Frame-> True]

Other Design Goals

Verify overall damping and settling time:

{ζ,ts}=With[{p2=Most[Eigenvalues[Normal[ToContinuousTimeModel@sscl][[1]]]][[2]]},{ζ1, 4.6/(ζ1 wn1)}/.Flatten@Solve[{-ζ1 wn1==Re[p2],wn1 * Sqrt[1-ζ1^2]== Im[p2]},{ζ1, wn1}]]

LQG/LTR Design: Optimal Control of a Remotely Piloted Aircraft

Model

State-space representation of the longitudinal dynamics, with 2 feedback inputs, 2 stochastic inputs and 2 sensor outputs.

a=(({{-0.025, -40, -20, -32}, {0.0001, -2, 1, -0.0007}, {0.01, 12, -2.5, 0.0009}, {0, 0, 1, 0}}));b=(({{3.25, -0.8, 3.25, -0.8}, {-0.2, -0.005, -0.2, -0.005}, {-32, 22, -32, 22}, {0, 0, 0, 0}}));c=(({{0, 1, 0, 0}, {0, 0, 0, 1}}));ss=StateSpaceModel[{a,b,c}];

Objectives

i) Regulate a perturbed flight condition {0,0.1,0,0} to zero within 10 s with |(Subscript[u, i])|≤ (5^°)
ii) Attempt to recover LQR robustness, lost due to Kalman filtering, by selecting noise covariance matrices

Procedure

i) Design pure LQR by selecting weighting matrices {q,r} to achieve desired performance.
ii) Select noise covariance matrices {w,v} = {ρ v, v} to retrieve LQR robustness

System Properties

System is bounded-input-bounded-output unstable,

Eigenvalues[a]//TableForm

but controllable and observable

Through@{ControllableModelQ,ObservableModelQ}@ss

LQG/LTR Design: Optimal Control of a Remotely Piloted Aircraft

Regulator Design for Deterministic System

q= 0.01IdentityMatrix[4];r=IdentityMatrix[2];
κ= LQRegulatorGains[{ss,{1,2}},{q,r}]

The system has been stabilized:

sscl= SystemsModelStateFeedbackConnect[ss,κ,All,{1,2}];
ListPlot[Map[{Re[#],Im[#]}&,Eigenvalues[Normal[sscl][[1]]]],PlotRange-> All,AxesOrigin-> {0,0},PlotMarkers-> {"X",Medium},PlotStyle-> Red]

Meets the design goals:

x=StateResponse[{sscl, {0,0.1,0,0}},0,{t,0,10}];
GraphicsGrid[{{Plot[x,{t,0,10},PlotRange-> All,FrameLabel-> Style["States",Bold,FontSize-> 14],Frame-> True],Plot[-κ.Transpose[x]*180/Pi,{t,0,10},PlotRange-> All,FrameLabel-> Style["Controls",Bold,FontSize-> 14],Frame-> True]}},ImageSize-> 500]

LQG/LTR Design: Optimal Control of a Remotely Piloted Aircraft

Regulator Design for Stochastic System

Design two LQG regulators for two sets of process-noise covariance matrices

v=IdentityMatrix[2];ρ1=5;ρ2=100;

{lqg1,lqg2}=Map[LQGRegulator[{ss,All,{1,2}},{# v,v},{q,r}]&,{ρ1,ρ2}]

Loop-transfer recovery

How does the LQG return ratio -((κ(s (Subscript[I, n])-a+b κ + l c))^(-1))l c (((s (Subscript[I, n])-a))^(-1))b compare with the pure LQR return ration -((κ(s (Subscript[I, n])-a))^(-1)) b for two estimator gains?

sysk= StateSpaceModel[{a,b,-κ}];
{sysKG1,sysKG2}= Map[SystemsModelSeriesConnect[#,SystemsModelExtract[ss,{1,2}]]&,{lqg1,lqg2}]

LQG/LTR Design: Optimal Control of a Remotely Piloted Aircraft

<<PlotLegends`;
p1=SingularValuePlot[sysk,{(10^(-2)),(10^2)},PlotStyle-> Red];
p2=SingularValuePlot[sysKG1,{(10^(-2)),(10^2)},PlotStyle-> Blue];
p3=SingularValuePlot[sysKG2,{(10^(-2)),(10^2)},PlotStyle->Cyan];
ShowLegend[Show[{p1,p2,p3},Frame-> True,PlotLabel-> Style["Return-ratio singular values",Bold,FontSize-> 14]],{Transpose[{{Graphics[{Red,Thickness[.05],Line[{{-2,0},{2,0}}]}],Graphics[{Blue,Thickness[.05],Line[{{-2,0},{2,0}}]}],Graphics[{Cyan,Thickness[.05],Line[{{-2,0},{2,0}}]}]},{Style["LQR",Bold,FontSize-> 14],Style["ρ = 5",Bold,FontSize-> 14],Style["ρ = 100",Bold,FontSize-> 14]}}],LegendShadow->None,LegendSize->{.5,.20},LegendPosition->{-0.80,.05}}]

Design with ρ = 100 offers better loop-transfer recovery, but performance goals not met

x1=StateResponse[{SystemsModelFeedbackConnect[ss,lqg1,{1,2},{1,2}],{0,0.1,0,0,0,0,0,0}},0,{t,0,10}]//Flatten;
x2= StateResponse[{SystemsModelFeedbackConnect[ss,lqg2,{1,2},{1,2}],{0,0.1,0,0,0,0,0,0}},0,{t,0,10}]//Flatten;
With[{KX1=-κ.Transpose[{x1[[5;;8]]}],KX2=-κ.Transpose[{x2[[5;;8]]}],KX=-κ.Transpose[x]*180/Pi},GraphicsGrid[{{Plot[{(KX1*180/Pi)[[1]],(KX2*180/Pi)[[1]],(KX)[[1]]},{t,0,2},PlotRange-> All,PlotStyle-> {Blue,Cyan,Red},PlotLabel-> Style["(Subscript[u, 1])",Bold,FontSize-> 14]],Plot[{(KX1*180/Pi)[[2]],(KX2*180/Pi)[[2]],(KX)[[2]]},{t,0,2},PlotRange-> All,PlotStyle-> {Blue,Cyan,Red},PlotLabel-> Style["(Subscript[u, 2])",Bold,FontSize-> 14]]}},ImageSize-> 700,Frame-> True]]

Conclusion : Design with ρ = 5 recovers pure LQR performance, and provides larger roll - off at higher frequencies (lower sensitivity to measurement noise)

Questions?