Towards incremental control structure optimisation for process control

Towards Incremental Control Structure Optimisation for Process Control

Wolfgang Birk

Abstract— This paper deals with the optimisation of control structures for multivariable process control systems. Incremen- tal optimisation of control structures is the task of improving the performance of an existing control structure by incrementally introducing or removing dynamic connections in the control structure. A preliminary method for the detection of an im- proving structural change in an existing control structure is presented. The method is applied to the control of a pulverized coal injection process to show its validity and has proven to give good results.

I. I NTRODUCTION

Large scale industrial process plants are characterized by a high degree of interaction between process variables. These interactions are caused by the physical plant construction and also by the additional connections that are created by the plants control systems. The latter is usually provided by the plant vendors and tuned during commissioning of the plant.

Interaction in multivariable control system has been inves- tigated since the 1960s. First measures were the Rijnsdorp interaction measure [1] and the relative gain array (RGA) [2]. These measures are derived from steady-state gain information of the process and indicated how a decentralized controller for a process should be structured, i.e. which controlled/manipulated variables pairs should chosen. Later the relative gains concept was extended to a dynamic mea- sure [3] and also to a measures that can consider richer controller structures than the decentralized control structure, [4], [5]. Apart from the relative gains concept, gramian based measures for interaction and controller structure selection were suggested recently, [6], [7].

A drawback of the above interaction measures is that stability of the closed loop system is not guaranteed. On the contrary, the Niederlinski index [8] gives a steady state stability condition for decentralized controllers. Later the µ- interaction measure [9] was introduced as a tool providing information on closed loop stability and performance loss due to control structure selection. Recently, [10] used the difference between the complementary sensitivity functions of the ideal and the actual closed loop system to provide a performance based structure selection process. Still, the given methods do not provide insight for more richer control structures.

Although, the above methods enable the design of new control structures, which in many cases requires subsequent

W. Birk is with the Control Engineering Group at the Department of Computer Science and Electrical Engineering, Lule˚a University of Technol- ogy, SE-971 87 Lule˚a, Sweden; wolfgang.birk@ltu.se

Fig. 1. 1-DOF multivariable control system

design of controllers, an existing control structure is not considered. From the perspective of the process industry, i.e. steel and papermaking, control systems are rarely re- designed. Even if the complete control system is upgraded, the implemented controllers usually remain unchanged. This is also valid, if physical changes in the plant occur, inten- tional or unintentional. Such changes can lead to deterio- rated plant performance which can not be compensated for by tuning of the existing control scheme. Only structural changes in the control system with subsequent controller design and tuning can restore plant performance or improve it. At the same time, plant owner acceptance for such changes is usually increased if the changes are minor.

An alternative to re-design of the current control structure, is its improvement, where structural changes are introduced incrementally. An incremental change of the control structure means that a dynamic connection in the control system is introduced or removed. The aim of this paper is to determine if a structural change which would yield a performance improvement can be detected for a given control structure, i.e. a decentralized control structure.

II. M ULTIVARIABLE C ONTROL S YSTEMS

In this paper a 1-DOF multivariable control system with process G(s) and controller K(s) is considered, see Fig. 1.

There G(s) has n actuators and m measurement outputs.

Thus, G(s) and K(s) can be given in transfer function matrix notation:

G(s) =







G 11 (s) · · · G 1n (s) .. . . .. .. . G m1 (s) · · · G mn (s)







K(s) =







K 11 (s) · · · K 1 m (s) .. . . .. .. . K n1 (s) · · · K nm (s)







For the sake of simplicity, the Laplace operator s will be

dropped in the sequel. Furthermore, G ij denotes the transfer

function from input u j to output y i and similarly, K ji

New Orleans, LA, USA, Dec. 12-14, 2007

Fig. 2. Simplified multivariable control system for a 2 × 2 process under decentralized control

(a) (b)

Fig. 3. Effects of interaction in multivariable systems

denotes the transfer function from control error e i to the actuator u j . The index j and i refer to the corresponding elements in the vectors u, y and e.

For the closed loop system the complementary sensitivity transfer matrix T and the sensitivity transfer matrix S are defined as

T = (I + GK) ⁻¹ GK (1)

S = I − T = (I + GK) ⁻¹ (2)

Both T and S are transfer matrices in the closed loop system depicted in Fig. 1 and reflect the relative sensitivity of the closed loop system to changes in G. The following relationships can be easily verified

y = T r − T n + SGd I + Sd O (3)

e = Sr − Sn − SGd ^I − Sd ^O (4)

