Mathematica 8 is now available

Advanced Numerical Methods
An Industrial Application: Controlling the Drum Boiler

Design and analysis of a control system for the drum boiler illustrates a practical application of concepts presented in Advanced Numerical Methods. First, the stability of the model is investigated by computing its poles and time-domain responses. After verifying the controllability of the model, a family of LQR controllers is designed. The characteristics of each LQR controller are plotted to identify the controller that provides a good compromise between improving the damping factor of the model and the size of the feedback gains. Next a controller is designed via the feedback stabilization method using as guidelines the closed-loop poles of the already-designed LQR controller. A state estimator is then constructed and the estimated states are used to implement the controller.

The State-Space Model of a Drum Boiler

Make sure the application is loaded.

<< NumericalControlMethods`

Load the collection of test examples.

<< NumericalControlMethods`Examples`

This state-space model of a drum boiler has a maximum steam flow of about 350 t/h and a drum pressure of 140 bar around an operating point of 90% full load.

ss = DrumBoiler[]

( -0.129` 0.` 0.0396` 0.025` 0.0191` 0.` 0.00139` ) _ •^÷‚ 0.0329` 0.` -0.000078` -0.000122` -0.621` 0.` 0.0000359` 0.0718` 0.` -0.1` 0.000887` -3.851` 0.` -0.00989` 0.0411` 0.` 0.` -0.0822` 0.` 0.0000249` 0.` 0.000361` 0.` 0.000035` 0.0000426` -0.0743` 0.` -5.43`*^-6 1.` 0.` 0.` 0.` 0.` 0 0 0.` 1.` 0.` 0.` 0.` 0 0

The states of the system represent, respectively, the drum pressure (bar), the drum liquid level (m), the drum liquid temperature (degrees Celsius), the riser wall temperature (degrees Celsius), and the steam quality (percent). The inputs are the heat flow to the risers (kJ/s) and the feedwater flow (kg/s). The measured quantities (output of the system) are the drum pressure (bar) and the drum liquid level (m).

System Responses, Stability, and Poles

Start with a quick review of the system’s time-domain behaviors. First the open-loop poles of the system are computed.

These are the open-loop poles of the system.

poles = Eigenvalues[ss [[ 1 ]]]

{-0.18027421301419977`, -0.08581592307994242`, -0.05970493195292894` + 0.017152436138336342` i, -0.05970493195292894` - 0.017152436138336342` i, 0.`}

This is a measure of the damping factor of the system.

Max[Re[poles]]

0.`

The system is marginally stable with damping factor α = 0. Therefore, the input disturbances would not necessarily decay and, at resonance, the output response may become unbounded. This is verified by computing the time responses.

Oscillations in the heat flow to the risers near the nominal value are used to compute the output responses of the open-loop system that is initially at rest. The oscillation frequency corresponds to one of the stable poles.

u0 = {Sin[Im[poles [[ 3 ]]] t], 0}

{Sin[0.017152436138336342` t], 0}

These are the resulting variations in the drum pressure (solid line) and in the drum liquid level (dashed line). Note that the disturbances do not decay over time.

SimulationPlot[ss, u0, {t, 1000}, PlotStyle -> {Thickness[.001], Dashing[{.015}]}] ;

[Graphics:HTMLFiles/anmexamples_13.gif]

The resonance occurs when the system is excited with a signal whose frequency corresponds to the marginally stable eigenvalue, which in the system is λ = 0. Therefore, it has no imaginary part and the resonance frequency is 0 (a constant disturbance).

This adds a constant 0.01 kg/s disturbance to the feedwater flow in the previous input signal.

u = u0 + {0, 0.01} ;

The system resonates and the drum liquid level grows without a bound.

SimulationPlot[ss, u, {t, 1000}, PlotStyle -> {Thickness[.001], Dashing[{.015}]}, PlotRange -> All] ;

[Graphics:HTMLFiles/anmexamples_17.gif]

Obviously, such behavior of the open-loop system is undesirable and needs to be improved by a proper control loop.

Testing the Controllability

To design a feedback controller, the controllability of the system is explored using the function ControllerHessenbergForm.

This is the controller-Hessenberg form of the system.

NumberForm[ControllerHessenbergForm[ss], 3]

( ) _ •^÷‚ -0.0822 0.00572 -0.0407 0.000193 0.0000713 -0.0000249 0 0.0026 -0.118 0.0668 0.0475 -3.81 0 -0.00999 -0.0249 0.0412 -0.114 0.00684 -0.514 0 0 0 0.0052 -0.0311 -0.00779 0.638 0 0 -6 0 0 9.56 × 10 0.000788 -0.0634 0 0 0. -0.139 0.99 -0.00471 -0.00173 0 0 0. -0.00359 -0.00528 -1. -0.0124 0 0

The system is controllable with the tolerance 10^(-4).

Controllable[ss, ControllabilityTest -> FullRankControllerHessenbergBlocks, Tolerance -> 10^(-4)]

True

The system, however, is not controllable if a greater tolerance must be attained.

Controllable[ss, ControllabilityTest -> FullRankControllerHessenbergBlocks, Tolerance -> 10^(-3)]

False

In general, a system is uncontrollable with tolerance tol if there exists a perturbation with a norm less than or equal to tol that makes the perturbed system uncontrollable. Larger feedback gains are necessary to stabilize a controllable system that is close to an uncontrollable state. In the example, the system becomes uncontrollable for relatively small tolerances somewhere between 10^(-3) and 10^(-4).

The Design of the LQR Controller

Since the damping of the open-loop system is low, a controller that increases the damping of the system is needed. With Advanced Numerical Methods, the feedback controller can be designed using several different approaches, such as solving the LQR problem, using constrained feedback stabilization, or implementing the pole assignment approach. This section illustrates the LQR design. For the experiment, Q is a symmetrical matrix constructed from the output matrix of the system and R is an identity matrix multiplied by a scalar coefficient ρ. Evaluation of how different choices of the coefficient ρ affect the performance of the closed-loop system completes the process.

These are the weighting matrices needed for the LQR design.

q = Transpose[ss [[ 3 ]]] . ss [[ 3 ]]

( 1.` 0.` 0.` 0.` 0.` ) 0.` 1.` 0.` 0.` 0.` 0.` 0.` 0.` 0.` 0.` 0.` 0.` 0.` 0.` 0.` 0.` 0.` 0.` 0.` 0.`

r = ρ IdentityMatrix[CountInputs[ss]]

( ρ 0 ) 0 ρ

The LQR gains with control weight ρ range from ρ = 10^(-6) (cheap control) to ρ = 10^2 (minimum energy control). The underlying Riccati equation is solved using Newton method.

rhos = Table[10^i, {i, -6, 2, 0.5}] ;

lqcontrollers = (LQRegulatorGains[ss, q, r /. ρ -> #1, SolveMethod -> Newton] &) /@ rhos ;

This displays the poles of the closed-loop system with varying ρ.

Block[{$DisplayFunction = Identity, len = Length[lqcontrollers]}, Table[ListPlot[ {Re[#], Im[#]} & /@ Flatten[Eigenvalues [ss [[ 1 ]] - ss [[ 2 ]] . lqcontrollers [[ r ]]]], PlotStyle -> {PointSize[0.02], Hue[r/len]}, AxesLabel -> {TraditionalForm[σ], TraditionalForm[j ω]}, PlotRange -> {{-0.28, 0.01}, Automatic} ], {r, len}] ] // Show ;

[Graphics:HTMLFiles/anmexamples_41.gif]

Make sure LogLinearListPlot is available.

<< Graphics`Graphics`

Varying ρ changes the damping factor of the closed-loop system. Increasing ρ improves the stability of the closed-loop system. However, changes in the damping factor are not significant when ρ is either very low or very high.

LogLinearListPlot[Transpose[{rhos, (Max[Re[Eigenvalues[ss [[ 1 ]] - ss [[ 2 ]] . #1]]] &) /@ lqcontrollers}], AxesLabel -> {"ρ", "α(ρ)"}, PlotJoined -> True, Frame -> True] ;

[Graphics:HTMLFiles/anmexamples_45.gif]

This graph shows how varying ρ changes the norm of the feedback gain matrix. The LQR gains are high when ρ is small.

LogLogListPlot[Transpose[{rhos, (MatrixNorm[#1, 2] &) /@ lqcontrollers}], AxesLabel -> {"ρ", {|k(ρ) |} }, PlotJoined -> True] ; 2

[Graphics:HTMLFiles/anmexamples_49.gif]

These are the closed-loop system output responses for various feedback matrices. For convenience, the original open-loop system is represented as a closed-loop system with a dummy (zero) feedback.

responses = (OutputResponse[StateFeedbackConnect [ss, #1], u, {t, 1000}] &) /@ Append[lqcontrollers, Table [0, {CountInputs[ss]}, {CountStates[ss]}]] ;

In the following plots, the closed-loop responses are represented by colored graphs and the open-loop responses are represented by thick black graphs.

len = Length[lqcontrollers] ;

pstyles = Append[Table[Hue[c/len], {c, 1, len}], Thickness[0.008]] ;

These are the output responses of the original and closed-loop systems. From the spectrum, the designer can choose a feedback that brings about suitable values of the drum pressure and the drum liquid level.

Plot[Evaluate[First /@ responses], {t, 0, 1000}, PlotRange -> All, PlotStyle -> pstyles, PlotLabel -> "Drum Pressure"] ;

[Graphics:HTMLFiles/anmexamples_54.gif]

Plot[Evaluate[Last /@ responses], {t, 0, 1000}, PlotRange -> All, PlotStyle -> pstyles, PlotLabel -> "Drum Liquid Level"] ;

[Graphics:HTMLFiles/anmexamples_56.gif]

•The Controller Design Using Constrained Feedback Stabilization

Selected for further experiments is a controller in the middle of the list lqcontrollers that adequately improves the damping factor and yet avoids extremely large feedback gains. In effect, this choice provides a reasonable compromise between the settling time of the system and the sensitivity of the feedback to the noise.

This selects the controller in the middle of the list.

klqr = lqcontrollers [[ Floor[Length[rhos]/2] ]]

( 0.024840806243955807` 0.06796975271227104` 0.00896799734324236` 0.00653443064662637` -0.9358005978063891` ) -6.106959080860623` -17.782664202237918` -2.555784938349549` -1.8732233605264705` 260.3470658531158`

This is a norm of the feedback gains.

MatrixNorm[klqr]

288.66569743509035`

A measure of the noise sensitivity of the feedback is computed.

MatrixNorm[klqr] * MatrixConditionNumber[Eigenvectors[ss [[ 1 ]] - ss [[ 2 ]] . klqr]]

149047.14543712992`

These are the closed-loop poles for the selected feedback matrix.

clplqr = Eigenvalues[ss [[ 1 ]] - ss [[ 2 ]] . klqr]

{-0.18022749559736578`, -0.08604316794507905`, -0.06203239646626457` + 0.018219976182751164` i, -0.06203239646626457` - 0.018219976182751164` i, -0.009900663937787085`}

This is a measure of the damping factor of the closed-loop system.

Max[Re[%]]

-0.009900663937787085`

This finds a state feedback controller such that the closed-loop system will have the same damping factor as the system obtained by the LQR design.

k = StateFeedbackGains[ss, DampingFactorRegion[%], Method -> PartialLyapunovShift]

( 0.017729626054166343` 0.057791767287260036` 0.006687083100384839` 0.004951020844482244` -0.8250606816599704` ) -5.024374243238092` -16.37752912228426` -1.8950432450926347` -1.4030629598589786` 233.81280718363902`

This feedback gain norm is smaller than the one provided by the LQR design with a comparable damping factor.

MatrixNorm[k]

258.512816754113`

The sensitivity of the feedback is also smaller.

MatrixNorm[k] * MatrixConditionNumber[Eigenvectors[ss [[ 1 ]] - ss [[ 2 ]] . k]]

137235.4988149718`

Here are the closed-loop poles with the computed feedback.

clp = Eigenvalues[ss [[ 1 ]] - ss [[ 2 ]] . k]

{-0.1802742130141999`, -0.08581592307994229`, -0.05970493195292887` + 0.017152436138336388` i, -0.05970493195292887` - 0.017152436138336388` i, -0.009900663937787005`}

•The Observer Design

Any feedback design requires knowledge of state variables. Thus it is necessary to estimate the state variables (design an observer) using the inputs, measured outputs, and the system matrices. The first step in designing an observer is to determine whether the system is observable. This is accomplished by using the function ObserverHessenbergForm.

This is the observer-Hessenberg form of the system.

ObserverHessenbergForm[ss]

( -0.129` 0 0.05057637788533299` 0 0 0.` -0.00139` ) _ •^÷‚ 0.0329` 0 -0.23463995043111743` 0.5749653160232213` 0 0.` -0.0000359` 0.07666968775010903` 0 -1.2303169305738535` 2.7848885802533174` 0.008331789876273336` -0.000012308117465654321` 0.0077456656522175696` 0.03088008204748434` 0 -0.47161523265418165` 1.0616929251031983` 0.0033808286682340503` -5.017586611998325`*^-6 0.0031537539816038312` -0.0035720909900808276` 0 0.7835170132658195` -1.9004947307542248` -0.08787599452934525` -0.000021055737200173025` -0.005279264050270166` -1.000000000000071` 0 0 0 0 0 0 0 -1.000000000000071` 0 0 0 0 0

Since the superdiagonal blocks of the observer-Hessenberg form are nonsingular, the system is observable.

Observable[ss, ObservabilityTest -> FullRankObserverHessenbergBlocks, Tolerance -> 10^(-3)]

True

The first two outputs of the system permit direct measurement of two state variables in the system. Thus a reduced-order state estimator of order 3 is constructed. The state estimator poles are chosen in such a way that the state estimator performance is no worse than that of the closed-loop system with feedback obtained earlier.

This chooses the poles of the state estimator.

observerpoles = {-0.06, -0.05, -0.011} ;

The reduced-order state estimator of order 3 is constructed using the block version of the recursive bidiagonal scheme.

observer = ReducedOrderEstimator[ss, observerpoles, Method -> RecursiveBlockBidiagonal] ;

NumberForm[Chop[observer], 3]

( -8 ) _ •^÷‚ -0.06 0.0389 0 9.7 × 10 -0.0000381 0 0 -7 0 -0.0389 0.000432 1.45 × 10 -0.0000449 -0.000265 0.000168 -7 -6 0 0.717 -0.0221 1.01 × 10 5.95 × 10 0.00683 0.0108 0 0 0 0 0 1. 0 0 0 0 0 0 0 1. -18.1 -17.4 42.5 0 0 -7.29 -41.4 124. 124. -50.9 0 0 6.42 46.9 0.603 -0.397 -0.00777 0 0 -0.000393 0.00774

The observer is connected to the system. An implementation of the general approach, shown in the figure below, includes a dummy identity block I that allows construction of the closed-loop system using only the elementary interconnecting functions.

[Graphics:HTMLFiles/anmexamples_85.gif]

The reduced-order estimator in the feedback loop.

The system has two inputs.

p = CountInputs[ss]

2

This connects the all-pass block in parallel to the system ss.

ss1 = ParallelConnect[TransferFunction[IdentityMatrix[p]], ss, Table[{i, i}, {i, p}], {}]

( -0.129` 0.` 0.0396` 0.025` 0.0191` 0.` 0.00139` ) _ •^÷‚ 0.0329` 0.` -0.000078` -0.000122` -0.621` 0.` 0.0000359` 0.0718` 0.` -0.1` 0.000887` -3.851` 0.` -0.00989` 0.0411` 0.` 0.` -0.0822` 0.` 0.0000249` 0.` 0.000361` 0.` 0.000035` 0.0000426` -0.0743` 0.` -5.43`*^-6 0 0 0 0 0 1 0 0 0 0 0 0 0 1 1.` 0.` 0.` 0.` 0.` 0 0 0.` 1.` 0.` 0.` 0.` 0 0

The controller consists of the serially connected observer and the feedback gains block.

controller = SeriesConnect[observer, TransferFunction[k]] ;

NumberForm[%, 3]

( -8 ) _ •^÷‚ -0.06 0.0389 0. 9.7 × 10 -0.0000381 0. 0. -7 0. -0.0389 0.000432 1.45 × 10 -0.0000449 -0.000265 0.000168 -7 -6 0. 0.717 -0.0221 1.01 × 10 5.95 × 10 0.00683 0.0108 -0.00221 0.824 0.0388 0. 0. 0.00109 0.00709 0.626 -233. -11. 0. 0. -0.309 -2.01

This closes the feedback loop and deletes the outputs of the dummy system.

css = DeleteSubsystem[FeedbackConnect[ss1, controller], {}, Range[p]] ;

The eigenvalues of the combined system are the union of the eigenvalues of the closed-loop system and the state estimator. This property is known as separation property.

The separation property is verified.

Sort[Eigenvalues[css [[ 1 ]]]] - Sort[Join[Eigenvalues[ss [[ 1 ]] - ss [[ 2 ]] . k], Eigenvalues[observer [[ 1 ]]]]] // Chop

{0, 0, 0, 0, 0, 0, 0, 0}

Finally, compare the performance of closed-loop systems that use exact and approximated knowledge of the state variables. The approximation is obtained by the reduced-order state estimator.

This is the Bode plot for the second input of the closed-loop systems. The exact and reduced-order designs are represented as solid and dashed lines, correspondingly. The systems possess comparable frequency responses.

DisplayTogetherGraphicsArray[BodePlot[Subsystem [StateFeedbackConnect[ss, k], {2}]], BodePlot[Subsystem[css, {2}], PlotStyle -> {{Dashing[{.02}]}, {Dashing[{.01, .025}], Thickness[.01]}}]] ;

[Graphics:HTMLFiles/anmexamples_97.gif]

This chooses the random initial conditions.

x0 = Table[Random[], {CountStates[ss]}] ;

These are the output responses of the open-loop system for the input signal and the above initial conditions.

SimulationPlot[ss, u, {t, 500}, InitialConditions -> x0, PlotRange -> All] ;

[Graphics:HTMLFiles/anmexamples_100.gif]

These are the output responses of the closed-loop system with the exact (solid lines) and reduced-order state estimator (dashed lines) control law designs.

DisplayTogether[SimulationPlot[StateFeedbackConnect[ss, k], u, {t, 500}, InitialConditions -> x0, PlotRange -> All], SimulationPlot[css, u, {t, 500}, InitialConditions -> Join[x0, Table[0, {CountStates[observer]}]], PlotRange -> All, PlotStyle -> {{Dashing[{.02}]}, {Dashing[{.01, .025}]}}]] ;

[Graphics:HTMLFiles/anmexamples_102.gif]

These plots show that the responses obtained by using exact states and estimated states match closely as time increases (that is, as t approaches infinity).