u = KSr − KSn − KSGd ^I − KSd ^O (5) When the multivariable system is under decentralized control, i.e. the controller K is diagonal and the process G is square, the cause of interaction are the off-diagonal elements in the process G. Fig. 2 shows a simplified block diagram for a 2 × 2 process. Suppose, the controllers K ¹¹ and K 22 are designed to achieve a certain behavior of y 1 and y 2 , respectively, and the elements G 12 and G 21 are neglected during the design. Then the channels (e 1 , y 1 ) and (e 2 , y 2 ) are used to specify the performance of the control system and the design yields two sub-systems A and B.

Depending on the structure of G, the sub-systems are connected with each other and interaction of the sub-systems

occurs. Then, two case have to be considered:

1) There is a path from A to B, see Fig. 3a.

Then A acts as a disturbance on B, but not vice versa.

Usually, B can be designed so that the effect of A is attenuated, see [11].

2) There exists a loop from A via B back to A, see Fig. 3b.

Then the dynamics of A are altered, which can have an effect on the stability of the closed loop system. This effect can first be studied after A and B are designed, which complicates controller synthesis.

It can be verified that the process transfer function matrix has a triangular structure if there are no loops. Then the sequential loop closing method [12], [13] yields good design results for decentralized controller structures.

III. O PTIMISATION OF CONTROL STRUCTURES

One of the most common control structures in the industry is the decentralized control structure. It is the one of least complexity and apprehended as the least vulnerable structure.

Additionally, classic control theory can be applied for the design and tuning of controllers, and the implementation is facilitated by the availability of pre-defined building blocks in current process control systems.

Despite these obvious advantages, decentralized control neglects large parts of the dynamics of interconnected sys- tems and therefore can not reach the maximum achiev- able performance. The maximum achievable performance is bounded by the fundamental limitation.

When the performance of a decentralized control system needs to be improved, structural changes that increase the controller complexity slightly can increase the achievable performance for the control system. This increase of achiev- able performance comes with a higher development cost for the control algorithm. It is desirable that this trade-off is quantified, as the optimal control structure can be found using a cost function that considers performance increase and development costs. Thus, the optimization of control structures needs to consider the initial control structure, the performance and the development costs.

In a first step towards incremental optimisation of control structures it is focused on the detection of structural changes that yield an increase in the achievable performance. The development costs will not be considered. For the sequel the performance is assessed and required through properties of the complementary transfer matrix T . For the required performance m scalar transfer functions are defined, one for each measurement output. These m transfer function will be collected into a diagonal transfer matrix T m . A quantification of the performance is derived from T and T m using following three properties:

• Bandwidth ω b defined by σ(T (jω b )) = σ(T (0))/ √ 2

• Steady state tracking: [σ(T (0)); σ(T (0))]

• Peak performance: max ω (σ(T (jω))

Clearly, these properties of T are assumed to be as close as

possible to the required properties of T m and perfect steady

state tracking is achieved if σ(T (0)) = σ(T (0)) = 1.

IV. S TRUCTURAL MISMATCH

When decentralized control, i.e. K is diagonal, is applied to a multivariable process G the controller design is based on the model

G ˆ , diag(G ¹¹ , . . . , G nn )

It implies that the off-diagonal elements of G are small or zero, and thus have little effect on the closed loop performance of the control system. The structural mismatch between G and ˆ G is given by

G ˜ , G − ˆ G

The sensitivity function matrices ˆ T and ˆ S are a function of K and ˆ G. Clearly, ˆ T and ˆ S coincide with T and S only if ˜ G = 0. Moreover, ˆ T and ˆ S can be interpreted as the performance requirements for the closed loop system and the controller K is designed according to them.

Consequently, the mismatch between ˆ T and T is caused by ˜ G or by uncertainties that perturb ˆ G. If ˆ G is perturbed by uncertainties then ˆ T is also perturbed.

Firstly, it is assumed that there are no uncertainties in G. Then, the sensitivity function matrices T and S of the ˆ closed loop system can be derived in terms of ˆ S, ˆ T and the mismatch ˜ G.

Proposition 1: The sensitivity function matrix S can be derived in terms of ˆ S and a multiplicative uncertainty as

S = (I + ∆) ⁻¹ S ˆ where ∆ = ˆ S ˜ GK. Proof:

S = (I + GK) ⁻¹

= ³

I + ˜ GK + ˆ GK ´ −1

= n³

I + ˆ GK ´ ³ ˆ S + ˆ S ˜ GK + ˆ T ´o −1

= ³

I + ˆ S ˜ GK ´ −1

S ˆ

= (I + ∆) ⁻¹ S ˆ

Proposition 2: The complementary sensitivity function matrix T can be expressed in terms of ˆ T using the uncertainty term ∆ = ˆ S ˜ GK

T = (I + ∆) ⁻¹ ³

∆ + ˆ T ´

(6) Proof: The proof follows immediately from T = I − S and ˆ T = I − ˆ S.

Clearly, ∆ gives a description for the mismatch between T and ˆ T , and is denoted as the structural mismatch. Moreover,

∆ has a nice structural property, which is due to the special structure of ˆ S, ˜ G and K. Since both ˆ S and K are diagonal transfer function matrices and ˜ G has zero diagonal elements, it can be shown that ∆ has the same structure as ˜ G. The elements of ∆ are given by

[∆] _ij =

½ 0, i = j

S ˆ ii G ij K jj , otherwise (7)

Hence, each off-diagonal element in ∆ is a filtered version of the corresponding off-diagonal element in G

[∆] _ij = G ij

K jj

1 + G ii K ii

, with i 6= j (8) Thus, ∆ quantifies the effect of off-diagonal elements in a decentralized control system. If ∆ is zero, then T coincides with ˆ T . Otherwise, a minimization of the elements in ∆ minimizes the interaction in the closed loop system.

Furthermore, the relationship (6) can be rewritten to

∆S = ˆ T − T (9)

The output y of the desired closed loop system can be ˆ computed from y = ˆ ˆ T r. Consequently, the output of ∆S is the residual y = ˆ ˜ y −y and ∆ represents the transfer function matrix from the control error e to the residual ˜ y. Thus, ∆ can be directly retrieved from process data by applying system identification algorithms.

Until now, uncertainties were not considered. Suppose, ˆ T is uncertain due to modelling errors in ˆ G, then (9) cannot be used to describe ∆. If an upper bound for the uncertainty in ˆ T is given by an additive diagonal uncertainty term ∆ T , then (9) can be modified to

∆S = ˆ T − T + ∆ T (10)

Consequently, the transfer function matrix ˜ T = ˆ T − T contains information on the uncertainty in ˆ T .

V. E VALUATION OF THE MISMATCH

The transfer matrix ∆ quantifies the mismatch between ˆ T and T , and ˜ T contains a description for the uncertainty in T . The magnitude of the diagonal elements can be used as ˆ an indicator for the uncertainty that is present in ˆ T .

Naturally, the following indicator κ can be introduced κ = max

i

°

°

° T ˜ ii

°

°

° _N (11)

where N denotes a system norm. Thereby, the indicator κ considers the overall dynamics of the diagonal elements and can be used as a minimal threshold for the evaluation of the structural mismatch.

When the off-diagonal elements of ∆ are evaluated, then the following can be suggested. A structural mismatch ap- pears as a non-zero ∆ as soon as ˜ G is non-zero. Due to (7), if the element in ˜ G ij has significant impact on the output, then the element ∆ ij has large magnitude. Still, it has to noted that the magnitude of the ∆ is not invariant to input or output signal scaling. Thus, a proper scaling of the inputs and outputs is a pre-requisite.

In terms of system norms this means that if k ˜ G ij k ^N is large then k∆ ^ij k ^N is large. Thus, the element G ij should be considered in the controller design. This motivates the following rule.

Decision rule. The performance of the closed loop system can be improved by considering the dynamics from input j to output i in the controller design, if

k∆ ^ij k N > τ, i 6= j

where N is a system norm and τ a threshold. ¤ The quantification of the possible performance improve- ment of the closed-loop system in terms of the magnitude of ˜ T ij can be used to choose a threshold τ . This unsolved problem has yet to be addressed. However, it is assumed that a lower bound for τ is given by κ, as κ is an upper bound for the uncertainty in ˆ T . Thus, τ can be chosen so that uncertainties and disturbances are not interpreted as a structural mismatch.

For evaluation purposes the following structure array can be defined.

Definition. Let P be an n × m transfer function matrix.

The structure array S ^N (P ) with system norm N is a real n × m matrix, where the elements are defined as

[S ^N (P )] _ij = kP ^ij k N

¤ Now it can be concluded that from inspection of S ^N (∆), the following conclusions can be drawn:

• If κ is large, then the uncertainties in ˆ T have to be reduced by re-designing the individual controllers or by modifying the performance specification ˆ T .

• If k∆ ^ij k ^N > κ then the dynamics of G ij should be modelled and a new dynamic connection in the control system has to be designed. Only the element with the largest magnitude is considered for the first increment.

As stated in [14], [15], the effect of interaction on a control loop is largest in the frequency region where control is active, namely in the cross-over region. Therefore, the evaluation of ∆ and ˜ T should focus on this region. In [16] weighted system norms are suggested as a tool where certain frequency regions should be emphasized in the analysis of system norms. Hence, the structure array is preferably derived with weighted system norms.

It is still an open question which system norm should be preferred in the evaluation of ∆ and ˜ T . Thence, for the sequel, the H ² -norm is chosen as system norm for the structure array.

VI. I DENTIFICATION OF ∆ ^AND T ˜

When control structure evaluation is applied on-line, a model of ∆ and ˜ T can be obtained by applying system identification algorithms to process data. The theory of system identification is well-developed and many methods for model estimation are available, see e.g. [17] and [18].

Moreover, implementations of these algorithms in Matlab are given in [19].

For the identification of ∆, the input signal vector and the output signal vector are chosen as control error e and residual y, respectively. Obviously, this is a closed-loop ˜ identification problem. In contrary to the identification of ∆, the identification data sets for ˜ T have open-loop character, as the input signal vector becomes r instead of e. Both the open and closed-loop identification problem are studied thoroughly and algorithms are available. A further discussion is not pursued here.

Fig. 4. Schematic drawing of a coal injection vessel

VII. C ASE STUDY

In this case study it is investigated if the control of a coal injection vessel can be improved by a structural change in the controller. First, the coal injection process is shortly introduced and decentralized control of the process is pre- sented. Second, the control structure evaluation methodology is applied to a process model and to process data. Finally, the control structure is improved according to the indications and the re-designed closed loop system is then analyzed.

A. The process

The coal injection vessel is a pressurized multivariable tank system which is discussed in [20], [21] and [22]. There, models for the process are derived and successfully used for the design and analysis of multivariable controllers and a gas-leakage detection system.

In Fig. 4 a schematic drawing of a coal injection vessel is given. The process can be described shortly as follows.

During the injection phase of the vessel, the pressure control valve u N and the flow control valve u C are used to release a constant coal flow from the vessel. The coal flow cannot be measured directly. The available measurements are the net weight m T of the vessel, which is the sum of the nitrogen weight m N and the fine coal weight m C , and the pressure in the vessel p. A linear physical model for the coal injection vessel for the injection phase is given by

˙x(t) = 10 ⁻³

· −0.3234 0.3604 0.2128 −0.2963

¸ x(t) +

· −0.4878 0.6816 0.5721 −0.3043

¸

u(t) (12a) y(t) =

· −0.8379 0.6477

−0.0237 −0.0198

¸

x(t) (12b)

where u(t) = [u N (t) u C (t)] ^T and y(t) = [p(t) m T (t)] ^T . It

is important to note, that the model uncertainty is large for

frequencies below 10 ⁻⁴ rad/sec.

0 500 1000 1500 2000 0

10 20 30

pressure (kPa)

time (s)

0 500 1000 1500 2000

−0.5 0 0.5 1

uN (units)

time (s)

0 500 1000 1500 2000

−0.8

−0.6

−0.4

−0.2 0

mass (t)

time (s)

0 500 1000 1500 2000

−0.2 0 0.2 0.4 0.6

uC (units)

time (s)

Fig. 5. Simulation of step responses of the coal injection vessel with a decentralized PI controller. Reference signal (dotted), output of ˆ T (dashed), output of T (solid)

B. Control structure evaluation

The conventional way to control the vessel with a decen- tralized controller is to stabilize the pressure by means of the pressure control valve and to retain a constant flow using the flow control valve. The performance requirements for the closed loop system are: perfect tracking at steady state, a bandwidth of 0.039 rad/sec and a peak performance of at most 3.7 dB.

A decentralized PI controller that should achieve these requirements is given as

K(s) = 1 s

· 0.04s + 0.001 0

0 −(2s + 0.05)

¸

(13) In Fig. 5, a simulation of coal injection vessel with the controller K(s) is displayed. First, a step in both pressure and mass is applied and later only a step in the mass is applied. Clearly, the control system exhibits slight overshoots and retains a zero steady state control error. Furthermore, it can be seen that pressure loop is largely affected by control action in the mass loop. Assessing the achieved performance yields: perfect tracking at steady state, a bandwidth of 0.0026 rad/sec and a peak performance of 38 dB. This performance is far from what was required. Still, it is a worst case performance and does not need to be exhibited by the system in real life.

A Bode magnitude plot of ∆ and ˜ T is depicted in Fig. 6.

Clearly, the diagonal elements in ∆ are zero, which is not the case for ˜ T . It can be observed, that the off-diagonal elements of ˜ T and ∆ coincide for frequencies above 0.03 rad/sec.

Now the structure arrays for both ∆ and ˜ T using the H ² - norm are derived

S ² (∆) =

· 0.0000 8.4029 0.0003 0.0000

¸

S ² ( ˜ T ) =

· 0.0026 8.0261 0.0000 0.0038

¸

In both cases the off-diagonal element (1, 2) is large com- pared to all the other elements. Moreover, the value for κ

Frequency (rad/sec)

Magnitude (dB)

−150

−100

−50 0 50

From: p

To: p

10⁻⁴ 10⁻³ 10⁻² 10⁻¹ 10⁰ 10¹

−150

−100

−50 0 50

To: m

From: m

10⁻⁴ 10⁻³ 10⁻² 10⁻¹ 10⁰ 10¹

Fig. 6. Bode magnitude plot of ∆ (dashed) and ˜ T (solid).

is 0.0038 which is very small. Consequently, the dynamics from u C to p should be considered in the controller design.

Intuitively, opening the flow control valve introduces a gas flow out of the vessel, which discharge the pressure control loop should compensate for.

From the analysis of ˆ S and ˆ T it can be seen, that the cross-over region of the scalar control loops is between 0.01 rad/sec and 0.1 rad/sec. In order to emphasize the evaluation on the cross-over region, the following weighting function is chosen

F (s) = 0.032s ² (s + 0.057)(s + 0.017) (s + 0.195)(s − 0.1) ² (s − 0.01) ² (s + 0.005) Now the structure arrays can be re-computed with the weighted H ² -norm

S ^{2 F} (∆) =

· 0.0000 10.440 0.0001 0.0000

¸

S ² F ( ˜ T ) =

· 0.0026 9.0882 0.0000 0.0037

¸

The usage of the weighted system norms yields the same result as above, namely the off-diagonal element (1, 2) has to be considered. Already from Fig. 6, it can be concluded that all elements have nearly the same break frequencies, namely about 0.03 rad/sec. Thus, the usage of the weighted system norms lead to the same conclusions in this case.

C. Improved control structure

In order to consider the dynamics from the flow control valve to the pressure, a dynamic compensation that considers mass flow controller and process dynamics is introduced.

There are two options, either sequential loop closing is chosen as design method and the two controllers are re- designed or an additional compensator link is introduced.

In the latter, the control structure become triangular:

K ² (s) = 1 s

· 0.04s + 0.001 −(1.44s + 0.036)

0 −(2s + 0.05)

¸

(14)

0 500 1000 1500 2000 0

10 20 30

pressure (kPa)

time (s)

0 500 1000 1500 2000

−0.5 0 0.5 1 1.5

uN (units)

time (s)

0 500 1000 1500 2000

−0.8

−0.6

−0.4

−0.2 0

mass (t)

time (s)

0 500 1000 1500 2000

−0.2 0 0.2 0.4 0.6

uC (units)

time (s)

Fig. 7. Triangular controller versus decentralized PI controller. Reference signal (dotted), triangular controller (dashed-dotted), decentralized controller (solid)

A new simulation with the triangular controller (14) is per- formed and displayed in Fig. 7. Obviously, the performance of the mass loop has not changed and the performance in the pressure loop is improved. A drawback of the improved control structure is a more intensive control action in the pressure control valve. In terms of the performance proper- ties the bandwidth has improved from 0.0026 rad/sec to 0.0054 rad/sec and the peak performance has improved from 38 dB to 27 dB. Still, this performance is far from what is expected.

VIII. C ONCLUSIONS

Based on the structural mismatch between desired control loop design and achieved control loop design a method for the evaluation and incremental improvement of a decentral- ized control structure in multivariable control is derived.

It is shown that the structural mismatch can be described as a transfer function matrix that directly relates to the ignored off-diagonal elements in the process model. The system norm of the elements in the transfer function matrix are arranged in a so-called structure array. An interpretation of the structure array is discussed and a decision rule is given. In a case study it is shown that the method can be used to improve the control structure of a real-life multivariable control systems.

Although the results are very promising, further inves- tigations and development need to be conducted before a framework for the incremental optimization of control structure is derived.

First of all, it has to be investigated how uncertainties in ˆ T , disturbances, non-linearities and the system dimension effect the modelling of ∆ and ˜ T . Especially, when ∆ and ˜ T are derived from process data together with a-priori knowledge of ˜ G. It also has to be studied how different system norms effect the decision process and how reliability is affected.

Additionally, in the current setup the decision rule is not

scaling invariant, which means that the signal scaling of inputs and outputs affects the values in the structure array.

Thus, the method has to be developed further to achieve scaling invariance.

Finally, it has to studied how thresholds for the structure evaluation are related to performance improvements. Such a relationship would facilitate the actual work of process engineer, that need to trade-off performance improvements against development costs.

R EFERENCES

[1] J. Rijnsdorp, “Interaction on two-variable control systems for distilla- tion columns,” Automatica, vol. 1, no. 1, pp. 15–28, 1965.

[2] E. Bristol, “On a new measure of interaction for multivariable process control,” IEEE Transactions on Automatic Control, vol. 11, pp. 133–

134, January 1966.

[3] M. Hovd and S. Skogestad, “Simple frequency-dependent tools for control systems analysis, structure selection and design,” Automatica, vol. 28, no. 5, pp. 989–996, 1992.

[4] V. Manousiouthakis, R. Savage, and Y. Arkun, “Synthesis of decen- tralized process control structures using the concept of block relative gain,” AIChE Journal, vol. 32, no. 6, pp. 991–1003, 1986.

[5] K. E. H¨aggblom, “Control structure analysis by partial relative gains,”

in Proc. of the 36th Conference on Decision & Control, San Diego, USA, pp. 2623–2624, December 1997.

[6] B. Wittenmark and M. E. Salgado, “Hankel-norm based interaction measure for input-output pairing,” in Proc. of the 2002 IFAC World Congress, Barcelona, 2002.

[7] M. E. Salgado and A. Conley, “Mimo interaction measure and con- troller structure selection,” International Journal of Control, vol. 77, no. 4, pp. 367–383, 2004.

[8] A. Niederlinski, “A heuristic approach to the design of linear multi- variable interacting control systems,” Automatica, vol. 7, pp. 691–701, 1971.

[9] P. Grosdidier and M. Morari, “The µ interaction measure,” Ind. Eng.

Chem Res., vol. 26, pp. 1193–1202, 1987.

[10] J. Lee and T. F. Edgar, “Dynamic interaction measures for decentral- ized control of multivariable processes,” Ind. Eng Chem. Res., vol. 43, pp. 283–287, 2004.

[11] A. Medvedev and G. Hillerstr¨om, “An external model control system,”

Control Theory Adv. Tech., vol. 10, no. 4, part 4, pp. 1643–1665, 1995.

[12] J. Maciejowski, Multivariable Feedback Design. Addison-Wesley, 1989.

[13] S. Skogestad and I. Postlethwaite, Multivariable Feedback Control - Analysis and Design. John Wiley & Sons, 1996.

[14] M. Witcher and T. McAvoy, “Interacting control systems: Steady-state and dynamic measurement of interaction,” ISA Transactions, vol. 16, no. 3, pp. 35–41, 1977.

[15] J.-P. Gagnepain and D. Seborg, “Analysis of process interactions with applications to multiloop control system design,” I&EC Process Design & Development, vol. 21, no. 5, 1982.

[16] W. Birk and A. Medvedev, “A note on gramian-based interaction measures.” Paper no. 4 in Wolfgang Birk, ”Industry Applications of Multivariable Control´´, PhD Thesis, Lule˚a University of Technology, September 2002.

[17] T. S¨oderstr¨om and P. Stoica, System Identification. Prentice Hall, 1989.

[18] L. Ljung, System Identification: Theory for the User. Prentice Hall, 2 ed., 1999.

[19] L. Ljung, System Identification Toolbox: For use with Matlab. The MathWorks Inc., 2000.

[20] A. Johansson and A. Medvedev, “Model based leakage detection in a pulverized coal injection vessel,” IEEE Transactions on Control Systems Technology, vol. 7, pp. 675–682, November 1999.

[21] W. Birk and A. Medvedev, “Pressure and flow control of a pulverized coal injection vessel,” IEEE Transactions on Control Systems Tech- nology, vol. 8, pp. 919–929, November 2000.

[22] W. Birk and A. Medvedev, “Sensitivity analysis of an LQ optimal mul-

tivariable controller for a fine coal injection vessel,” IEEE Transactions

on Industry Applications, vol. 36, pp. 871–876, May/June 2000.