Robust Control Of An Industrial Distillation Column

  • Uploaded by: H.E. Musch
  • 0
  • 0
  • May 2020
  • PDF

This document was uploaded by user and they confirmed that they have the permission to share it. If you are author or own the copyright of this book, please report to us by using this DMCA report form. Report DMCA


Overview

Download & View Robust Control Of An Industrial Distillation Column as PDF for free.

More details

  • Words: 46,153
  • Pages: 227
Diss. ETH No. 10666

Robust Control of an Industrial High-Purity Distillation Column

A dissertation submitted to the SWISS FEDERAL INSTITUTE OF TECHNOLOGY ZURICH for the degree of Doctor of Technical Sciences presented by HANS-EUGEN MUSCH Dipl. Chem.-Ing. ETH born June 19, 1965 citizen of Germany accepted on the recommendation of Prof. M. Steiner, examiner Prof. Dr. D. W. T. Rippin, co-examiner 1994

3

To my grandparents

4

5

Acknowledgments

This Ph. D. thesis was written during my years as a research and educational assistant of the Measurement and Control Laboratory at the Swiss Federal Institute of Technology (ETH) at Zurich. I would like to take this opportunity to thank the numerous persons who have supported this project. First of all I express my gratitude to Prof. M. Steiner. He arranged this project and helped to overcome many difficulties with the industrial environment. Many thanks are also due to him and to Prof. D. W. T. Rippin for the critical examination of this thesis, which essentially improved its clarity. The numerous discussions with my colleagues and their uncountable suggestions gave rise to important contributions to this work. In this context, E. Baumann, U. Christen, and S. Menzi must be specially mentioned. Last but not least I should emphasize the support of B. Rohrbach. She never lost her patience with my never ending questions concerning the English language. Without her willingness to correct the manuscript, the choice of the English language for this thesis would have been impossible.

6

7

Content

Abstract .........................................................................................

15

Kurzfassung ..................................................................................

17

1

Introduction .........................................................................

19

1.1 “Modern Control: Why Don’t We Use It?” .............................

19

1.2 Scope and significance of this thesis ......................................

21

1.2.1 Distillation as a unit operation example ....................

21

1.2.2 Earlier research ...........................................................

21

1.2.3 Robust control and nonlinear plants...........................

22

1.2.4 Contributions of this thesis .........................................

22

1.3 Structure of the dissertation ..................................................

23

1.4 References ...............................................................................

26

2

The Distillation Process — An Industrial Example .......................................................

29

2.1 Introduction ............................................................................

29

2.2 Column design and operation ................................................

29

2.3 Steady-state behavior .............................................................

32

2.4 Composition dynamics ............................................................

35

2.5 Control objectives and configurations ...................................

37

2.5.1 The 5x5 control problem ..............................................

39

2.5.2 Control design steps.....................................................

40

2.6 Tray temperatures as controlled outputs ..............................

42

Pressure-compensated temperatures 42

8

2.6.2 Temperature measurement placement ......................

44

2.7 References ..............................................................................

45

3

A Rigorous Dynamic Model of Distillation Columns ..........................................................

47

3.1 Introduction ............................................................................

47

3.2 Conventions ............................................................................

48

3.3 The objective of modelling .....................................................

48

3.4 Simplifying assumptions .......................................................

48

3.5 Balance equations ..................................................................

51

3.5.1 Material balances ........................................................

51

3.5.2 Energy balance equations ...........................................

52

3.6 Fluid dynamics .......................................................................

55

3.6.1 Liquid flow rates .........................................................

55

3.6.2 Pressure drop...............................................................

57

3.7 Phase equilibrium ..................................................................

59

3.7.1 Vapor phase composition ...........................................

59

3.7.2 Boiling points...............................................................

60

3.8 Volumetric properties ............................................................

60

3.8.1 PVT relations ..............................................................

61

3.8.2 Density .........................................................................

61

3.9 Enthalpies ..............................................................................

62

3.10 Numerical solution ................................................................

63

The dependent variables and the equation system 63 3.10.2 Formal representation of the DAE .............................

66

3.10.3 The index .....................................................................

66

3.10.4 Solution methods and software ..................................

67

3.11 Notation ................................................................................

71

3.12 References ............................................................................

74

9

4

Linear Models ......................................................................

77

4.1 Introduction ............................................................................

77

4.2 How to linearize the rigorous model? ....................................

78

4.2.1 The state, input, and output vectors ..........................

78

4.2.2 Handling of the algebraic equation system ................

80

4.3 Linearization of a simplified nonlinear model ......................

80

4.3.1 The simplified model ...................................................

80

4.3.2 Analytical linearization ...............................................

84

4.4 Linearization of the rigorous model .......................................

86

4.4.1 Model modifications ....................................................

86

4.4.2 Numerical linearization...............................................

88

4.5 Comparison of the linear models ...........................................

89

4.5.1 Open loop simulations ................................................

89

4.5.2 Singular values ............................................................

92

4.6 Order reduction .......................................................................

94

4.7 Summary .................................................................................

96

4.8 Appendix: Model coefficients ..................................................

97

4.9 Notation ...................................................................................

101

4.9.1 Matrices and Vectors ..................................................

101

4.9.2 Scalar values ................................................................

102

4.10 References ..............................................................................

103

5

A Structured Uncertainty Model .....................................

105

5.1 Introduction ............................................................................

105

5.2 Limits of uncertainty models .................................................

106

5.3 Input uncertainty ...................................................................

107

5.4 Model uncertainty ...................................................................

110

5.4.1 Column nonlinearity ...................................................

110

5.4.2 Unmodelled dynamics..................................................

117

10

5.5 Measurement uncertainty .....................................................

118

5.6 Specification of the controller performance ..........................

119

5.7 Summary ................................................................................

120

5.8 References ..............................................................................

122

μ-Optimal Controller Design ...........................................

123

6.1 Introduction ............................................................................

123

6.2 The structured singular value ...............................................

124

6.2.1 Representation of structured uncertainties...............

124

6.2.2 Definition of the structured singular value ...............

126

6.2.3 Robustness of stability and performance ...................

128

6.3 The design model ...................................................................

130

6.4 Controller design with μ-synthesis .......................................

133

6.4.1 Synthesis algorithms ..................................................

134

6.4.2 Applying the DK-Iteration..........................................

137

6.4.3 Applying the μK-Iteration...........................................

137

6.5 Design of controllers with fixed structure ............................

148

6.5.1 Diagonal PI(D) control structures .............................

149

6.5.2 PI(D) control structures with two-way decoupling....

156

6.5.3 PID control structures with one-way decoupling ......

162

6.6 Summary ................................................................................

165

6.7 References ..............................................................................

167

6

7

Controller Design for Unstructured Uncertainty — A Comparison ......................................................................

171

7.1 Introduction ............................................................................

171

7.2 Diagonal PI-control ................................................................

172

7.2.1 The BLT method..........................................................

172

11

7.2.2 Sequential loop closing ................................................

174

7.2.3 Optimized robust diagonal PI-control.........................

176

7.3 PI-control with decoupling .....................................................

178

7.4 H∞ optimal design ...................................................................

184

7.5 Summary .................................................................................

189

7.6 References ...............................................................................

189

8

Feedforward Controller Design .......................................

191

8.1 Introduction ............................................................................

191

8.2 The design problem ................................................................

192

8.2.1 The design objective .....................................................

192

8.2.2 One-step or two-step design?.......................................

192

8.3 H∞-minimization .....................................................................

194

8.4 Optimization approach ...........................................................

198

8.5 Summary .................................................................................

201

8.6 References ...............................................................................

202

9

Practical Experiences ........................................................

205

9.1 Introduction ............................................................................

205

9.2 Controller implementation .....................................................

206

9.3 Composition estimators ..........................................................

209

9.4 Controller performance ..........................................................

210

9.5 Economic aspects ....................................................................

216

9.6 Summary .................................................................................

216

10 Conclusions and Recommendations ...............................................................

219

10.1 Introduction ..........................................................................

219

12

10.2 Controller synthesis .............................................................

220

10.3 State-space or PID control? .................................................

221

10.4 How many temperature measurements? ............................

222

10.5 Column models .....................................................................

223

10.6 Recommendations ................................................................

223

10.6.1 Academic research ......................................................

223

10.6.2 Decentralized control systems ....................................

224

10.6.3 Cooperation industry—university ..............................

225

Curriculum vitae ..................................................................

227

13

Symbols δ

Uncertainty scalar value

Δ

Uncertainty matrix or deviation from nominal operating point

θ

Parameter vector

κ

Condition number, κ = σ max ⁄ σ min

λ

Eigenvalue

μ

Structured singular value

ρ

Spectral radius

σ

Singular value

B

Bottom product stream (mol/s)

D

Distillate stream (mol/s) or diagonal scaling matrix

d

Disturbance signals

e

Control error

F

Feed flow rate (mol/s)

Fl

Lower fractional transformation

G(s)

Transfer function

G u → y Transfer function from control signals u to output signals y I

Identity matrix

K(s)

Controller

L0

Reflux (mol/s)

M

Joint weighted plant and controller, M ( P, K ) = F l ( P, K )

P

Weighted plant

p

Pressure (N/m2)

14

r

Reference signals

Se(s)

Sensitivity function at e, S e ( s ) = [ I + G ( s )K ( s ) ] –1

Su(s)

Sensitivity function at u, S u ( s ) = [ I + K ( s )G ( s ) ] – 1

T

Temperature (°C)

T r → y Transfer function from reference signals to output signals u

Control signals

V51

Boilup (mol/s)

W(s)

Diagonal matrix of weighting transfer functions

w(s)

Weighting transfer function

xD

Top product composition (mol/mol)

xB

Bottom product composition (mol/mol)

xF

Feed composition (mol/mol)

y

Output signals

15

Abstract

It is well known that high-purity distillation columns are difficult to control due to their ill-conditioned and strongly nonlinear behavior. Usually distillation columns are operated within a wide range of feed compositions and flow rates, which makes a control design even more difficult. Nevertheless, a tight control of both product compositions is necessary to guarantee the smallest possible energy consumption, as well as high and uniform product qualities. This thesis discusses a new approach for the dual composition control design, which takes the entire operating range of a distillation column into account. With the example of an industrial binary distillation column, a structured uncertainty model is developed which describes quite well the nonlinear column dynamics with several simultaneous model uncertainties. This uncertainty model forms the basis for feedback controller designs by μ-synthesis or μ-optimization. The resulting controllers are distinguished by a high controller performance and high robustness guaranteed for the entire operating range. This method enables the synthesis of state-space controllers as well as the μ-optimal tuning of advanced PID control structures. The already satisfactory compensation of feed flow disturbances can be improved even further by use of feedforward control. Even for the design of the feedforward controllers the basic ideas of the feedback controller design can be employed. A simultaneous feedforward controller design for two column models representing the extreme column loads yields outstanding results. Similar to the feedback controller design, a design of state-space controllers by H ∞ -minimization or an optimal tuning of simple feedforward control structures by parameter optimization is possible. Control engineers working in an industrial environment are conscious of the high effort needed for the implementation of state-space control-

16

lers in a distributed control system. Therefore a controller design based on PID or advanced PID control structures is of high relevance for the industrial practice. Usually, the performance of these PID control structures is expected to lag significantly behind the performance of highorder state-space controllers. However, comparing the performances of the state-space controllers with those of the advanced PID controllers, merely slight advantages of the state-space controllers are detected. This surprising result, achieved with an unconventional tuning of the PID control structures, allows the simple implementation of advanced PID control structures in a decentralized control system without a significant loss of controller performance. The good robustness properties and the high performance of the control schemes are confirmed by the implementation of an advanced PID control scheme on a real industrial distillation column. An estimation of the economic benefits made by this project much more than justifies the effort expended.

17

Kurzfassung

Bekanntermaßen sind Rektifikationskolonnen mit hohen Produktreinheiten wegen ihres schlecht konditionierten und stark nichtlinearen Verhaltens schwierig zu regeln. Häufig werden sie in einem weiten Bereich unterschiedlicher Zulaufkonzentrationen und -mengen betrieben, was den Entwurf von Regelungen zusätzlich erschwert. Dennoch ist eine gute Regelung beider Produktkonzentrationen notwendig, um einerseits einen möglichst kleinen Energieverbrauch und andererseits hohe und einheitliche Produktqualitäten sicherzustellen. Diese Arbeit beschreibt einen neuen Ansatz für den Entwurf von Konzentrationsregelungen, der den gesamten Arbeitsbereich einer Rektifikationskolonne berücksichtigt. Am Beispiel einer industriellen binären Rektifikationskolonne wird ein strukturiertes Unsicherheitsmodell entwickelt, welches das nichtlineare dynamische Verhalten der Rektifikationskolonne durch mehrere Modell-Unsicherheiten gut beschreibt. Dieses Unsicherheitsmodell bildet die Basis für den Entwurf von Reglern mittels μ-Synthese oder μ-Optimierung. Die resultierenden Regler zeichnen sich durch eine – über den gesamten Betriebsbereich garantierte – hohe Regelqualität bei sehr grosser Robustheit aus. Dieses Vorgehen erlaubt sowohl den Entwurf von Zustandsregelungen als auch die Berechnung μ-optimaler Einstellungen für erweiterte PID-Regelstrukturen. Die bereits zufriedenstellende Unterdrückung von Störungen der Zulaufmenge wird durch den Einsatz einer Störgrößenaufschaltung noch verbessert. Auch für ihren Entwurf können ähnliche Konzepte verwendet werden. Ein Entwurf von Störgrössenaufschaltungen, bei dem gleichzeitig zwei Modelle der Rektifikationskolonne berücksichtigt werden, welche die extremen Kolonnenbelastungen wiedergeben, führt zu hervorragenden Ergebnissen. Vergleichbar mit dem Regelungsentwurf können sowohl Störgrössenaufschaltungen mit der Struktur

18

von Zustandsregelungen (durch Minimierung der H ∞ -Norm) als auch Störgrößenaufschaltungen mit einfacher Struktur (durch Parameteroptimierung im Zeitbereich) berechnet werden. In der industriellen Praxis tätige Regelungstechniker sind sich der Schwierigkeiten, die mit der Realisierung von Zustandsregelungen auf dezentralen Prozeßleitsystemen verbunden sind, sicherlich bewußt. Daher ist der Regelungsentwurf auf der Grundlage von PID- oder erweiterten PID-Regelstrukturen von hoher praktischer Relevanz. Meist bleibt die mit solchen Strukturen erzielbare Regelgüte hinter der von Zustandsregelungen deutlich zurück. In dieser Arbeit werden die entworfenen Zustandsregelungen und die optimal eingestellten fortgeschrittenen PID-Regelstrukturen verglichen. Dabei zeigt sich, daß auch mit einfachen Regelstrukturen, die entsprechenden unkonventionellen Regler-Einstellungen vorausgesetzt, eine Regelqualität erzielt wird, die der von Zustandsregelungen nahekommt. Dieses überraschende Resultat erlaubt die einfache Implementierung von erweiterten PIDRegelstrukturen in dezentralen Prozeßleitsystemen ohne wesentlichen Verlust an Regelgüte. Die Erprobung eines Regelungsentwurfs auf der Grundlage fortgeschrittener PID-Strukturen an der industriellen Rektifikationskolonne bestätigt die große Robustheit und die hohe Regelgüte in der Praxis. Dabei zeigt eine Abschätzung der Wirtschaftlichkeit, daß der bei einem solchen Projekt notwendige Aufwand mehr als gerechtfertigt ist.

1.1 “Modern Control: Why Don’t We Use It?”

19

Chapter 1 Introduction

1.1 “Modern Control: Why Don’t We Use It?” “Modern Control: Why Don’t We Use It?” is the title of a paper written by R. K. Pearson in 1984 [1.4]. In the first section of that paper Pearson states: “Advanced control systems utilizing multivariable strategies based on process models can outperform traditional designs in broad classes of application. Yet, in spite of market forces demanding better process performance and ample evidence showing that the improvements can be achieved, the gap between theory and practice in the industrial sector is not narrowing appreciably.” Ten years later the situation has not changed. The modern control theories provide the process control engineer with increasingly sophisticated tools for a robust, model-based controller design. The advantages of these controllers over the PID control structures which are usually tuned on-line, have been shown in numerous publications. Nevertheless, more than 90% of all control loops in the process industry use PID control, while only a few applications of the modern control theories can be reported [1.10]. Therefore the mismatch between theory and practice is still evident. Some of the reasons for this situation are discussed below.

20

1 Introduction

Distributed Control Systems For a control engineer in the process industry, process control in the first place is a hardware problem. His perspective is the installation and configuration of a Distributed Control System (DCS) [1.1]. Even the modern DCS are often limited to PID and advanced PID control. For the DCS, an implementation of modern state space controllers requires either the coupling with an external computer or the programming of software modules. Both ways are troublesome and expensive. The university research pays little attention to this situation. The design of robust controllers with fixed structures (e.g., PID control structures) is a largely unexplored field. Dynamic Models Linear dynamic models are the foundation of a modern, robust controller design. However, no general dynamic models are available for unit operations. For each plant linear dynamic models must be developed, based on either linearization of nonlinear models or on system identification methods. Both ways are often expensive and very timeconsuming ([1.5], [1.6]). Furthermore, most plants in the process industry show a strongly nonlinear dynamic behavior, which is unsatisfactorily described by a single linear model. Economic benefits The economic benefits of improved control tend to be significantly underestimated. A benchmark study by ICI “indicated that the effective use of improved process control technology could add more than one third to the worldwide ICI Group’s profits” [1.1]. Another study shows smaller, but still massive benefits [1.2]. Of course it is not necessary to replace all PID-controllers by modern advanced control structures. Most control problems in the process industry are handled well with simple PID control. However, strongly nonlinear or/and ill-conditioned plants require advanced control techniques for a high controller performance.

1.2 Scope and significance of this thesis

21

1.2 Scope and significance of this thesis 1.2.1 Distillation as a unit operation example Distillation is one of the most widely used unit operations in the process industry. In the simplest case, a distillation column separates a feed of two components into a top product stream (with a high fraction of the low-boiling component) and a bottom product stream (with a high fraction of the high-boiling component). In an industrial setting, the feed flow rate and the feed composition may vary within a wide range of operating conditions. This separation consumes a huge amount of energy. A minimization of the energy consumption and an economic optimal operation usually require (1) a tight control of both product compositions (dual composition control) and (2) often small fractions of impurities in the product streams (high purity distillation). However, the strongly nonlinear and ill-conditioned behavior makes high-purity distillation columns difficult to control. Therefore high-purity distillation columns have become an interesting test case for robust control design methods. 1.2.2 Earlier research Without any doubt the distillation process is the most studied unit operation in terms of control. Skogestad estimates that new papers in this field appear at a rate of at least 50 each year [1.7]. It is practically impossible to give a review of all these publications. The interested reader is advised to consult the reviews of Tolliver and Waggoner [1.8], Waller [1.9], MacAvoy and Wang [1.3], and the recent review of Skogestad [1.7]. If we focus our interest on the design of linear, time-invariant controllers, we must state that all the well-known model-based and robust control design methods (LQG/LTR, H∞, Normalized Coprime Factorization, μ-synthesis, etc.) have been applied to distillation columns. However, all these publications discuss the controller design for just one operating point. The problem designing a robust controller which maxi-

22

1 Introduction

mizes the controller performance for the entire operating range has not been addressed as yet. 1.2.3 Robust control and nonlinear plants The well-known robust control design methods like H∞ -minimization or LQG/LTR are based on the assumption of an unstructured, frequency dependent uncertainty at one location in the plant. Such an unstructured uncertainty may be a multiplicative uncertainty at plant input or output, or an additive uncertainty. A controller design for the entire operating range of a distillation column using one of these well-known methods has two inherent problems: • Due to the high nonlinearities an estimation of unstructured uncertainty bounds will lead to very large bounds, prohibiting any acceptable controller design. • A controller design using any arbitrary, smaller uncertainty bound guarantees robust performance (RP) and robust stability (RS) for the actual operating point, but not for the entire operating range. 1.2.4 Contributions of this thesis This thesis presents a new approach for the composition control design of a binary distillation column (Figure 1.1). The design concept is based on a structured uncertainty model which describes the column dynamics for the entire operating range quite well. The resulting controller designs using μ-synthesis (for state-space controller) or μ-optimization (for controllers with fixed structure), respectively, lead to results which guarantee robust performance and robust stability for the entire operating range of the distillation column. Special emphasis is placed on the optimal tuning of easy-to-realize PID-control structures. It will be shown that extraordinary controller performance can be achieved even with these relatively simple controller structures.

1.3 Structure of the dissertation

Standard approaches

Linear model for a single operating point

23

Improved approach Uncertainty model describing column dynamics for entire operating range

Robust control design

μ-synthesis

H∞, LQG/LTR, ........

μ-optimization

Weak point:

Advantage:

RS & RP at other

RS & RP guaranteed for whole operating range

operating points

Figure 1.1: Robust control design approaches

1.3 Structure of the dissertation A robust, model-based controller design for a distillation column consists of several steps. A typical course is illustrated in Figure 1.2. The results and methods of each step influence all the following steps. The consideration of just one of these design steps, disengaged from all others, neglects the conceptional coherence. Therefore all of the design steps are discussed within this thesis. The sequence orients itself to the natural course of the controller design.

24

1 Introduction

Nonlinear Model

Linear Models

Uncertainty structure

Parameter adaptations

Controller synthesis

Nonlinear simulations

Tests on plant

Implementation in DCS

Figure 1.2: Steps of a model based controller design

The following chapter consists of three parts: The first part describes the design and operating data of the distillation column, followed by an overview of the steady-state and dynamic column behavior. The second part discusses the control objectives and control configuration for this column, while the third part describes the use of pressure-compensated temperatures as controlled outputs. Rigorous nonlinear dynamic models are the basis for simulation studies and for linearization. They are discussed in Chapter 3.

1.3 Structure of the dissertation

25

The main subject of Chapter 4 is the derivation of linear models. Two different methods are presented which lead to linear models which neglect and include flow dynamics, respectively. A structured uncertainty model which describes the nonlinear behavior of the distillation column for the entire operating range is developed in Chapter 5. Based on that structured uncertainty model, controllers can be designed within the framework of the structured singular values. In the first part of Chapter 6 the theoretical background of the structured singular value μ is summarized. While the second part of that chapter presents the μoptimal design of state-space controllers, the third part is dedicated to the μ-optimal design of PID control structures. Simulation studies confirm the theoretical results. In Chapter 7 the results of the μ-optimal controller design are compared with results obtained by more common design methods, based on an unstructured uncertainty description. Usually the feed flow rate is a measured disturbance input to a distillation column. Therefore, feedforward control can significantly improve the compensation of feed flow disturbances, which is discussed in Chapter 8. A controller design should yield a satisfactory control quality not only in dynamic simulations but also in the real plant. The results of the practical implementation are presented in Chapter 9. The conclusions and the recommendation for further research in Chapter 10 complete this thesis. The literature references and, if necessary, the special notations are given at the end of each chapter.

26

1 Introduction

1.4 References [1.1] Brisk, M.L.: “Process Control: Theories and Profits,” Preprints of the 12th World Congress of the International Federation of Automatic Control, Sydney, July 18-23, 7, 241-250 (1993) [1.2] Marlin, T. E., J. D. Perkins, G. W. Barton, and M. L. Brisk: “Benefits from process control: results of a joint industry−university study,” J. Proc. Cont., 1, 68-83 (1991) [1.3] McAvoy, T. J. and Y. H. Wang, “Survey of Recent Distillation Control Results,” ISA Transactions, 25, 1, 5-21 (1986) [1.4] Pearson, R. K.: “Modern Control: Why Don’t We Use It?,” InTech, 34, 47-49 (1984) [1.5] Schuler, H., F. Algöwer, and E. D. Gilles: “Chemical Process Control: Present Status and Future Needs ⎯ The View from European Industry,” Proceedings of the Fourth International Conference on Chemical Process Control, South Padre Island, Texas, February 17-22, 29-52 (1991) [1.6] Schuler, H.: “Was behindert den praktischen Einsatz moderner regelungstechnischer Methoden in der Prozess-Industrie,” atp, 34, 3, 116-123 (1992) [1.7] Skogestad, S.: “Dynamics and Control of Distillation Columns – a Critical Survey,” Preprints of the 3rd IFAC Symposium on Dynamics and Control of Chemical Reactors, Distillation Columns and Batch Processes, April 26-29, College Park, Maryland, 1-25 (1992) [1.8] Tolliver, T. L. and R. C. Waggoner: “Distillation Column Control; a Review and Perspective from the CPI,” Advances in Instrumentation, 35, 1, 83-106 (1980) [1.9] Waller, K. V.: “University Research on Dual Composition Control of Distillation: A Review”, Chemical Process Control 2, Sea Island, Georgia, January 18-23, 395-412 (1981)

1.4 References

27

[1.10] Yamamoto, S. and I. Hashimoto: “Present Status And Future Needs: The View from Japanese Industry,” Proceedings of the Fourth International Conference on Chemical Process Control, South Padre Island, Texas, February 17-22, 1-28 (1991)

28

1 Introduction

2.1 Introduction

29

Chapter 2 The Distillation Process — An Industrial Example

2.1 Introduction A distillation column is not just any mass-produced article such as a toaster or a washing-machine. Each distillation column is a unique process unit, specially designed for the separation of a particular substance mixture. Nevertheless, the thermodynamic principles and basic dynamics are always the same. Therefore it is possible to demonstrate ideas for the controller design by the example of one column without extensive loss of generality. First in this chapter, the design and operating data of the industrial distillation column are outlined, followed by a brief description of the composition dynamics. The further two sections outline the control objectives, the control structures, and the use of tray temperatures as controlled outputs. The literature references terminate the chapter.

2.2 Column design and operation The distillation column described in this thesis is an industrial binary distillation column. A synopsis of the most important data for this distil-

30

2 The Distillation Process — An Industrial Example

lation column is given in Table 2.1. The distillation column (Fig. 2.1) is equipped with 50 sieve trays, a total condenser, and a steam-heated reboiler. The subcooled feed F enters the column on tray 20 (counted from the top) and for the greater part consists of a mixture of two substances. Because of the small fraction of impurities, these are neglected and the distillation column is considered to be a binary distillation column. The desired product compositions are 0.99 mol/mol (low boiling component) for the top product D and 0.015 mol/mol for the bottom product B. As these product purities are relatively high, this distillation column can be classified as a “high purity distillation column.” Table 2.1: Steady-state data Column data No. of trays

50

Column diameter (m)

0.8

Feed tray

20

Murphree tray efficiency

≈0.4

Relative volatility α

1.61

Operating data Top composition xD (mol/mol)

0.99

Bottom composition xB (mol/mol)

0.015

Feed composition xF (mol/mol)

0.7-0.9

Feed flow rate F (mol/min)

20-46

Top pressure (mbar)

60

Nominal operating point Feed composition (mol/mol)

0.8

Feed flow rate (mol/min)

33

Reflux L0 (mol/min)

65

Boilup V51 (mol/min)

104

2.2 Column design and operation

31

Vacuum Condenser

VT

Top product (Distillate) 1 2 3 4

D, xD L0

Reflux accumulator

Reflux

Feed

20

F, xF

47 48 49 50

Boilup V51 Reboiler B, xB Bottom product

Figure 2.1: The industrial distillation column

32

2 The Distillation Process — An Industrial Example

Feed disturbances The distillation column is connected in series following two other distillation columns, which operate in parallel. The bottom product streams of these two columns are buffered by a tank and fed into the column considered here. The level of the buffer tank is measured periodically (typical period: 2 hours) and the feed of the column is set to keep the tank level within specified bounds. Therefore, the feed flow is varied not continuously but stepwise. In contrast to that, the variations of the feed composition are always smooth. Even a shutdown of one of the other two columns cannot cause a sudden increase of the buffer tank’s composition. Top pressure control The boiling points of the entering substances are high at standard atmospheric pressure. Because of a thermal decomposition of the light component at higher temperatures, the column is operated under vacuum. Correspondingly, the cooling water flow rate for the condenser is kept constant and the top pressure is controlled by a vacuum pump. Top level control The reflux accumulator level is controlled by overflow. Hence the top product flow rate D is not available as a manipulated variable for a composition control system.

2.3 Steady-state behavior Let us assume a composition control scheme with integrating behavior, e.g., one PI controller which controls the top composition by manipulating the reflux and one which controls the bottom composition by manipulating the boilup. Then, in steady-state, the product compositions are kept perfectly at their set-points, and an S-shaped composition profile is developed within the distillation column. Figure 2.2 shows the simulated composition profiles for different feed flow rates and compo-

2.3 Steady-state behavior

33

0

10

20 Tray TrayNo. No.

xF = 0.7 mol/mol xF = 0.8 mol/mol xF = 0.9 mol/mol

30

F = 20 mol/min F = 33 mol/min F = 46 mol/min

40

50 0.0

0.2

0.4 0.6 Composition Composition (mol/mol)

0.8

1.0

Figure 2.2: Simulated composition profiles for the industrial distillation column

34

2 The Distillation Process — An Industrial Example

sitions. While these steady-state profiles are nearly independent of the feed flow rate, they depend essentially on the feed composition. This has a high significance for a controller design: If we want to keep the product compositions close to their setpoints, we must allow profile variations in the middle of the column. Consequently, we cannot control any composition in the middle of the column. The internal flow rates can be illustrated in a similar manner. Figure 2.3 shows the simulated liquid and vapor flow rates for the nominal operating point. As previously mentioned, the reflux as well as the feed are subcooled, i.e. they enter the column at a temperature below the boiling point. A fraction of the vapor flow is condensed at the trays where these two streams are fed into the distillation column. The two 0

10

Liquid flow Vapor flow

Tray No.

20 Figure 2.3: Simulated vapor and liquid flow rates at nominal operating point

30

40

50 60 80 100 120 Flow rate (mol/min)

2.4 Composition dynamics

35

discontinuities of the vapor flow profile at trays 1/2 and 20/21 result therefrom. The reason for the slopes of the two profiles within the stripping and rectifying section of the column is the different heat of evaporation of the two substances.

2.4 Composition dynamics The composition dynamics within a distillation column is effectively described by movements and shape alterations of the composition profile. In order to illustrate this, let us control the reboiler level of the distillation column by the bottom product flow rate B, and let us keep the reboiler heat duty constant. The simulated step responses of the composition profile to a 5% increase and a 5% decrease of the reflux flow rate are shown by Figure 2.4. An increase of the reflux (Fig. 2.4 a) raises the fraction of the light component in the column bottom. Consequently, the composition profile of the light component moves towards the column bottom, degrading the bottom product composition from 1.5% to more than 30% impurity. The opposite effect is observed for a decrease of the reflux flow rate (Fig. 2.4 b): The composition profile moves towards the column top, which improves the bottom product composition and debases the top product composition. These plots illustrate two important properties of the composition dynamics: • Column nonlinearity: The product compositions are a nonlinear function of the reflux, boilup, and the feed condition: A 5%increase of the reflux flow rate improves the top product composition by 0.007 mol/mol, but a 5% decrease degrades it by more than 0.2 mol/mol. • Strong interactions: A change of reflux or boilup alters both product compositions.

36

2 The Distillation Process — An Industrial Example

0

0

10

10

20

20

t

Tray No.

Tray No.

t = 20 h

Δt= 30 min

30

t=0h

30 Δt= 30 min t=0h

40

40

t t = 20 h

50 0.0

50 0.2 0.4 0.6 0.8 1.0 Composition Composition (mol/mol)

a)

0.0

0.2 0.4 0.6 0.8 1.0 Composition Composition (mol/mol)

b)

Figure 2.4: Simulated composition profiles (light component) for a step change of the reflux. Reboiler heat duty, feed flow rate and composition are kept at their nominal values (see Table 2.1) a) L0=1.05*L0,nom

b) L0=0.95*L0,nom

2.5 Control objectives and configurations

37

The interaction between both product compositions and reflux and boilup has a severe consequence for the composition dynamics, usually called • Ill-conditioned behavior. This is best explained by another two examples. If we like to increase both product purities simultaneously, we have to increase reflux and boilup by an exact quantity, for example the reflux by +26.5% and the boilup by +19% (Figure 2.5 a). This keeps the composition profile’s position constant, but it slowly intensifies the S-shape of profile. However a slightly smaller step size for the reflux completely alters the dynamic behavior (Fig. 2.5 b): The purity of the top product decreases, the purity of the bottom product increases, and the dynamic response is much faster. Therefore an exact direction of the input vector [ L, V ] T is required in order to achieve a simultaneous increase of both product purities. Consequently, even a small uncertainty of the input vector [ L, V ] T may lead to undesired results. High condition numbers σ max { G ( jω ) } κ = ---------------------------------σ min { G ( jω ) }

(2.1)

of the plant model G indicate such a behavior.

2.5 Control objectives and configurations The control of distillation columns has three objectives [2.2]: • Control of the material balance (inventory control) • Product quality control • Satisfaction of constraints The first objective includes the control of the vapor holdup (top pressure), the reflux accumulator level, and the reboiler level. Generally, these control objectives are easily achieved by simple PI controllers.

38

2 The Distillation Process — An Industrial Example

0

0

10

10

Δt = 4h

t Δt = 10h

20

t

Tray No.

Tray No.

20

30

t = 40 h

30

t=0h

t t = 40 h t=0h

40 Composition

Composition

40

0

Time (h)

0

40

50

Time (h)

40

50

0.0

0.2 0.4 0.6 0.8 1.0 Composition Composition (mol/mol)

a)

0.0

0.2 0.4 0.6 0.8 1.0 Composition Composition (mol/mol)

b)

Figure 2.5: Simulated composition profiles (light component) for a step change of the reflux and the reboiler heat duty. The feed is kept at nominal condition (see Table 2.1). a) L0=1.265*L0,nom V51=1.19*V51,nom

b) L0=1.260*L0,nom V51=1.19*V51,nom

2.5 Control objectives and configurations

39

The second objective is the most important objective. It is strongly related to the economic and ecological optimal operation of a distillation column. Tight control of both product qualities minimizes the energy consumption and the amount of products being off the specifications. It is not a simple task to keep both product compositions close to their setpoints, especially in the presence of disturbances such as variations of feed flow rate and feed composition. Tight composition control requires sophisticated control schemes. Their design is the main topic of this thesis. Reflux, boilup, and pressure drop are allowed to vary within a predefined range. Any operation of a distillation column outside of this range may cause insufficient separation or even damage of the column. Each control system must handle such constraints to enable safe operation. This topic is well discussed by Buckley et al. [2.2] and Shinskey [2.4]. 2.5.1 The 5×5 control problem A simple distillation column, such as the industrial example discussed here, presents a control problem with the five control objectives • • • • •

Top composition Bottom composition Reflux accumulator level Reboiler level Top pressure

and the five manipulated variables • • • • •

Reflux Boilup (indirectly controlled by reboiler duty) Top product flow rate Bottom product flow rate Cooling water flow rate (or vapor flow rate to vacuum)

This problem is often called the 5×5 control problem. As mentioned above, the top pressure is controlled by a vacuum pump and the reflux

40

2 The Distillation Process — An Industrial Example

Controlled outputs

Manipulated inputs

3×3 control problem Top product xD

Reflux L

Bottom product xB

Boilup V (Reboiler duty Q)

Reboiler level MB

Bottom product flow rate B

Condenser level MD

Top product flow rate D

Top pressure p

Overhead vapor VT (Cooling water flow rate, vacuum pump) 5×5 control problem

Figure 2.6: The distillation control problems

accumulator level by overflow. Thus the 5x5 control problem is reduced to a 3×3 control problem. These relations are illustrated in Figure 2.6. 2.5.2 Control design steps In principle, the design of a MIMO controller for the 5×5 or in this case the 3×3 control problem does not cause any particular difficulties. However, the failure of just one actuator or sensor disables all control loops. Due to the high sensitivity of MIMO controllers to sensor or actuator failure, the inventory control and the composition control usually are independently designed, thus improving the robustness of the control system and simplifying the controller design. The corresponding design approach consists of three steps [2.5]:

2.5 Control objectives and configurations

41

1. Choosing the control configuration In a first step the two manipulated variables for the composition control are to be chosen. This choice names the control configuration. For example, if the top composition xD is controlled by reflux L and the bottom composition xB is controlled by boilup V, the control configuration is called L,V control configuration. After the choice of the manipulated variables for composition control, the remaining three manipulated variables are available for level and pressure control. The choice of the control configuration is often based on configuration selection methods such as Relative Gain Array (RGA), Niederlinski Index, or Singular Value Decomposition (SVD). The application of these indices may lead to very different results (see [2.1], [2.6]), and the reliability seems to be low. One reason for the limited reliability may be the neglect of inventory control: Yang et. al. [2.9] point to the substantial influence of inventory control on the composition control dynamics. Most indices for control configuration selection are based on steadystate gains. Consequently, perfect inventory control is assumed and dynamic effects due to the interaction of inventory and composition control are neglected. The most common control configuration in the chemical industry is the L,V configuration [2.7]. This control structure is rather independent of inventory control dynamics [2.9] and has shown good results within an experimental comparison of different control structures [2.8]. 2. Inventory control design In general, tight inventory control can be achieved with three simple PI controllers. Some distillation columns show an inverse response of the reboiler level to an increase of boilup. In this case, tight level control with boilup as manipulated variable may be difficult. 3. Composition control design A 2×2 controller for composition control is to be designed as a third step of the design. This step is discussed in chapters 5-8.

42

2 The Distillation Process — An Industrial Example

2.6 Tray temperatures as controlled outputs On-line composition analyzers are frequently used to determine product compositions. However, their investment and maintenance costs are prohibitive for distillation columns below a certain size. Provided that substances with a boiling point difference of at least 10 °C are separated and that the product purity specifications are not extremely stringent, pressure-compensated temperatures may substitute composition measurements ([2.2], [2.4]). 2.6.1 Pressure-compensated temperatures For binary mixtures a definite correlation exists between boiling temperature, pressure, and composition T = f ( p, x )

(2.1)

This correlation is illustrated in Figure 2.7 for the two components entering the industrial distillation column. A substitution of the composition measurements by temperature measurements requires a compensation for the effect of pressure variations. If the pressure variations are small, the temperature measurement can be compensated by a linear function. The nominal pressure and composition are denoted by the index N. ∂T T Compensated = T + -----∂p

( p – pN )

(2.2)

N

In case of larger pressure variations, a second-order term has to be supplied:

T Compensated

∂T = T + -----∂p

2

1 ∂ ( p – p N ) + --- 2 T 2∂p N

( p – pN ) 2 N

(2.3)

2.6 Tray temperatures as controlled outputs

Figure 2.7: Boiling points of the two-component-mixture

43

44

2 The Distillation Process — An Industrial Example

Estimation of tray composition It is possible to infer the tray composition directly. By regression of {x, T, p} data, the coefficients of a simple polynomial expression can be calculated. An example is given by x = θ 1 + θ 2 ( T + T Corr ) + θ 3 p + θ 4 p 2

(2.4)

Such an equation in terms of the absolute temperature and pressure is simpler to implement in a distributed control system than an equation in terms of deviations from reference values x = θ1 + θ2 ( T – TN ) + θ3 ( p – pN ) + θ4 ( p – pN )2

(2.5)

One problem of the tray composition estimate is a potential bias of the temperature measurements. Practical experience has shown that a bias of up to 2 °C is to be expected due to heat transport phenomena. In (2.4) the bias is corrected by the parameter T Corr . In practice, however, this correction is difficult to estimate. In principle, it would be possible to include cross terms such as θTp in the regression model. However, errors in the absolute temperature may lead to incorrect numerical values of these cross terms. Therefore, in the regression model, cross terms should be avoided. Pressure compensation as well as the estimation of tray composition are easily implemented in a process control system. Without a pressure compensation, it is impossible to use tray temperatures in a vacuum column as controlled variables and expensive composition analyzers are necessary. For temperature measurements close to the column top, a linear compensation is usually sufficient. For trays close to the column bottom, we have to expect higher pressure variations, and a compensation with a second-order polynomial is recommended. 2.6.2 Temperature measurement placement The sensitivity of the tray temperatures near the ends of the column to changes of the product compositions is very small. To make the temper-

2.7 References

45

ature measurement sensitive enough, it has to be located at some distance from the column ends. Figure 2.2 shows simulated steady-state composition profiles for the industrial distillation column. These profiles illustrate the fact that the effect of a change of operating conditions increases with growing distance from the column ends. On the other hand, a deterioration of the correlation between product composition and tray temperature results from an increasing distance from the column ends. A compromise between correlation with product composition and sensitivity must thus be found. Kister discusses the most important rules and tools in [2.3]. In the case of the industrial distillation column, the temperatures on trays 10 and 44 are chosen as controlled outputs. Additionally, the temperature on tray 24 is measured. Since tray 24 is close to the feed tray, it is expected to be sensitive to any change of feed composition and, dynamically, to the feed flow rate.

2.7 References [2.1] Ariburnu, D., C. Ozge, and T. Gurkan: “Selection of the Best Control Configuration for an Industrial Distillation Column,” Preprints of 3rd IFAC Symposium on Dynamics and Control of Chemical Reactors, Distillation Columns and Batch Processes, April 26-29, 1992, College Park, MD, 387-392 (1992) [2.2] Buckley, P. S., W. L. Luyben, and J. P. Shunta: Design of Distillation Column Control Systems, Instrument Society of America, Research Triangle Park, NC (1985) [2.3] Kister, H. Z., Distillation Operation, McGraw-Hill, New York (1990) [2.4] Shinskey, F. G., Distillation control for productivity and energy conservation, 2. ed., McGraw Hill, New York (1984) [2.5] Skogestad, S., and M. Morari: “Control Configuration Selection for Distillation Control,” AIChE J., 33, 10, 1620-1635 (1987)

46

2 The Distillation Process — An Industrial Example

[2.6] Skogestad, S., P. Lundström, and E. W. Jacobsen: “Selecting the Best Distillation Control Configuration,” AIChE J., 36, 5, 753764 (1990) [2.7] Skogestad, S.: “Dynamics and Control of Distillation Columns ⎯ A Critical Survey,” 3rd IFAC Symposium on Dynamics and Control of Chemical Reactors, Distillation Columns and Batch Processes, April 26-29, 1992, College Park, MD, 1-25 (1992) [2.8] Waller, K. V., D. H. Finnerman, P. M. Sandelin, K. E. Häggblom, and S. E. Gustafsson, “An Experimental Comparison of Four Control Structures for Two-Point Control of Distillation,” Ind. Eng. Chem. Res., 27, 624-630 (1988) [2.9] Yang, D. R., D. E. Seborg, and D. A. Mellichamp: “The Influence of Inventory Control Dynamics on Distillation Composition Control,” Preprints of the 12th World Congress of the International Federation of Automatic Control, Sydney, 18-23 July 1993, 1, 7176 (1993)

3.1 Introduction

47

Chapter 3 A Rigorous Dynamic Model of Distillation Columns

3.1 Introduction The rigorous dynamic process simulation has become an accepted and widespread tool in process and even more so in controller design [3.11]. Increasing competition and environmental protection provisions require an optimization of process and control structures, which can be obtained only by a substantial knowledge of process dynamics. At the same time, dynamic experiments on a running plant are less and less desired. Rigorous dynamic modelling and simulation can replace such expensive and time-consuming measurements. This has special significance for high-purity distillation columns. Due to their long time constants and varying feed flow rates and feed compositions, reproducible operating conditions are difficult to guarantee. Therefore, new controllers are usually tested thoroughly by dynamic simulation for the full operating range of the distillation column. The rigorous models of distillation columns used for that purpose match the reality to a large extent [3.17]. In this chapter, a rigorous dynamic model for distillation columns is discussed. This model is used in all nonlinear dynamic simulations

48

3 A Rigorous Dynamic Model of Distillation Columns

within this thesis. In a special section, the numerical treatment of the resulting system of algebraic-differential equations is outlined. The modelling and control fields use very different notations. Therefore the notation used within this chapter is explained in section 3.11.

3.2 Conventions Figure 3.1 shows a schematic representation of a distillation column equipped with nt trays. The column top (condenser and reflux accumulator) is denoted by the index 0, the trays with the indices 1, 2, … nt, and the column bottom (including the reboiler) with the index nt+1. To simplify the formal mathematical description the reflux stream R is designated as liquid flow (L0). The feed of the industrial distillation column, as described in Chapter 2, is in liquid phase and subcooled. The top pressure is controlled by a vacuum pump and the condenser is operated with a constant cooling water flow rate. Flash calculations for the feed stream as well as dynamic models for the top pressure of the column are therefore not considered here. For other applications, the model presented is easily extended with appropriate model equations.

3.3 The objective of modelling The control or process engineer is interested in the dynamic behavior of various important process variables (e.g., tray temperatures, product compositions) as a function of the time-varying column inputs. The objective of a dynamic model is an approximation of the real process input/output behavior by a system of differential and algebraic equations. These model equations are based on material and energy balances as well as on thermodynamic and fluid dynamic correlations.

3.4 Simplifying assumptions Within a distillation column many different physical phenomena occur. Although it would be possible to include models for the fluid streams on

3.4 Simplifying assumptions

49

Q0

0 …………… …………

R

1

Condenser

Reflux accumulator

D

(=L0)

2 3 Sl,4

4

Fj

Vj j Lj

Fnt-2

nt-2

Sv,nt-1

nt-1 nt ………… ……

Vnt+1

Qnt+1 B

nt+1

Figure 3.1: Distillation column

Reboiler

50

3 A Rigorous Dynamic Model of Distillation Columns

the trays, for the dead time caused by the transport time of vapor flow from one tray to the next one above, or for the heat exchange with the environment, the resulting model would be of very high order. As mentioned earlier, the aim of modelling the distillation column dynamics is a sufficient description of the real macroscopic behavior. This means that we are interested primarily in the dynamics of tray compositions, temperatures, and pressures etc. rather than in the fluid streams on the trays. Experience shows that no substantial improvement can be achieved with models including effects with more microscopic characteristics. Hence the following assumptions are usually introduced in order to achieve a compromise between model accuracy and order ([3.3], [3.13], [3.17]): • The holdup of the vapor phase is negligible compared to the holdup of the liquid phase. • Liquid phase and vapor phase are each well mixed on all trays, i.e., the composition of the liquid and of the vapor phase are independent of the position on the tray. • The residence time of the liquid in the downcomer is neglected. • The variation of the liquid enthalpy on a tray can be neglected on all trays. (This assumption is not applicable to the evaporator.) In the literature so far, uniform liquid flows and constant holdups for all trays have often been assumed (equimolar overflow). This assumption is problematic because it implies a neglect of flow dynamics. Essential dynamic effects may remain unmodelled, e.g., a non-minimum phase behavior (inverse response) of the reboiler level and the tray compositions in the lower section of the column to an increase in reboiler heat supply.

3.5 Balance equations

51

3.5 Balance equations 3.5.1 Material balances The differential equations describing the dynamics of the holdup for each component on a tray are derived from a material balance for each component. The balance border is the single tray with its ingoing and outgoing streams (Figure 3.2). Lj-2 xk,j-2 . . . . .. . . .. . . .. . . . . . . . . . .. .

j-1 Vj yk,j

Fj xF,k,j

Vj-1 yk,j-1

Lj xk,j . . .. . . . . . .. . . . . . . . . . . . . .

j

Lj-1 xk,j-1 . . .. . .. . . .. . . . . .. . . .. . . .

Vj+1 yk,j+1

Sl,j-1 xk,j-1

Sl,j xk,j Sv,j+1 yk,j+1

j+1 Vj+2 yk,j+2

Lj+1 xk,j+1

Figure 3.2: Balance border for the material balances

Material balance for component k on tray j (k=1, …, nc; j=1, …, nt) dn k, j d ( n j x k, j ) - = F j x F, k, j + L j – 1 x k, j – 1 – ( L j + S l, j )x k, j -------------- = ----------------------dt dt + ( V j + 1 – S v, j + 1 )y k, j + 1 – V j y k, j

(3.1)

52

3 A Rigorous Dynamic Model of Distillation Columns

In the same way, the balance equations for the column top and the column bottom are formulated: Material balance for component k in condenser (k=1, …, nc) d ( n 0 x k, 0 ) dn k, 0 - = ( V 1 – S v, 1 )y k, 1 – ( L 0 + D )x k, 0 -------------- = ------------------------dt dt

(3.2)

Material balance for component k in the evaporator (k=1, …, nc) Usually the liquid phases in the column bottom and the reboiler are mixed either by natural convection or by a pump. Assuming perfect mixing we obtain dn k, nt + 1 d ( n nt + 1 x k, nt + 1 ) ------------------------- = --------------------------------------------dt dt

(3.3)

= L nt x k, nt – Bx k, nt + 1 – V nt + 1 y k, nt + 1 The total holdup on tray j equals the sum of the holdups of the individual substances: nc

nj =

∑ n k, j

(3.4)

k=1

3.5.2 Energy balance equations The vapor flows within a distillation column are calculated by an energy balance. The balance border is the same to the border in Figure 3.2, which was used for the material balance equations. Energy balance for tray j: d ------ ( n j h' j ) = F j h' F, j + L j – 1 h' j – 1 + ( V j + 1 – S v, j + 1 )h'' j + 1 dt – ( S l, j + L j )h' j – V j h'' j For the left-hand side of this equation the following holds

(3.5)

3.5 Balance equations

53

dn dh' d ------ ( n j h' j ) = h' j ---------j + n j ----------j dt dt dt

(3.6)

If in (3.6) we substitute the expression for the differential term dn j ⁄ dt according to dn ---------j = F j + L j – 1 + V j + 1 – S v, j + 1 – S l, j – L j – V j dt

(3.7)

the following energy balance equation holds dh' j n j ---------- = [ F j ( h' F, j – h' j ) + L j – 1 ( h' j – 1 – h' j ) dt

(3.8)

+ ( V j + 1 – S v, j + 1 ) ( h'' j + 1 – h' j ) – V j ( h'' j – h' j )] Usually, the assumption n j ( dh' j ⁄ dt ) ≈ 0 is permissible, except for cases with large temperature variations on the trays, a large heat of mixing, or a large tray holdup. With this assumption we can rewrite equation (3.8) as an algebraic expression for the vapor flow rate V j 1 V j = -------------------- [ F j ( h' F, j – h' j ) + L j – 1 ( h' j – 1 – h' j ) h'' j – h' j

(3.9)

+ ( V j + 1 – S v, j + 1 ) ( h'' j + 1 – h' j )] A similar balance equation is formed for the evaporator. Because of the large inventory, the derivative n j ( dh' j ⁄ dt ) cannot be neglected. Since an increase in vapor flow causes an increase in bottom pressure and consequently an increase of boiling temperature in the evaporator, the vapor flow follows any change in reboiler heat supply with a time lag. Hajdu et al. [3.9] present a model for this vapor flow lag. We can imagine that an energy stream Q supplied to the evaporator is subdivided into two fractions: One part causes an evaporation of liquid, the other increases the bottom temperature. Written as a differential equation we obtain the energy balance equation

54

3 A Rigorous Dynamic Model of Distillation Columns

ΔQ = ΔH v, nt + 1 ΔV nt + 1 dΔT nt + 1 + n nt + 1 ν' nt + 1 ρ' nt + 1 c' p, nt + 1 -----------------------dt

(3.10)

To achieve a first-order differential equation in ΔVnt+1, the differential term dΔTnt+1/dt has to be substituted by a differential term in ΔVnt+1. The increase of the pressure drop due to a changing vapor flow rate (assuming a constant total holdup on the tray) can be estimated with ∂Δp j Δ ( Δp j ) = ⎛ ----------------⎞ ⎝ ∂V j + 1⎠

ΔV j + 1

(3.11)

nj

Hence the pressure change in the evaporator can be approximated for a distillation column with nt trays by ∂Δp j Δ ( p nt + 1 ) = nt ⎛ ----------------⎞ ⎝ ∂V j + 1⎠

ΔV nt + 1

(3.12)

nj

The increase in boiling point temperature caused by the increase in bottom pressure can be calculated according to ∂T nt + 1 ΔT nt + 1 = ⎛ -------------------⎞ ⎝ ∂p nt + 1 ⎠

Δ ( p nt + 1 )

(3.13)

x k, nt + 1

Substituting (3.13) in equation (3.10), the following differential equation is obtained: ΔQ – ΔH v, nt + 1 ΔV nt + 1 = ∂T nt + 1 ∂Δp j dΔV nt + 1 n nt + 1 ν' nt + 1 ρ' nt + 1 c' p, nt + 1 ⎛ -------------------⎞ nt ⎛ ----------------⎞ -----------------------⎝ ∂p nt + 1 ⎠ ⎝ ∂V j + 1⎠ dt

(3.14)

Therefore, the vapor flow lag at an increase in reboiler heat supply can be described by the first-order lag

3.6 Fluid dynamics

55

dQ lag 1 = ----------- ( Q – Q lag ) --------------T lag dt

(3.15)

with the time constant

T lag

∂T nt + 1 ∂Δp j n nt + 1 ν' nt + 1 ρ' nt + 1 c' p, nt + 1 ⎛ -------------------⎞ nt ⎛ ----------------⎞ ⎝ ∂p nt + 1 ⎠ ⎝ ∂V j + 1⎠ = -----------------------------------------------------------------------------------------------------------------------------ΔH v, nt + 1

(3.16)

If we substitute the total bottom holdup balance equation in the energy balance equation dn nt + 1 h' nt + 1 -------------------- = L nt h' nt + Q lag – Bh' nt + 1 – V nt + 1 h'' nt + 1 dt

(3.17)

the following equation holds: Energy balance for the evaporator L nt ( h' nt – h' nt + 1 ) + Q lag V nt + 1 = ----------------------------------------------------------------h'' nt + 1 – h' nt + 1

(3.18)

The parameters ( ∂T nt + 1 ) ⁄ ( ∂p nt + 1 ) and ( ∂Δp j ) ⁄ ( ∂V j + 1 ) can be evaluated numerically or analytically from the appropriate equations (see sections 3.6.2 and 3.7.2)

3.6 Fluid dynamics In the previous sections, the equations describing composition and total holdup dynamics, as well as the vapor flow rates have been derived. Here the calculation of the liquid flow rates and of the pressure drop is discussed. 3.6.1 Liquid flow rates The volumetric liquid flow rate over the weir on tray j can be calculated according to the Francis weir formula ([3.16], [3.10]):

56

3 A Rigorous Dynamic Model of Distillation Columns

2 3/2 L V, j = μ 2g --- b W h LW ,j 3

(3.19)

For sharp-edged weirs μ ≈ 0.64 holds. Perfect mixing on the trays, including the liquid in the downcomers, is assumed. Nevertheless, if we calculate the effective liquid head h LW, j above the weir edge, we have to take the liquid phase fraction εj and the liquid volume in the downcomer into account (Figure 3.3). The liquid level in the downpipe is the sum of the liquid head on the tray and of the hydrostatic level due to the pressure drop according to pj – pj – 1 Hydrostatic liquid level in downcomer = ----------------------ρ' j g

(3.20)

The liquid head h L, j of the pure liquid on a tray (without a vapor phase fraction) is equal to the total liquid volume on the tray n j ν' j minus the

pj-1

pj h LW, j --------------εj h L, j ---------εj

hW L V, j

pj – pj – 1 ----------------------ρ' j g

AA

hL,j AB

Figure 3.3: Liquid levels on a tray

3.6 Fluid dynamics

57

liquid volume in the downcomer due to pressure A B ( p j – p j – 1 ) ⁄ ( ρ' j g ) , both divided by the total area A A + A B :

h L, j

pj – pj – 1 n j ν' j – ----------------------- A B ρ' j g = -----------------------------------------------AA + AB

drop

(3.21)

For the application of the Francis weir formula, we have to evaluate the liquid level of the pure liquid (liquid without vapor phase fraction). For that purpose, first the height of the two-phase layer is to be evaluated and second the liquid phase fraction εj must be taken into account. The effective liquid level becomes

h L, W, j

pj – pj – 1 n j ν' j – ----------------------- A B ρ' j g h L, j = ⎛ ---------– h W⎞ ε j = ------------------------------------------------ – εj hW ⎝ ε ⎠ + A A j A B

(3.22)

Substituting (3.22) into the Francis weir formula (3.19), we obtain the volumetric liquid flow rate of the two-phase mixture. The flow rate from tray j in molal units is calculated by: 3/2 pj – pj – 1 ⎛ n ν' – ---------------------⎞ AB ⎜ j j ⎟ ρ' j g 2 μ 2g --- b W ⎜ ------------------------------------------------- – ε j h W⎟ 3 AA + AB ⎜ ⎟ ⎝ ⎠ L j = -----------------------------------------------------------------------------------------------------------ν' j

(3.23)

In many industrial distillation columns, calming zones exist in front of the weir. For this case, ε j = 1 holds at the weir edge. Otherwise, we have to estimate the liquid phase fraction on the trays. The Stichlmair correlation is well suited for that purpose [3.18]. 3.6.2 Pressure drop A vapor flow through a tray in a distillation column suffers a pressure drop. Its amount depends on the vapor flow rate, the tray holdup, and

58

3 A Rigorous Dynamic Model of Distillation Columns

the geometry of the tray. Usually, the pressure drop is assumed to consist of three different parts ([3.7], [3.12]): • Dry pressure drop occurring at the flow through the tray without liquid ( Δp tr, j ) • Hydrostatic pressure drop due to liquid head and liquid density (Δp L, j ) • Pressure drop by bubble-forming due to surface tension of liquid (Δp σ, i ) The pressure drop by bubble-forming usually is insignificant and can be neglected. Dry pressure drop With sufficient accuracy, the dry pressure drop can be approximated by the following well-known expression: Δp tr, j

ρ'' j ⎛ V V, j + 1⎞ 2 = ξ ( Re ) ------- ------------------2 ⎝ A0 ⎠

(3.24)

The orifice coefficient ξ ( Re ) either can be evaluated by measurement on comparable trays, or it can be estimated with experimentally verified correlations. During the simulations, the following correlation for sieve trays is used [3.19]: A 2 s –1.27 ρ'' j ⎛ V V, j + 1⎞ 2 ------- ------------------Δp tr, j = ⎛ 1 – -------0 ⎞ + 0.21 ⎛ ------⎞ ⎝ ⎝ d 0⎠ 2 ⎝ A0 ⎠ A A⎠

(3.25)

Hydrostatic pressure drop The hydrostatic pressure drop results from the liquid head and the liquid density. We have to take the liquid volume in the downcomer into account (see 3.6.1).

3.7 Phase equilibrium

59

Δp L, j

pj – pj – 1 n j ν' j – ----------------------- A B ρ' j g = ------------------------------------------------- ρ' j g AA + AB

(3.26)

The total pressure drop consists of the sum of the two parts dry pressure drop and hydrostatic pressure drop: Δp j = p j + 1 – p j = Δp tr, j + Δp L, j

(3.27)

3.7 Phase equilibrium All equations we have discussed in the previous sections are explicitly or implicitly interrelated with the vapor phase composition. In this section, the most important correlations concerning the vapor phase compositions and boiling points are presented. 3.7.1 Vapor phase composition The liquid on each tray and in the evaporator is at boiling-point. Phase equilibrium thus can be assumed. At moderate pressures up to some few bar, the concentration of a substance in the vapor flow leaving tray j can be obtained according to y kEquilibrium ,j

γ k, j p k0, j = -------------------- x k, j = K k, j x k, j pj

(3.28)

If the substance mixture exhibits ideal behavior, the activity coefficient γ becomes one, and the vapor phase compositions are equal to the ratios of the partial pressures of the substances and the absolute pressure on the tray. 0

The vapor pressures of the pure substances p k can be calculated with a high level of accuracy by the Antoine equation (3.29). The parameters A, B, and C are listed in many tables of substance properties (e.g., [3.5]). 0 B log [ p k ] = A – ------------T+C

(3.29)

60

3 A Rigorous Dynamic Model of Distillation Columns

The calculation of the liquid phase activity coefficients γ k, j can be effected by one of the well known correlations (Wilson, NRTL, UNIQUAC etc.). Murphree tray efficiency In a distillation column only little contact time exists on each tray for the mass transfer between liquid and vapor phase. Therefore no perfect phase equilibrium can be achieved, and the tray efficiency will deviate from the unit value. This effect can be modelled by the Murphree tray efficiency for the vapor phase. y k, j – y k, j + 1 η = ------------------------------------------------– y k, j + 1 y kEquilibrium ,j

(3.30)

3.7.2 Boiling points The vapor phase composition according to (3.28) is a function of the tray temperature Tj. At boiling point, the sum of the vapor phase mole fractions calculated becomes one. Hence for a tray j, the following boiling point equation holds: nc ,j ∑ ykEquilibrium k=1

nc

=

∑ Kk, j ( Tj, pj, xk, j )xk, j

= 1

(3.31)

k=1

The Murphree tray efficiency is not considered for the boiling point calculation, because it relates to the mass transfer between vapor and liquid phase rather than to the equilibrium composition.

3.8 Volumetric properties The fluid dynamic models discussed are interrelated with the molar volumes of the vapor phase and of the liquid phase, and with the corresponding densities. Their calculation is the subject of this section.

3.8 Volumetric properties

61

3.8.1 PVT relations The molar volumes of the liquid phase ν' j or the vapor phase ν'' j depends on the pressure pj, the temperature Tj, and the actual compositions xk,j or yk,j. A great number of different equations of state has been developed to describe this behavior. They are extremely well documented ([3.5], [3.6]), and a discussion of their properties is not repeated here. The PVT behavior is described here by the Soave-Redlich-Kwong equation (SRK equation, [3.15], [3.6]) with the Peneloux correction. This correction improves the estimate of the molar volumes of the liquid phase, which is overrated by 10-15% using the SRK equation. If measurement data of the PVT behavior of the pure substances exist and their mixing behavior is nearly ideal, a different possibility has shown good results for the liquid phase: We can correlate the molar volumes measured with the temperature by a polynomial regression. The molar volume ν' j of the substance mixture can be approximated as a weighted sum of the individual molar volumes: nc

ν' j =

∑ xk, j ν'k, j

(3.32)

k=1

3.8.2 Density The densities of liquid and vapor phase can be computed from the molar volume, the molar mass, and the mole fractions. nc

∑ x k, j M k k=1 Liquid phase density: ρ' j = ----------------------------ν' j

(3.33)

nc

∑ y k, j M k k=1 Vapor phase density: ρ'' j = ----------------------------ν'' j

(3.34)

62

3 A Rigorous Dynamic Model of Distillation Columns

3.9 Enthalpies The quantity not discussed so far is the enthalpy of a substance mixture in liquid or vapor phase. The enthalpy of a real fluid is estimated by the real sum of an ideal part and the value of a departure function Δh T, P describing the deviation of the enthalpy from the enthalpy of the ideal gas state: T 0

h = h +



id

real

c p dT + Δh T, P

(3.35)

T0

The ideal part can be calculated by summing the ideal parts for each component: T

∫ T0

⎛ T ⎞ ⎜ c pid dT = ∑ x k ∫ c pid, k dT⎟ ⎜ ⎟ ⎝ ⎠ k=1 T nc

(3.36)

0

id

The ideal heat capacities c p are often approximated by a third-order polynomial for each component: c pid, k = A k + B k T + C k T 2 + D k T 3

(3.37)

The parameters for equation (3.37) are listed in many tables of substance properties, or they can be estimated with very high accuracy by Joback’s method ([3.15], p. 154-156). real

The real part Δh T, P describes the departure of a mixture from the ideal behavior. It can be evaluated using one of the well-known equations of state, e.g., the SRK equation ([3.15], [3.6]). If measurement data for the heat capacities and for the heat of vaporization are available, a simple solution is possible in a manner similar to that mentioned in section 3.8.1:

3.10 Numerical solution

63

T

j ⎛ ⎞ ⎜ Liquid phase enthalpy: h' j = ∑ x k, j ∫ c p, j, k dT⎟ + h 0 ⎜ ⎟ ⎝ ⎠ k=1 T

nc

(3.38)

0

T

⎛ j ⎞ ⎜ Vapor phase enthalpy: h'' j = ∑ y k, j ∫ c p, j, k dT + ΔH v, j, k⎟ + h 0 (3.39) ⎜ ⎟ ⎝ ⎠ k=1 T nc

0

3.10 Numerical solution The complete rigorous dynamic model for distillation columns, as introduced above, consists of a system of differential and algebraic equations (DAE). The complexity of the model is illustrated by Figure 3.4. It illustrates the interconnection of the model equations for three adjoining trays. The solution of the differential equations obviously depends on the solution of the algebraic equation system. Therefore an efficient numerical integration using standard integration methods is not possible. This requires special adapted integration algorithms, as outlined in section 3.10.4. 3.10.1 The dependent variables and the equation system As a first step for the numerical treatment, we have to decide which variables should form the vector of the dependent variables. This vector of dependent variables must at any time completely describe the state of a distillation column and should be of minimum size to avoid excessive computation times. The vapor phase composition is an illustrative example for the complete description of the distillation’s state: If we know the tray composition, the tray temperature, and the tray pressure, then the vapor phase composition in equilibrium is easily calculated by an explicit algebraic equation. Consequently, it is not necessary to insert the vapor phase composition into the vector of the dependent variables.

nk,j+1

xk,j+1

nj

xk,j

nk,j+1 ∑nk,j+1 nj+1

nk,j ∑nk,j

xk,j-1 nk,j-1 ∑nk,j-1 nj-1

×



k, j + 1 ----------------dt -



nk,j

nk,j-1

× BP

BP

Δpj-1

pj+2

Murphree

Murphree

Δpj+1

Eq. of state

Tj+1

Δpj

Eq. of state

y*k,j+1

pj+1

Tj

y*k,j

pj

Eq. of state

Tj-1

Murphree States

Flow dynamics

States

Flow dynamics

States

Flow dynamics

yk,j+1

yk,j

yk,j-1

×

dn

nc

k, j ---------dt

dn



k, j – 1 ---------------dt -

dn

BP

y*k,j-1

× Vj+1

× Vj+2

h’’j+1-h’j

Vj

h’’j-h’j-1

Vj-1

×

Figure 3.4: Model structure for three adjoining trays h’’j+2-h’j+1

Lj+1

h’’j+1-h’j+1

h’j-h’j+1

Lj

h’’j-h’j

h’j-1-h’j

Lj-1

h’’j-1-h’j-1

h’j-2-h’j-1

× /

× × /

Δpj-2

× /

pj-1

64 3 A Rigorous Dynamic Model of Distillation Columns

×

×

3.10 Numerical solution

65

As one vector which satisfies the requirements of a complete description and of minimum order, the following vector is proposed (as a modification of the vector proposed by Holland & Liapis [3.10]): y = [QKond, D, n1,0, …, nnc,0, T0, p0, L0, {Vj, n1,j, …, nnc,j, (Sl,j), (Sv,j), Tj, pj, Lj}j=1, 2, …, nt Qlag, Q, Vnt+1, n1,nt+1, …, nnc,nt+1, B, Tnt+1, pnt+1, States of the control system]T ( ): Value is inserted only if it physically exists The Jacobian matrix of the equation system (as described below) corresponding to these dependent variables has a numerically advantageous block diagonal dominant structure. For the calculation of these dependent variables y, the following equations are to be solved Differential equations nc material balance equations (3.1) Algebraic equations 1 equation for vapor flow rate (3.9) 1 equation for liquid flow rate (3.23) 1 equation for boiling point (3.31) 1 equation for pressure drop (3.27) Total: nc +4 equations per tray and in addition the equations for the evaporator, the condenser, and the control system. Considering industrial distillation columns which are often equipped with more than 50 trays, the resulting algebraic differential equations amount to several hundred equations. The model for the industrial binary distillation equipped with 50 trays gives an impression of these numbers: It consists of a system of 107 differential and 210 algebraic equations.

66

3 A Rigorous Dynamic Model of Distillation Columns

3.10.2 Formal representation of the DAE We can formally represent the entire dynamic model by the semiexplicit equation system dn ------- = f ( t, n ( t ), z ( t ) ) dt

n ( t0 ) = n0

(3.40a)

0 = g ( t, n ( t ), z ( t ) )

z ( t0 ) = z0

(3.40b)

The vector n consists of all tray holdups (for all components), while the vector z contains all other dependent variables. A different but equivalent formal representation is the implicit form: F ( t, y ( t ), y' ( t ) ) = 0

y ( t0 ) = y0

(3.41)

Here the vector y contains all the dependent variables. A simulation of the dynamic behavior requires a simultaneous solution of the whole equation system. 3.10.3 The index The index of a set of differential-algebraic equations (DAE) characterizes the integration problem. The higher the index, the more difficult is a solution of the DAE. The differential index is the most common definition: The differential index m of the system

F ( t, y ( t ), y' ( t ) ) = 0

imal number m such that the system of F ( t, y ( t ), y' ( t ) ) = 0

is the minand of the

analytical differentiations d m ( F ( t, y ( t ), y' ( t ) ) ) d ( F ( t, y ( t ), y' ( t ) ) ) - = 0 ------------------------------------------------ = 0, …, --------------------------------------------------m dt dt can be transformed by algebraic manipulations into an explicit ordinary differential system [3.8]. Consequently, a system of ordinary differential equation has an index of m=0.

3.10 Numerical solution

67

3.10.4 Solution methods and software The first general method for the numerical solution of semi-explicit DAE with index 1 was proposed by C. W Gear in 1971 [3.4] and was soon extended to the solution of implicit index 1 problems. The method is based on a special class of the linear multistep methods entitled the backward differentiation formulas, which are standard algorithms for the integration of stiff systems. The most important convergence results may be found in [3.1]. In theory, it is possible to solve problems of higher indices with the backward differentiation formulas. However, the necessary software is not available as yet. The apparently very frequently used integrators DASSL and LSODI are based on Gear’s method. These methods are distinguished for their effectiveness in solving continuous problems. However, the computational effort grows significantly for systems with discontinuities arising, for example, during the simulation of the response to several feed flow or feed composition step changes. For such cases, the one-step methods find more and more interest [3.11]. The one-step methods are extensions of the well-known Runge-Kutta, Rosenbrock, or extrapolation methods. An extensive discussion of the properties of these methods is found in [3.8]. However, the development of the integrators (RADAU5, LIMEX) is in an early stage, and no implementations are found in any of the widespread Fortran libraries. For the simulation studies the DASSL integrator, as implemented in the NAG Fortran Library is used with good success. The differential-algebraic equations (DAE) are solved in an implicit manner according to (3.41). The calculation sequence During the integration, the right-hand sides of the differential and algebraic equations repeatedly have to be evaluated for a given vector y of the dependent variables and for a given time t. The algebraic equations, and often the differential equations as well are solved in an implicit manner. The equation errors, which have to be supplied to the integration, are the difference between the right-hand sides of the equations

68

3 A Rigorous Dynamic Model of Distillation Columns

(that means the calculated vapor flow rates, liquid flow rates, etc.) and the corresponding values within the vector y. A correct calculation sequence evaluating these terms is stringent: If, for example, we calculate the right-hand side of the boiling point equation (3.31), we must know the vapor phase composition. Therefore we first have to calculate the vapor phase composition and subsequently the error of the boiling point equation. If this basic idea is applied to the whole model, the calculation sequence illustrated in Figure 3.5 results. The vector y, which is supplied by the integration routine (Step a), contains the component holdups in liquid phase, the tray temperatures, tray pressures, as well as liquid and vapor flow rates. With the data supplied, all vapor phase compositions in equilibrium can be calculated (Step b). In a next step (Step c), using the distribution coefficients K k, j obtained in the previous step, the errors of the boiling point equations are calculated. Since all the feed data are known, its enthalpy, molar volumes and densities are computed in step d. The vapor phase compositions deviate from the equilibrium compositions. Applying the Murphree tray efficiency, the effective vapor phase compositions are computed (Step e). Since for the computation of the effective vapor phase composition for a tray the effective vapor phase composition of the next lower tray must be known, the computation starts at the column bottom, assuming η = 1 for the reboiler. Now all necessary data are known to calculate the enthapies, molar volumes, and densities for all trays, the condenser and the evaporator (Step f). In step g, the energy balance equations for the trays are applied, and the differences between the resulting vapor flow rates and the flow rates supplied by the integration routine are calculated. In step h, the same is done for the evaporator. Similarly to the error of the vapor flow rates, the errors of the liquid flow rates and the tray pressures are computed in steps i and j. Using the flow rates and compositions calculated in the previous steps rather than the data supplied by the integration routine, the differential terms (left-hand sides) of the equations describing the vapor flow lag (Step k) and the holdup of the substances in liquid phase (Step l) are calculated. In a last step (Step m), all differential terms, the errors

3.10 Numerical solution

69

a)

Vector of dependent variables y

b)

Vapor phase composition for evaporator and trays (Equation (3.28))

Error for boiling point at evaporator and trays c)

1 – ∑ K k, j x k, j k

Calculation of the thermodynamic states d)

h', ν', ν'', ρ', ρ''

for the feed

Murphree tray efficiency for trays nt, nt-1, …, 1 e)

y k, j = η ( y kEquilibrium – y k, j + 1 ) + y k, j + 1 ,j

Calculation of the thermodynamic states f)

h' j, h'' j, ν' j, ν'' j, ρ' j, ρ'' j

for condenser, all trays, and evaporator

Error for vapor flows g)

F j ( h' F, j – h' j ) + L j – 1 ( h' j – 1 – h' j ) + ( V j + 1 – S v, j + 1 ) ( h'' j + 1 – h' j ) ------------------------------------------------------------------------------------------------------------------------------------------------------------------------- – V j h'' j – h' j

Figure 3.5: Calculation sequence Explanation: see text

70

3 A Rigorous Dynamic Model of Distillation Columns

Error for vapor flow leaving evaporator h)

L nt ( h' nt – h' nt + 1 ) + Q lag ------------------------------------------------------------------ – V nt + 1 h'' nt + 1 – h' nt + 1

Error for liquid streams i)

3/2 pj – pj – 1 ⎛ n ν' – ---------------------⎞ - AB j j ρ' j g ⎜ ⎟ 2 1 μ 2g --- b W ⎜ ------------------------------------------------- – ε j h W⎟ ------ – L j 3 ν' AA + AB ⎜ ⎟ j ⎝ ⎠

Error for pressure drop j)

p j + 1 – p j – Δp j

(Equation (3.27))

k)

Differential equation for vapor flow lag (Equation (3.15))

l)

Differential equations for holdup of substances (Equations (3.1), (3.2) and (3.3))

m)

Vector of differentials and errors

Figure 3.5 continued

3.11 Notation

71

between supplied and calculated flow rates and pressures, and the errors of the boiling point equations are combined in one vector and supplied back to the integration routine.

3.11 Notation A0

[m2]

Hole area in tray

AA

[m2]

Tray area without downcomer area

AB

[m2]

Downcomer area

bW

[m]

Length of weir

c pid

[J/mol·K]

Ideal gas heat capacity

cp,l

[J/kg·K]

Liquid heat capacity

d0

[m]

Diameter of holes of sieve tray

Fj

[mol/s]

Feed flow rate to tray j

h

[J/mol]

Molar enthalpy

h' j

[J/mol]

Molar enthalpy of liquid phase

h'' j

[J/mol]

Molar enthalpy of vapor phase

hL

[m]

Liquid level above upper edge of weir

hW

[m]

Weir height

ΔHv,k,j [J/mol]

Heat of evaporation of component k on tray j

ΔHv,j

[J/mol]

Heat of evaporation of liquid on tray j

Kk,j

[mol/mol]

Distribution coefficient for comp. k on tray j

Lj

[mol/s]

Liquid flow leaving tray j

LV,j

[m3/s]

Volumetric flow from tray j

Mk

[g/mol]

Molar mass of component k

nt

[–]

Number of trays in column

72

3 A Rigorous Dynamic Model of Distillation Columns

nj

[mol]

Total holdup on tray j

nk,j

[mol]

Holdup of substance k on tray j

nc

[–]

Number of components

pj

[N/m2]

Pressure on tray j

Δpj

[N/m2]

Pressure drop over tray j

p k0

[N/m2]

Steam pressure of pure component k

P

[N/m2]

Pressure

Q

[J/s]

Heat supply to evaporator

Qlag

[J/s]

“active” heat supply

s

[m]

Thickness of sieve tray

Sl,j

[mol/s]

Side product flow rate from tray j, liquid phase

Sv,j

[mol/s]

Side product flow rate from tray j, vapor phase

t

[s]

Time

T

[K]

Temperature

Tj

[K]

Temperature on tray j

Vj

[mol/s]

Vapor stream from tray j

VV,j

[m3/s]

Volumetric vapor stream from tray j

xk,j

[mol/mol]

Liquid phase mole fraction of component k on tray j

xF,k,j

[mol/mol]

Mole fraction of component k in feed to tray j

yk,j

[mol/mol]

Vapor phase mole fraction of component k above tray j

γk

[–]

Liquid phase activity coefficient of component k

3.11 Notation

73

εj

[m3/m3]

Liquid phase fraction on tray j

η

[mol/mol]

Murphree tray efficiency for vapor phase

ν

[m3/mol]

Molar volume

ν' j

[m3/mol]

Molar volume of liquid phase on tray j

ν'' j

[m3/mol]

Molar volume of vapor phase on tray j

ξ

[–]

Orifice coefficient

ρ' j

[kg/m3]

Liquid density on tray j

ρ'' j

[kg/m3]

Vapor density on tray j

74

3 A Rigorous Dynamic Model of Distillation Columns

3.12 References [3.1] Brenan, K. E., S. L. Campbell, and L. R. Petzold, Numerical solution of initial-value problems in differential-algebraic equations, North-Holland, New York (1989) [3.2] Byrne, G. D., P. R. Ponzi, Differential-Algebraic Systems, Their Application and Solution, Comp. Chem. Eng., 12, 5, 377-382 (1988) [3.3] Gani, R., C. A. Ruiz, and I. T. Cameron: “A Generalized Model for Distillation Columns,” Comp. Chem. Eng., 10, 3, 181-198 (1986) [3.4] Gear, C. W.: “Simultaneous Numerical Solution of DifferentialAlgebraic Equations,” IEEE Trans. on Circuit Theory, CT-18, 1, 89-95 (1971) [3.5] Gmehling, J. and U. Onken: “Vapor-Liquid Equilibrium Data Collection,” 1, Part 1, XI-XXII, DECHEMA, Frankfurt (1977) [3.6] Gmehling, J. and B. Kolbe: Thermodynamik, Georg Thieme Verlag, Stuttgart (1988) [3.7] Grassmann, P. and F. Widmer, Einführung in die thermische Verfahrenstechnik, 2nd ed., de Gruyter, Berlin (1974) [3.8] Hairer, E. and G. Wanner: Solving Ordinary Differential Equations II ⎯ Stiff and Differential-Algebraic Problems, Springer Verlag, Berlin (1991) [3.9] Hajdu, H., A. Borus, and P. Földes: “Vapor Flow Lag in Distillation Columns,” Chem. Eng. Sc., 33, 1-8 (1978) [3.10] Holland, C. D. and A. I. Liapis, Computer Methods for Solving Dynamic Separation Problems, Chapter 8, McGraw-Hill, New York (1983)

3.12 References

75

[3.11] Marquardt, W.: “Dynamic Process Simulation — Recent Progress and Future Challenges,” Fourth International Conference on Chemical Process Control, South Padre Island, Texas (1991) [3.12] McCabe, W. L., J. C. Smith, and P. Harriott: Unit Operations of Chemical Engineering, 4th ed., McGraw-Hill, New York (1985) [3.13] Najim, K. (Editor): Process Modeling and Control in Chemical Engineering, Marcel Dekker, New York (1989), Chapter III, 145211, S. Domenech, L. Pibouleau, “Distillation” [3.14] Petzold, L.: “Differential/Algebraic Equations are not ODE,” SIAM J. Sci. Stat. Comput, 3, 3, 367-384 (1982) [3.15] Reid, R. C., J. M. Prausnitz, and B. E. Poling: The Properties of Gases and Liquids, 4th ed., McGraw-Hill, New York (1988) [3.16] Retzbach, B.: “Mathematische Modelle von Destillationskolonnen zur Synthese von Regelungskonzepten,” Fortschritt-Berichte VDI, Reihe 8: Mess-, Steuerungs- und Regelungstechnik, Nr. 126, VDI Verlag (1986) [3.17] Rovaglio, M., E. Ranzi, G. Biardi, and T. Faravelli: “Rigorous Dynamics and Control of Continuous Distillation Systems — Simulation and Experimental Results,” Comp. Chem. Eng., 14, 8, 871887 (1990) [3.18] Stichlmair, J.: Grundlagen der Dimensionierung des Gas/Flüssigkeit-Kontaktapparates Bodenkolonne, Verlag Chemie, Weinheim (1978) [3.19] Weiss, S. et. al.: Verfahrenstechnische Berechnungsmethoden, Teil 2: “Thermisches Trennen”, VCH Verlagsgesellschaft, Weinheim (1986)

76

3 A Rigorous Dynamic Model of Distillation Columns

4.1 Introduction

77

Chapter 4 Linear Models

4.1 Introduction Robust controllers are designed on the basis of linear process models. Therefore the elaboration of linear dynamic models for the distillation column is a central part of control system synthesis. These models should describe the dynamic behavior of the process within a wide frequency range. They can be obtained in two ways: • System identification • Linearization of a nonlinear model It is a big advantage of the system identification that it avoids a complicated and expensive nonlinear model. Nevertheless, this approach has some severe drawbacks, for example: • The time-constants of the composition dynamics are large. A recording of input/output data for the real plant is very timeconsuming. • Due to the high sensitivity of distillation columns to changes of the internal flow rates, even for small magnitudes of the input variation (e.g., 5% of the steady-state value) the response may far exceed the linear region.

78

4 Linear Models

• Each experiment causes undesired disturbances of the product qualities. • It is practically impossible to obtain models for the entire operating range of the distillation column These disadvantages and some other fundamental problems of the identification itself (see Jacobsen et al. [4.5]) lead to a strong recommendation of the second method (Skogestad, [4.12]) that means the linearization of nonlinear column models. Two linear models are evaluated within this chapter, which are based on the linearization of different nonlinear column models. The first linear model is obtained by an analytical linearization of a simplified nonlinear model neglecting flow dynamics. The other model is obtained by a numerical linearization of the rigorous model presented in Chapter 3. In further sections the accuracy of these linear models and the role of the flow dynamics are discussed. Different mathematically order reduction methods are compared at the end of this chapter. The notation is listed in section 4.9 on page 101, the literature references are collected in section 4.10.

4.2 How to linearize the rigorous model? 4.2.1 The state, input, and output vectors The complete rigorous dynamic model as discussed in Chapter 3 consists of a high-order system of coupled differential and algebraic equations (DAE). A linearization of this large system would be possible in principle. However, the resulting linear state-space model would be of the same order as the DAE. Such a high order causes high computation times for a controller design or even for a model reduction. Consequently, a compromise between model order and accuracy must be sought. This means we have to decide, which dependent variables are very important for the composition dynamics and should be included in the state vector x of the linear model. Most important, of course, are the tray compositions themselves. Because flow dynamics have a high influ-

4.2 How to linearize the rigorous model?

79

ence on the composition dynamics in the high-frequency range, the tray holdup is a candidate as well. Assuming a perfect level control of the reboiler and the reflux accumulator, it is not necessary to include their holdup. The corresponding candidates for the state vector of the linear model are

dx 0 dx 1

dx 0 x =

dx 1 … dx 51

x = n

or

… dx 51

(4.1)

dn 1 … dn 50

The dynamics of the distillation column are stimulated by the manipulated variables (reflux L0 and boilup V51) and the several disturbance sources. Most important disturbances are variations of the feed composition xF and the feed flow rate F. Other disturbances such as variations of the reflux temperature or the feed temperature have significantly less influence and can be neglected for the composition control design. Hence we define the input vector according to dx F d = u

dF

(4.2)

dL 0 dV 51

The output vector y follows directly from the temperature measurements. It represents the deviations of the pressure compensated temperatures on tray 10, 44, and 24:

80

4 Linear Models

dT P 10 y = dT P 44

(4.3)

dT P 24 4.2.2 Handling of the algebraic equation system The algebraic equation system of the rigorous column model defines dependent variables such as tray pressures, vapor flows, and liquid flows, which are not included in the state vector x of the linear model. Nevertheless the algebraic equation system represents algebraic constraints for the composition and holdup dynamics. These equations can be handled in two ways: • elimination by idealizing assumptions, or • numerical solution during linearization The first method allows an analytical linearization of the resulting model. This has the advantage that merely one steady-state data set must be supplied, which can be calculated by one of the common flowsheeting programs (e.g., PROCESSTM, ASPEN PLUSTM). In contrast to that, the second method requires a numerical linearization of the rigorous model, which is discussed in section 4.4.

4.3 Linearization of a simplified nonlinear model 4.3.1 The simplified model In this section we will derive a simplified nonlinear column model neglecting the holdup and thus the flow dynamics. For that purpose idealizing assumptions are formulated which allow to dispense with all flow dynamics and with most of the energy balance equations. However, the subcooling of reflux and feed have a significant influence on the internal flow rates (see section 2.3) and are explicitly taken into account.

4.3 Linearization of a simplified nonlinear model

81

Idealizing assumptions The algebraic constraints of the rigorous model and the holdup dynamics can be eliminated by the following idealizing assumptions: • constant pressure drop • constant and equal enthalpies on all trays • constant total holdup on all trays (equimolar overflow) Of course, all these assumptions do not agree with the real conditions. The first assumption means a neglect of the correlation between tray pressures, holdups, and boilup rates. The second assumption implies uniform vapor flows within the stripping section and within the rectifying section of the column. The assumption of a constant tray holdup contains a neglect of flow dynamics. The error in the high frequency range introduced by that is discussed in section 4.5. It has to be emphasized here that these assumptions concern only the simplified nonlinear model as a basis for an analytical linearization. The steady-state operating points must be calculated using a model, which includes the energy balance equations, as well as the flow dynamic models. The composition dynamics For a column separating a binary mixture, it is sufficient to formulate the material balance equation for the light component of the substance mixture. If we assume constant total holdup, the following material balance equation holds for the tray j: dn ---------j = L j – 1 + V j + 1 + F j – L j – V j = 0 dt

(4.4)

Similar balances are obtained for the reboiler and the condenser. Substituting these balance equations in the material balances (3.1)-(3.3), we obtain the following differential equations describing the composition dynamics:

82

4 Linear Models

Condenser dx 0 1 --------- = ------ V 1 ( y 1 – x 0 ) n0 dt

(4.5)

Trays (Feed is liquid phase, j = 1…50) dx 1 --------j = ----- [ L j – 1 ( x j – 1 – x j ) + V j + 1 ( y j + 1 – x j ) nj dt

(4.6)

– V j ( y j – x j ) + F j ( x F, j – x j ) ] Evaporator dx 51 1 ----------- = -------- [ L 50 ( x 50 – x 51 ) – V 51 ( y 51 – x 51 ) ] n 51 dt

(4.7)

Effect of subcooled reflux and feed Feed and reflux of the distillation column are subcooled. A portion of the vapor flow is condensed at the trays where these two streams enter the column. The effect of an additional condensation of the vapor stream caused by increasing the flow rates of these streams must be considered to avoid large model errors. The two energy balance equations for the reflux tray (tray 1) and for the feed tray (tray 20) become part of the nonlinear model: 0 = L 0 h' 0 – L 1 h' 1 + V 2 h'' 2 – V 1 h'' 1

(4.8)

0 = L 19 h' 19 – L 20 h' 20 + V 21 h'' 21 – V 20 h'' 20 + Fh' F

(4.9)

Tray temperatures as model outputs The model outputs are the deviations of the pressure compensated tray temperatures. These temperatures are correlated with the tray composition by the boiling point equation

4.3 Linearization of a simplified nonlinear model

83

2

∑ ( yk* ) – 1

= 0

(4.10)

k=1

Consequently the boiling point equation for tray 10, 44, and 24 are part of the simplified column model. The vapor phase composition The differential equations for the composition dynamics are in terms of the vapor phase composition yj, too. Usually, the tray efficiency is smaller than one and the vapor phase compositions deviate from the equilibrium compositions. As described in Chapter 3, this can be modelled by the Murphree tray efficiency η y j = ( 1 – η )y j + 1 + ηy * j

(4.11)

Primarily the vapor phase composition in equilibrium y * j is a function of the liquid phase composition x j . For example, if we assume a constant relative volatility α j on tray j, the vapor phase composition in equilibrium is calculated by αj xj y * j = --------------------------------1 + ( α j – 1 )x j

(4.12)

To simplify the analytical linearization, it is convenient to substitute the vapor phase composition y j by a correlation exclusive in terms of the equilibrium compositions y * . Then the calculation of the derivatives ∂ ( … ) ⁄ ∂x for each these terms cause no particular problem. Such an equation in terms of the equilibrium compositions y * is derived by subsequent substitutions of the vapor phase compositions in (4.11) from the evaporator up to actual tray. For example (η=1 assumed for the reboiler) y 50 = ( 1 – η )y * 51 + ηy * 50

(4.13)

84

4 Linear Models

y 49 = ( 1 – η )y 50 + ηy * 49 = ( 1 – η ) 2 y * 51 + η ( 1 – η )y * 50 + ηy * 49 y 48 = ( 1 – η )y 49 + ηy * 48 = ( 1 – η ) 3 y * 51 + η ( 1 – η ) 2 y * 50 + η ( 1 – η )y * 49 + ηy * 48

(4.14)

(4.15)

For a binary distillation column with nt trays the following generalized formula (with vapor/liquid equilibrium (η=1) in the evaporator) is obtained: y j = ( 1 – η ) nt + 1 – j y * nt + 1 nt + 1 – j

+



η ( 1 – η ) nt + 1 – j – n y * nt + 1 – n

(4.16)

n=1

This equation demonstrates the strong influence on the composition of the trays below the actual one, presuming the tray efficiency η is substantially smaller than one. 4.3.2 Analytical linearization Let us formally represent the simplified nonlinear model as the vector functions f and g: dx˜ ( t ) ˜ ( t ), d˜ ( t ) ] --------------- = f [ x˜ ( t ), u dt

(4.17a)

y˜ ( t ) = g [ x˜ ( t ) ]

(4.17b)

Then the matrices A, B, and C of the linear state space model dx = Ax + B d -----dt u

(4.18a)

y = Cx

(4.18b)

4.3 Linearization of a simplified nonlinear model

85

are evaluated as partial derivatives of the vector functions f and g at a steady-state operating point (OP). The OP is calculated either solving the equation system of the complete rigorous model for steady state, or using a steady-state flowsheeting program. The following relations hold for the partial derivatives: ∂f ∂f ∂f 0 -------0- -------0- … ---------∂x 0 ∂x 1 ∂x 51 ∂f A = -----∂x˜

= OP

∂f ∂f -------1- -------1- … … ∂x 0 ∂x 1

(4.19)

… … … ∂f 51 ∂f 51 --------- … … ---------∂x 0 ∂x 51

OP

∂f ∂f ∂f ∂f 0 --------0- ------0- --------0- ----------∂x F ∂F ∂L 0 ∂V 51 ∂f B = ---------˜ ∂ d u˜

∂f ∂f ∂f ∂f 1 --------1- ------1- --------1- ----------= ∂x F ∂F ∂L 0 ∂V 51 OP

… ∂f 51 --------∂x F

… ∂f 51 --------∂F

… ∂f 51 --------∂L 0

… ∂f 51 ----------∂V 51

(4.20)

OP

∂g 1 ∂g 1 ∂g 1 -------- -------- … ---------∂x 0 ∂x 1 ∂x 51 g C = ∂----∂x˜

OP

∂g ∂g ∂g 2 = -------2- -------2- … ---------∂x 0 ∂x 1 ∂x 51 ∂g 3 ∂g 3 ∂g 3 -------- -------- … ---------∂x 0 ∂x 1 ∂x 51

(4.21)

OP

86

4 Linear Models

Combining the idealizing assumptions, we can conclude that a deviation of reflux L0 or boilup V51 causes the same deviation of the liquid flow and vapor flow rates within the whole column: dL j = dL 0 dV j = dV 51

(4.22)

The resulting coefficients of the matrices A, B, and C are listed in the appendix of this chapter (page 97). Although important interactions in the column model are suppressed by the idealizing assumptions it turns out that this model coincides within acceptable bounds with the rigorous nonlinear model. This aspect is discussed in detail in section 4.5 below.

4.4 Linearization of the rigorous model In this section the linearization of the rigorous dynamic column model is discussed. This model includes the dynamics of the tray holdups and thus flow dynamics. Hence it describes the high-frequency dynamics much better than the simplified nonlinear model (section 4.3.1). 4.4.1 Model modifications The desired outputs of the linear model are the deviations of pressure compensated tray temperatures, which are functions of the tray compositions, but not of a component’s holdups. Therefore the material balance equations (3.1)-(3.3) are replaced by the following, equivalent differential equations: Condenser dx 0 1 --------- = ------ V 1 ( y 1 – x 0 ) n0 dt

(4.23)

4.4 Linearization of the rigorous model

87

Tray holdups: dn ---------j = L j – 1 – L j + V j + 1 – V j + F j dt

(4.24)

Tray composition (Feed is liquid phase, j = 1…50) dx j 1 -------- = ----- [ L j – 1 ( x j – 1 – x j ) + V j + 1 ( y j + 1 – x j ) nj dt

(4.25)

– V j ( y j – x j ) + F j ( x F, j – x j ) ] Evaporator dx 51 1 ----------- = -------- [ L 50 ( x 50 – x 51 ) – V 51 ( y 51 – x 51 ) ] n 51 dt

(4.26)

The second modification of the rigorous model follows from the assumption of a perfect level control in the reflux accumulator and in the evaporator. This assumption is justified by the fact that the level control loops for the condenser and the evaporator can be tuned much faster than the composition control loops. With perfect level control and therefore constant holdup for the reflux accumulator and the evaporator, the top and bottom product streams are calculated according to D = V1 – L0

(4.27)

B = L 50 – V 51

(4.28)

and

The whole equation system consists of the differential equations describing the composition dynamics for the light component (4.5)-(4.7) and the holdup dynamics (4.24), and the algebraic equations for the pressure drop (3.27), the vapor flows (3.9), liquid flows (3.23), and the boiling points (3.31). These algebraic equation form an equation system whose solution is a constraint for the differential equations. All these equations can be comprised by the following vector equations:

88

4 Linear Models

dx˜ ( t ) = f [ x˜ ( t ), n ˜ ( t ), u ˜, u ˜ ( t ), z˜ ( t ), v˜ ( x˜ , n ˜ , z˜ ) ] --------------dt

(4.29a)

˜ (t) dn ˜, u ˜ ( t ), u ˜ , z˜ ) ] ˜ ( t ), z˜ ( t ), v˜ ( x˜ , n --------------- = l [ x˜ ( t ), n dt

(4.29b)

y ( t ) = g [ x˜ ( t ) ]

(4.29c)

˜, u ˜ , z˜ ) ] 0 = k [ v˜ ( x˜ , n

(4.29d)

˜, u ˜ , z˜ ) represents the solution of the algebraic equaThe vector v˜ ( x˜ , n tion system k and consists of the tray pressures, the vapor flow rates, the liquid flow rates, and the boiling points. 4.4.2 Numerical linearization The matrices of the linear state space model including flow dynamics d ------ x = A x + B d dt n n u

(4.30a)

y = Cx

(4.30b)

can be numerically evaluated column by column using a finite difference approach. After solution of the whole equation system for a steady-state operating point (OP), each composition or tray holdup can be varied by a small increment, the algebraic equation system k can be solved, and each column of the state dynamic matrix A can be calculated according to

a i=1…102, j + 1 =

⎛ dx ⎞ dx – ------⎜ ------⎟ ⎝ dt x + Δx dt OP⎠ j j ----------------------------------------------------Δx j ⎛ dn ⎞ – dn ------⎜ ------⎟ ⎝ dt x + Δx dt OP⎠ j j -----------------------------------------------------Δx j

for j = 0…51

(4.31a)

4.5 Comparison of the linear models

a i=1…102, j + 52 =

⎛ dx ⎞ – dx ------⎜ ------⎟ ⎝ dt n + Δn dt OP⎠ j j -----------------------------------------------------Δn j ⎛ dn ⎞ dn – ------⎜ ------⎟ ⎝ dt n + Δn dt OP⎠ j j ------------------------------------------------------Δn j

89

for j = 1…50

(4.31b)

In a similar manner the input matrix B can be evaluated.

4.5 Comparison of the linear models 4.5.1 Open loop simulations Most important for the control design is a good representation of the dynamic behavior in the mid-frequency range. The steady-state behavior as well as the high-frequency behavior are less important. Some idea of a linear model’s quality is obtained by a simple qualitative comparison of the various models. Nevertheless, a definitive judgement requires a comparison of control designs based on the different models. A simple method to compare the two linear models with the complete rigorous model is the simulation of step responses to the model inputs (reflux L0, boilup V51, feed composition xF, and feed flow rate F). These are shown by the Figures 4.1- 4.4. During the nonlinear simulations, the bottom level was controlled by the bottom product flow rate B. Except for the denoted input, all other column inputs are kept constant at their steady-state values. The changes of flow rates and feed composition are very small to maintain the column close to the steady-state and to avoid large nonlinearities. The coincidence of the step responses with the rigorous nonlinear model is acceptable for both linear models. However, the linear model obtained by a linearization of the rigorous model is distinguished by a somewhat better representation of the low-frequency gains.

90

4 Linear Models

Tray 44 0

-0.05

-0.2

Temp. dev. (K)

Temp. dev. (K)

Tray 10 0

-0.1 -0.15 -0.2 0

2

4 6 Time (h)

8

10

-0.4 -0.6 -0.8 0

2

4 6 Time (h)

8

10

Temp. dev. (K)

Tray 24 0

0

-0.1

-0.2

-0.2

-0.4

-0.3

-0.6

-0.4

-0.8

Nonlinear model Anal. linearized model

-0.5 0

2

4 6 Time (h)

8

-1 0

10

Num. linearized model

2

4

6

8

10

Figure 4.1: Step response to a 0.3 mol/min (0.46%) increase in reflux Tray 10

Tray 44 0.5

0.15

Temp. dev. (K)

Temp. dev. (K)

0.2

0.1 0.05 0 0

0.4 0.3 0.2 0.1

2

4 6 Time (h)

8

10

0 0

2

4 6 Time (h)

8

10

Temp. dev. (K)

Tray 24 0.5

0

0.4

-0.2

0.3

-0.4

0.2

-0.6

0.1

-0.8

Nonlinear model Anal. linearized model

0 0

2

4 6 Time (h)

8

10

-1 0

Num. linearized model

2

4

6

8

Figure 4.2: Step response to a 0.3 mol/min (0.29%) increase in boilup

10

4.5 Comparison of the linear models

91

Tray 10

Tray 44 0.1

-0.02

Temp. dev. (K)

Temp. dev. (K)

0

-0.04 -0.06 -0.08 -0.1 0

2

4 6 Time (h)

8

0 -0.1 -0.2 -0.3 0

10

2

4 6 Time (h)

8

10

Temp. dev. (K)

Tray 24 0

0

-0.05

-0.2

-0.1

-0.4

-0.15

-0.6

-0.2

-0.8

Nonlinear model Anal. linearized model

-0.25 0

2

4 6 Time (h)

8

-1 0

10

Num. linearized model

2

4

6

8

10

Figure 4.3: Step response to a 0.005 mol/mol increase in feed composition Tray 10

Tray 44 0 -0.1

Temp. dev. (K)

Temp. dev. (K)

-0.02 -0.04 -0.06 -0.08 -0.1 -0.12 0

-0.2 -0.3 -0.4 -0.5

2

4 6 Time (h)

8

10

-0.6 0

2

4 6 Time (h)

8

10

Tray 24

Temp. dev. (K)

0

0 -0.2

-0.1

Nonlinear model

-0.4 Anal. linearized model

-0.2 -0.6 -0.3 -0.4 0

Num. linearized model

-0.8 2

4 6 Time (h)

8

10

-1 0

2

4

6

8

10

Figure 4.4: Step response to a 0.3 mol/min (0.91%) increase in feed flow rate

92

4 Linear Models

Surprising is the high coincidence for the analytically linearized model as well. Apparently the influence of the algebraic constraints on the composition dynamics is substantially smaller than the interactions within the composition dynamics themselves. 4.5.2 Singular values An important difference between the two linear models is the high frequency behavior due to unmodelled and modelled flow dynamics, respectively. Best suited for a comparison are the singular values of the transfer functions G d → y ( jω ) and G u → y ( jω ) , shown in Figure 4.5.

Magnitude

10

10

10

Disturbance inputs

3

0

-3

-6

10 -4 10

10

-3

-2

10

10

-1

0

10

10

1

Frequency (rad/min)

Magnitude

10

10

10

Control inputs

3

0

-3

-6

10 -4 10

10

-3

-2

10

10

-1

Frequency (rad/min) Figure 4.5: Singular values of the linear models Upper plots: G d → y ( jω ) , lower plots: G u → y ( jω ) Solid lines: Analytically linearized model Dashed lines: Numerically linearized model

0

10

10

1

4.5 Comparison of the linear models

93

Both models show the typical course of the singular values for a high purity distillation column. In the low frequency range, the maximum and minimum singular values of the transfer functions G u → y ( jω ) are very different and the condition numbers σ max { G u → y ( jω ) } κ ( jω ) = ---------------------------------------------σ min { G u → y ( jω ) }

(4.32)

are high. With increasing frequency, the maximum and minimum singular values approach and the corresponding condition number decreases, but never falls below 20. This is illustrated in more detail by Figure 4.6 for the transfer functions from the control inputs to the pressure compensated temperatures on tray 10 and 44. The large condition numbers indicate a high sensitivity of the column outputs to the direction of the control inputs u. Consequently, the performance of a control system can be very sensitive to uncertainty at the control inputs. Significant for the numerically linearized model is the double as big condition number in the low-frequency range and the completely different course in the high-frequency range. These differences of the high-frequency range between the models can be explained from the structure of the nonlinear models they are calculated from: As mentioned above the analytical model neglects the flow dynamics. Thus the high frequency behavior is determined only by the first-order equations of the composition dynamics. Therefore, the singular values in the higher frequency domain (above 0.1 rad/min) are dominated by a negative slope of one magnitude per decade. The numerically linearized model takes the flow dynamics into account. Thus, considering the reflux as column input, additional lags for the composition dynamics are introduced, and for the minimum singular value a negative slope of several magnitudes per decade for the frequency range above 1 rad/min results therefrom. The effect of the flow dynamics considering the boilup as the column input is different: An increase of the boilup increases the vapor fraction on the tray, which

94

4 Linear Models

Singular values

3

Magnitude

10

0

10

-3

10

-6

10 -4 10

Magnitude

10

10

10

-3

10

-2

10

10

-1

Frequency (rad/min)

10

0

1

10

Condition number

3

2

1

0

10 -4 10

-3

10

-2

10

10

-1

Frequency (rad/min)

10

0

1

10

Figure 4.6: Upper plot: Singular values of the transfer function from the control inputs u to the tray temperature T10 and T44 Lower plot: Condition number κ of the same transfer function Solid lines: Analytically linearized model Dashes lines: Numerically linearized model

causes higher liquid flow rates leaving the trays. Because the composition of the light component is higher in the upper part of the column, more of the light component is transported down, and the expected decrease of the light component’s composition is retarded in the frequency range between 0.2 and 1 rad/min.

4.6 Order reduction The orders of the linear models developed above are 52 for the analytical model and 102 for the numerically evaluated model. During modern

4.6 Order reduction

95

robust control synthesis procedures such as H∞ or μ-synthesis, the order of the model is enlarged by frequency-dependent weights for the model inputs and outputs. Since the computation time for the controller design strongly depends on the model order, order reduction is of utmost importance. Many methods exist to approximate the state-space representation of a linear system with a lower-order state-space approximation [4.13]. Most of the mathematical methods which are available in MATLABTM or MATRIXXTM toolboxes are based on computing the Hankel Norm singular values and subsequent removing of states corresponding to relatively small Hankel Norm singular values. Jacobsen et al. [4.5] compared the following four methods, with a reduction of a column model of 82 states to 2 states. These methods are implemented in one of the MATLABTM toolboxes: I

Balanced Truncated Approximation [4.7] (Robust Control Toolbox [4.2] and μ-Tools [4.1])

II

Balanced Truncated Approximation without balanced minimal realization [4.9] (Schur method, in Robust Control Toolbox)

III Hankel Norm Approximation [4.3] (μ-Tools [4.1]) IV Optimal Hankel Norm Approximation without balancing [4.8] (Robust Control Toolbox) Jacobsen et al. conclude that the methods II and IV gave significantly better models than the other two methods. These results have to be considered carefully: It is not necessary to reduce the column models to such an extremely low order. Models of an order 10-15 are absolutely suitable for control synthesis and show a very good coincidence with the full-order linear model. As an example, a numerically evaluated linear model of order 102 was reduced to an order 10 using each of the four methods mentioned above. All step responses to the different inputs and calculated with the different reduced-order models have shown a perfect coincidence with the full-order linear model. This fact is supported by the singular-value

96

4 Linear Models

3

10

0

Magnitude

10

Full order model

Full order model

-3

Methods I + II

10

Methods I+II

Method III

Method III

Method IV

Method IV

-6

10 -4 10

-3

10

-2

10

10

-1

0

10

1

10

Frequency (rad/min) Figure 4.7: Singular value plots of the transfer functions G u → y ( jω ) of the full order model and the different reduced order models

plots of the models. Figure 4.7 shows that all reduced order models approximate the low and medium frequency behaviors up to 1 rad/min very well. However, in the high frequency range the singular values are best approximated by the models derived with a Balanced Truncated Approximation (Method I or II).

4.7 Summary This chapter presented two methods to obtain linear models for the industrial distillation column. The first model is derived by an analytical linearization of a simplified nonlinear model neglecting flow dynamics and most of the energy balance equations. The second linear

4.8 Appendix: Model coefficients

97

model is obtained by a numerical linearization of the complete rigorous model. Both linear models exhibit an acceptable approximation of the process dynamics. The singular value plots indicate a high coincidence of the linear models in the mid-frequency range, but significant differences in the low and high-frequency range. Comparing step responses with those of the rigorous nonlinear model, a slightly better representation of the column dynamics by the numerically evaluated model is demonstrated. The relatively high order of the linear model (52 for the analytically, and 102 for the numerically linearized models) can be reduced essentially by one of the well known order reduction methods. All tested methods yielded a nearly perfect approximation of the linear model of order 102 up to a frequency of 1 rad/min by a model of order 10.

4.8 Appendix: Model coefficients For all coefficients the following holds: kj =

∂ * y ∂ xj j

(4.33)

Analytically differentiating the equation (4.12), the actual numerical values of k j may be calculated. A-Matrix Condenser (k=1…50) ∂f –V a 1, 1 = -------0- = ---------1n0 ∂x 0

a 1, k + 1

∂f 0 η ( 1 – η ) k – 1 V 1 k k = --------- = --------------------------------------------∂x k n0

∂f 0 ( 1 – η ) 50 V 1 k 51 a 1, 52 = ---------- = ------------------------------------∂x 51 n0

(4.34)

(4.35)

(4.36)

98

4 Linear Models

Trays (j=1…50, k = j…50) ∂f j Lj – 1 = ----------a j + 1, j = -------------∂x j – 1 nj

(4.37)

∂f – ( L j + V j + 1 – V j + V j ηk j + F j ) a j + 1, j + 1 = -------j = -----------------------------------------------------------------------------nj ∂x j

(4.38)

( V j + 1 – ( 1 – η )V j ) ∂f a j + 1, k + 1 = --------j- = η ( 1 – η ) k – j – 1 k k ----------------------------------------------nj ∂x k

(4.39)

( V j + 1 – ( 1 – η )V j ) ∂f j a j + 1, 52 = ---------- = ( 1 – η ) 50 – j k k ----------------------------------------------nj ∂x 51

(4.40)

Evaporator ∂f 51 L 50 a 52, 51 = ---------- = -------∂x 50 n 51

(4.41)

∂f 51 – ( B + V 51 k 51 ) a 52, 52 = ---------- = -----------------------------------n 51 ∂x 51

(4.42)

B-Matrix We have to consider the portion of the vapor flow which is condensed by an increase of the feed flow rate because of the subcooled feed. For the decrease of the vapor flow in the rectifying section of the distillation column h' 20 – h' F dV = – ⎛ --------------------------⎞ dF ⎝ h'' 20 – h' 20⎠

(4.43)

The liquid flow rate in the stripping section of the column is increased by the same amount.

4.8 Appendix: Model coefficients

99

Condenser ∂f b 1, 1 = --------0- = 0 ∂x F

(4.44)

h' 20 – h' F y 1 – x 0 ∂f b 1, 2 = ------0- = – ⎛ --------------------------⎞ ----------------⎝ h'' 20 – h' 20⎠ n 0 ∂F

(4.45)

Rectifying section (j = 1 … 19) ∂f b j + 1, 1 = --------j- = 0 ∂x F

(4.46)

h' 20 – h' F y j + 1 – y j ∂f b j + 1, 2 = ------j = – ⎛ --------------------------⎞ ----------------------⎝ h'' 20 – h' 20⎠ nj ∂F

(4.47)

Feed tray (20) ∂f 20 F 20 = -------b 21, 1 = --------∂x F n 20

b 21, 2

h' 20 – h' F x F – x 20 + -------------------------- ( y 20 – x 20 ) h'' 20 – h' 20 ∂f 20 = --------= --------------------------------------------------------------------------------∂F n 20

(4.48)

(4.49)

Stripping section (j = 21 … 51) ∂f b j + 1, 1 = --------j- = 0 ∂x F

(4.50)

∂f h' 20 – h' F ⎞ x j – 1 – x j ---------------------b j + 1, 2 = ------j = ⎛ 1 + -------------------------nj ∂F ⎝ h'' 20 – h' 20⎠

(4.51)

If the reflux is subcooled, the vapor stream in the column is condensed partially at the first tray. The liquid stream leaving the first tray and all

100

4 Linear Models

trays below is thereby increased by the same amount. From (4.8) we obtain h' 1 – h' 0 dV 1 = – ⎛ ----------------------⎞ dL 0 ⎝ h'' 1 – h' 1⎠

(4.52)

h' 1 – h' 0 ⎞ dL 1 = ⎛ 1 + --------------------- dL ⎝ h'' – h' ⎠ 0

(4.53)

1

1

With (4.52) and (4.53) there follows for the elements of the B matrix: Condenser

b 1, 3

h' 1 – h' 0 --------------------- ( x – y1 ) ∂f 0 h'' 1 – h' 1 0 = ------- = --------------------------------------------∂L n0

∂f y1 – x0 b 1, 4 = ------0- = ---------------∂V n0

(4.54)

(4.55)

Tray 1

b 2, 3

h' 1 – h' 0 x 0 – x 1 + ( y 1 – x 1 ) ---------------------h'' 1 – h' 1 ∂f = ------1- = -------------------------------------------------------------------∂L n1

∂f y2 – y1 b 2, 4 = ------1- = ---------------∂V n1

(4.56)

(4.57)

Trays (j = 2 … 50) xj – 1 – xj ∂f h' 1 – h' 0 ⎞ b j + 1, 3 = ------j = ---------------------- ⎛ 1 + --------------------⎝ nj ∂L h'' 1 – h' 1⎠

(4.58)

∂f yj + 1 – yj b j + 1, 4 = -------j = ---------------------∂V nj

(4.59)

4.9 Notation

101

Evaporator ∂f 51 x 50 – x 51 ⎛ h' 1 – h' 0 ⎞ = ---------------------- 1 + --------------------b 52, 3 = --------n 51 ⎝ ∂L h'' 1 – h' 1⎠

(4.60)

∂f 51 – ( y 51 – x 51 ) b 52, 4 = --------= -----------------------------∂V n 51

(4.61)

C-Matrix The coefficients of the measurement output matrix are numerically evaluated by solving the boiling point equation for small increments in tray composition.

4.9 Notation 4.9.1 Matrices and Vectors A

State dynamic matrix

B

Control input matrix

C

Measurement output matrix

G

Transfer function matrix

G u → y Transfer function matrix from control signals u to output signals y n

Vector of holdup deviations from operating point (OP) n = [ dn 1, dn 2, …, dn 50 ] T

˜ n

Vector of holdups ˜ = [ n , n , …, n ] T n 1 2 50

u

Vector of the manipulated variables (L0, V51) deviations u = [ dL 0, dV 51 ] T

˜ u

Vector of column inputs ˜ = [ L , V ]T u 0 51

102

4 Linear Models

x

Vector of composition deviations from OP x = [ dx 0, dx 1, …, dx 50, dx 51 ] T



Vector of tray compositions x˜ = [ x 0, x 1, …, x 51 ] T

y

Vector of the deviations of the pressure-comp. temperatures y = [ dT P , dT P , dT P ] T 10



44

24

Vector of pressure compensated tray temperatures y˜ = [ T P , T P , T P ] T 10 44 24

d

Vector of the disturbance input deviations from OP z = [ dx F, dF ] T



Vector of disturbance inputs z˜ = [ x F, F ] T

4.9.2 Scalar values B

[mol/s]

Bottom product flow rate

D

[mol/s]

Top product flow rate

F

[mol/s]

Feed flow rate

h' j

[J/mol]

Molar enthalpy of liquid phase on tray j

h'' j

[J/mol]

Molar enthalpy of vapor phase on tray j

L0

[mol/s]

Reflux

Lj

[mol/s]

Rate of liquid flow leaving tray j

nj

[mol]

Holdup on tray j

nt

[–]

Number of trays in column

T

[K]

Temperature

TPj

[K]

Pressure compensated temperature

V

[mol/s]

Boilup

Vj

[mol/s]

Rate of vapor flow leaving tray j

4.10 References

103

xj

[mol/mol]

Liquid phase composition on tray j

xB

[mol/mol]

Composition in column bottom

xF

[mol/mol]

Feed composition

yj

[mol/mol]

Vapor phase composition on tray j

y k*

[mol/mol]

Equilibrium vapor phase composition on tray j

α

[–]

Relative volatility

η

[–]

Murphree tray efficiency

σ

[–]

Singular value

κ

[–]

Condition number

4.10 References [4.1] Balas, G. J., J. C. Doyle, K. Glover, A. Packard, and R. Smith: μAnalysis and Synthesis Toolbox (μ-Tools), The MathWorks, Inc., Natick, MA (1991) [4.2] Chiang, R. Y., and M. G. Safonov: Robust Control Toolbox, The Mathworks, Inc., Natick, MA (1988) [4.3] Glover, K.: “All optimal Hankel-norm approximations of linear multivariable systems and their L∞ error bounds,” Int. J. Control, 36, 1115-1193 (1984) [4.4] Häggblom, K. E.: “Modeling of Flow Dynamics for Control of Distillation Columns,” Proc. 1991 American Control Conference, Boston, USA (1991) [4.5] Jacobsen, E. W., P. Lundström, and S. Skogestad: “Modelling and Identification for Robust Control of Ill-Conditioned Plants — a Distillation Case Study,” Proc. 1991 American Control Conference, Boston, USA (1991)

104

4 Linear Models

[4.6] Kapoor, N., and T. J. McAvoy: “An Analytical Approach to Approximate Dynamic Modeling of Distillation Towers,” IFAC Control

of

Distillation

Columns

and

Chemical

Reactors,

Bournemouth, UK (1986) [4.7] Moore, B.C.: “Principal Component Analysis in Linear Systems: Controllability, Observability and Model Reduction,” IEEE Trans. Automatic Control, 32, 115-122 (1981) [4.8] Safonov, M.G., R. Y. Chiang, and D. J. N. Limebeer: “Hankel Model Reduction without Balancing — A Descriptor Approach,” Proc. IEEE Conf. on Decision and Control, Los Angeles, CA, Dec. 9-11 (1987) [4.9] Safonov, M. G. and R.Y. Chiang: “Schur Balanced Model Reduction,” Proc. American Control Conference, Atlanta, GA, June 1517 (1988) [4.10] Skogestad, S. and M. Morari: “The Dominant Time Constant for Distillation Columns,” Comp. Chem. Eng., 11, 6, 607-617 (1987) [4.11] Skogestad, S. and M. Morari: “Understanding the Dynamic Behavior of Distillation Columns,” Ind. Eng. Chem. Res., 27, 18481862 (1988) [4.12] Skogestad, S.: “Dynamics and Control of Distillation Columns — A Critical Survey,” 3rd IFAC Symp. on Dynamics and Control of Chemical Reactors, Distillation Columns, and Batch Processes, April 26-29, College Park, MD, USA (1992) [4.13] Troch, I., P. C. Müller, and K.-H. Fasol: “Modellreduktion für Simulation und Reglerentwurf,” at, 40, 2, 45-53 (1992)

5.1 Introduction

105

Chapter 5 A Structured Uncertainty Model

5.1 Introduction Each linear or nonlinear dynamic model can only approximately describe the behavior of a real distillation column. While a nonlinear model may be valid for a wide range of operating conditions, the error of a linear model rapidly increases with the distance from its steady-state design point due to process nonlinearity. Since stochastic effects influence the process behavior as well, the error of a linear model compared to the real process can never be exactly determined. Lacking an exact error description, the error between the process model and the process itself is modelled as a single frequency-dependent uncertainty bound (unstructured uncertainty) or as several frequency-dependent uncertainty bounds (structured uncertainties). Typical sources of uncertainty for a distillation column are measurement errors, limited actuator speed, unmodelled high-frequency dynamics, and process nonlinearity. All these sources of uncertainty occur simultaneously and can be classified into three different groups: • Uncertainty of the manipulated variables (input uncertainties) • Model uncertainty due to process nonlinearity and unmodelled high-frequency dynamics

106

5 A Structured Uncertainty Model

• Uncertainty of the temperature measurements (output uncertainties) This grouping corresponds to the principle that uncertainty should be modelled where it physically occurs. In this chapter an uncertainty model for the industrial distillation column is developed. The complete uncertainty model covers not only a single operating point but the entire operating range of the column. It is the basis for the analysis and synthesis of controllers using the structured singular value μ.

5.2 Limits of uncertainty models Before we start to model the uncertainty in the frequency domain, we must be conscious of its limits: An uncertain model in the frequency domain is a model, which is time-invariant, but uncertain in its coefficients. This statement is best explained by an example: Let us model a ±10% uncertainty at the input of any plant and design a controller which guarantees closed-loop stability and a certain performance for all plants within the specified bounds. Then the stability and performance properties of the controller are not guaranteed for a time-varying plant, that means e.g. for variations of the input error between -10% and +10%. Consequently, using uncertainty models in the frequency domain, the excitation of the controller by the time-variation of the plant is not taken into account. If time-varying uncertainties are assumed, nonlinear simulations must be used for a validation of the robustness properties. However, the experience shows that for most cases uncertainty descriptions with frequency dependent and hence time-invariant uncertainty bounds are sufficient. This holds especially for our distillation column: The main disturbances are step changes of the feed flow rate. Each step change alters the steady state operating point and defines a new linear model describing the dynamic behavior up to the next step change. Each of these linear

5.3 Input uncertainty

107

models is one of the models within the set of all models. This set is defined by the specified uncertainty bounds.

5.3 Input uncertainty The actual values of the manipulated variables reflux and boilup will never match exactly the values requested by the control system. The error between the setpoints for the boilup or the reflux and the true streams will be frequency dependent. The main causes are • static and dynamic measurement errors of reflux and reboiler duty • changing heat of evaporation due to pressure and temperature variations • reboiler lags • actuator lags • effects of sampling The bounds for the relative errors of the column inputs u can be modelled by a multiplicative uncertainty description with the frequency-dependent error bound w u for the reflux L and the error L bound w u V for the boilup V. These bounds are combined in the diagonal matrix W u . As illustrated in Figure 5.1 the following uncertainty model holds: u˜ ( jω ) = { I + Δ u ( jω )W u ( jω ) }u ( jω )

Wu u

(5.1)

Δu u˜

Figure 5.1: Multiplicative uncertainty description for column input

108

5 A Structured Uncertainty Model

with Δ u ( jω )

W u ( jω ) =



≤1

(5.2)

w u ( jω )

0

0

w u ( jω )

L

(5.3)

V

The frequency-dependent complex matrix Δ u ( jω ) is limited in magnitude. It shapes only the spatial direction of the error and is chosen to be the worst case during μ-analysis (see Chapter 6). Therefore the phase behavior of the individual uncertainty bounds w u i is not significant. They should be chosen to be stable and minimum phase. If we assume that the reflux and the boilup errors are independent, the matrix Δ u ( jω ) becomes a diagonal matrix with two single perturbations δ yielding the following uncertainty model: ⎧ ⎫ 0 δ u ( jω ) ⎪ ⎪ L W u ( jω ) ⎬u ( jω ) u˜ ( jω ) = ⎨ I + 0 δ u ( jω ) ⎪ ⎪ V ⎩ ⎭

(5.4)

with δ u ( jω ) i



≤1

(5.5)

Both models have been used for μ-synthesis with very similar results. For two reasons, the model (5.1) is preferred in this study: • The number of uncertainty blocks (Δj or δi) is reduced by one compared to (5.4). This simplifies the μ-synthesis. • Any change in reflux may cause a change of the vapor flow rate within the column and vice versa. The interactions due to flow dynamics and to the energy balance are to be considered here.

5.3 Input uncertainty

109

Shaping the input error bounds It has been shown by Skogestad et al. [5.4] that the controlled system’s performance for a high purity distillation column is very sensitive to errors in the manipulated variables. For controller design or analysis the error bounds W u should be estimated as exactly as possible. This holds especially for the low-frequency range, where the condition number of the column models is high. Otherwise potential controller performance is given away in case of an overestimation. In the lower-frequency range the errors of the manipulated variables at the plant input are strongly dominated by flow measurement errors and parameter variations. As an example for a parameter variation we consider the heat of evaporation in the reboiler. The boilup is controlled indirectly by steam flow rate. Therefore a change in heat of evaporation will cause an error in vapor flow rate leaving the reboiler. Skogestad and Morari [5.4] assume a conservative 20% error in steady state, which is fairly high. If all flow measurements are carefully calibrated the error bounds should be less. An error bound of 10% for the lower frequency range is assumed to be large enough. The effects of reboiler lags, actuator lags, dynamic measurement errors, and sampling time concern the higher frequency range. The errors caused by these uncertainty sources increase with the frequency and easily exceed more than 100% of the nominal value for frequencies above 0.5-1 rad/min. The steady-state error, together with the high-frequency error, is well approximated by the first order lead/lag transfer function 1 + s ⁄ ωN G ( s ) = K -----------------------1 + s ⁄ ωD

with ω N < ω D

(5.6)

The gain K represents the steady-state error. The cut-off frequencies are typically chosen according to ω D > 10ω N .

110

5 A Structured Uncertainty Model

5.4 Model uncertainty 5.4.1 Column nonlinearity The highly nonlinear behavior of distillation columns is observed at varying operating points (varying feed flow rate and feed composition) and at transients during disturbance rejection. If we consider the simplified composition dynamics (without feed or side-product stream) of a tray (see (4.6)) dx j n j ⎛ --------⎞ = L j – 1 ( x j – 1 – x j ) + V j + 1 ( y j + 1 – x j ) – V j ( y j – x j ) ⎝ dt ⎠

(5.7)

we recognize that the composition dynamics and thus the nonlinear behavior depend on • the varying internal flow rates (L and V), and • the composition profile within the distillation column (represented by the liquid and vapor phase compositions) Effect of varying operating points Any control system for a distillation column must exhibit large gains in the low-frequency range to achieve small control errors at steady-state. Therefore, at steady-state both product compositions (or the temperatures on trays 10 and 44) can be kept at their setpoints. Thus transients have no significant influence in the low-frequency range and the internal vapor and liquid flow rates as well as the composition profile within a column become a function of feed flow rate and composition only. However, the dynamic behavior of a distillation column depends substantially on the actual composition profile and on the actual internal vapor and liquid flow rates. Normally the operating range of a distillation column can be bounded with a maximum and a minimum feed flow rate and composition. If we consider the whole operating range defined in this way, we can observe the largest internal flow rates for the

5.4 Model uncertainty

111

smallest feed composition and largest feed flow rate and, vice versa, smallest internal flow rates for the largest feed composition and smallest feed flow rates. The composition profiles for these two steady states bound the domain of all steady state composition profiles (see Figure 2.2). Hence we can conclude that the low-frequency behavior of a binary highpurity distillation column is bounded by the models for maximum and minimum column load. As a basis for further discussion the following three linear models are introduced: Model N

column at nominal load

Model I

column at maximum feed flow rate and minimum feed composition (increased load)

Model R

column at minimum feed flow rate and maximum feed composition (reduced load)

The feed data of the different models are listed in Table 5.1.

Table 5.1: Operating conditions for design purposes Operating point

Feed flow rate (mol/min)

Feed composition (mol/mol)

OP-N

33

0.8

OP-I

46

0.7

OP-R

20

0.9

The simplest way to represent the column nonlinearity due to varying operating points would be by a multiplicative output uncertainty. If we assume that the uncertainty for each model output is independent of the actual value of the other two model outputs, the following form for the output uncertainty holds (Figure 5.2):

112

5 A Structured Uncertainty Model

δy

1

0

0 δy

Wy

0

0 2

0

0 δy

3



y

Figure 5.2: Multiplicative uncertainty at output

⎧ ⎫ δ y ( jω ) 0 0 ⎪ ⎪ 1 ⎪ ⎪ W y ( jω ) ⎬y ( jω ) y˜ ( jω ) = ⎨ I + 0 0 δ y ( jω ) 2 ⎪ ⎪ ⎪ ⎪ 0 0 δ y ( jω ) 3 ⎩ ⎭

(5.8)

with δy

i



≤1

and y ( jω ) = G N ( jω ) d ( jω ) u ( jω ) The transfer matrix Wy is a diagonal matrix with the uncertainty bounds for each output ( w y , w y , w y ) on the main diagonal. An upper 1 2 3 bound for these uncertainties can be obtained by a calculation of the standardized errors ΔG I ( jω ) and ΔG R ( jω ) for each channel u i → y j or d i → y j of the models G I ( jω ) and G R ( jω ) , respectively: – 1 ( jω ) ΔG I ( jω ) = [ G I ( jω ) – G N ( jω ) ]G N

(5.9)

5.4 Model uncertainty

113

– 1 ( jω ) ΔG R ( jω ) = [ G R ( jω ) – G N ( jω ) ]G N

(5.10)

The upper bound for the uncertainty weights w y is the maximum of all j standardized errors for the output y j . In earlier papers it has already been recognized that column nonlinearity is not well represented by simple multiplicative uncertainty bounds at model output (McDonald [5.2]). This fact is confirmed by the uncertainty bounds for the two numerically evaluated linear models GI and GR which include the flow dynamics (Figure 5.3). The multiplicative output uncertainty exceeds 80% (for G u → y ) in the low-frequency range. It is significantly smaller in the medium frequency range, but increases sharply for frequencies above 0.1 rad/min, where the flow dynamics influence the dynamic behavior. An uncertainty description with such a high multiplicative uncertainty in the low-frequency range is prohibitive for any control design. Fortunately the errors are highly correlated: The variation of the steady-state operating points causes a simultaneous increase or decrease of the singular values of the transfer functions from the control signals u (L and V) to the model outputs y (T10, T44, T24). This is illustrated by the Nyquist plots for the individual channels u i → y j (Figure 5.4). It clearly shows that the variation of the column load causes a simultaneous increase or decrease of the open-loop gains in the lowfrequency domain. Thus we can assume that the dynamic behavior of the distillation column must lie “somewhere between OP-I and OP-R.” This can be represented by a linear combination of the two column models G I ( jω ) and G R ( jω ) (Figure 5.5) G I ( jω ) – G R ( jω ) G I ( jω ) + G R ( jω ) G ( jω ) = ------------------------------------------- + δ G ( jω ) -------------------------------------------2 2 with δG



≤1

δ G ∈ C or δ G ∈ R

(5.11)

114

5 A Structured Uncertainty Model

Standardized error Standardized errorofofT-10, T10,OP-I GI(s)

1

0.5

Magnitude

Magnitude

1

0 -0.5

Standardized error Standardized errorofofT-44, T44,OP-I GI(s)

0.5 0 -0.5

-1 -5 10

-4

10

-3

10

-2

10

-1

10

10

-1 -5 10

0

Frequency (rad/min) 1

-4

10

-3

10

10

-2

-1

10

0

10

Frequency (rad/min)

Standardized error Standardized errorofofT-24, T24,OP-I GI(s)

Legend 1

Magnitude

xF 0.5

0.5

0

0

-0.5

-0.5

F L V -1 -5 10

-4

10

-3

10

-2

10

-1

10

10

-1 -5 10

0

0

10

Frequency (rad/min) Standardized error Standardized errorofofT-44, T44,OP-R GR(s)

1

0.5

0.5

Magnitude

Magnitude

Standardized error error ofofT-10, Standardized T10,OP-R GR(s)

1

0 -0.5

0 -0.5

-1 -5 10

-4

10

-3

10

-2

10

-1

10

10

0

-1 -5 10

-3

10

10

-2

-1

10

0

10

Frequency (rad/min)

Frequency (rad/min) 1

-4

10

Standardized error GR(s) Standardized errorofofT-24, T24,OP-R

Legend Legend 1

Magnitude

x

xFF 0.5

0.5 FF

0

0

-0.5

-0.5

LL VV -1 -5 10

-4

10

-3

10

-2

10

-1

10

10

0

-1 -5 10

0

10

Frequency (rad/min)

Figure 5.3: Standardized model errors at operating points OP-I and OP-R

5.4 Model uncertainty

115

G L →L-T10 T ( jω )

GV → T ( jω ) V-T10

10

10

0

-100

0

0

200

-100

0

G L →L-T44 T ( jω )

GV → T ( jω ) V-T44

44

44

0

-400

0

0

700

-300

0

G L →L-T24 T ( jω )

0

600

GV → T ( jω ) V-T24

24

0

-250

200

24

0

500

-250

0

Figure 5.4: Nyquist plots for different column loads ⎯ solid lines: Model N; dash-dotted lines: Model I; dashed lines: Model R; ×: ω = 1×10-3 rad/min;

o: ω = 1×10-4 rad/min

500

116

5 A Structured Uncertainty Model

d

GI

1/2 + δG I3 –

GR

+

+ +

y

1/2

u Figure 5.5: Uncertainty model due to nonlinearity in the low-frequency range

The uncertainty parameter δ G may be either complex or real. If we define it to be complex we allow a phase shift for all models between G R and G I , that means the points of all models in the defined set in the Nyquist plots and for a fixed frequency are not required to be on a straight line. In this way we generate a plant which covers the properties of the distillation column at low and at high feed compositions, and at different feed flow rates without introducing additional conservatism. It is impossible to model such a behavior with an unstructured uncertainty description! Effect of transients While an appropriate uncertainty model for different operating points requires a highly structured uncertainty description, the effect of transients is rather unstructured: During disturbance rejection, the compositions on tray 10, tray 44, and tray 24 as well as the product compositions will deviate from their steady-state values, caused by a movement of the composition profile toward one column end. Due to the nonlinear vapor/liquid equilibrium, the singular values of the transfer functions G u → y may change in different directions, e.g., towards higher j singular values of G u → T and lower singular values of G u → T . 10

44

Due to the high controller performance in the low-frequency range, transients do not affect the low-frequency range. However, they cause nonlin-

5.4 Model uncertainty

117

earity in the middle and higher frequency range, which can be described with a multiplicative uncertainty description as in equation (5.8). The uncertainty weights w y are chosen to have large singular values in the i higher-frequency range and low singular values in the low-frequency range. It is not possible to calculate these uncertainty bounds exactly. Each disturbance input will cause a variation of the operating point, but the magnitude of the deviation from a steady-state operating point cannot be predicted. The selection of appropriate transfer functions is discussed in Chapter 6. 5.4.2 Unmodelled dynamics It has been shown in Chapter 2 that flow dynamics affect the highfrequency behavior of distillation columns. If linear models which neglect the flow dynamics are used for control design, an appropriate uncertainty model is necessary. Most authors treat the effect of flow dynamics in the same way as the effect of an input time delay τ with 0 < τ < 1 minute ([5.1], [5.3], [5.4], [5.5]). The corresponding input uncertainty is often modelled with a multiplicative uncertainty, using a first order Padé approximation for the uncertainty bound ([5.4], [5.5]): 1 – τs ----2 – τs e ≅ --------------τs 1 + ----2

(5.12)

This uncertainty can be combined with the other input uncertainties (Chapter 5.3). Lundström et al. [5.1] point to the fact that some combinations of gain uncertainty and time delay uncertainty are not represented using simple uncertainty weights. They developed new and more complicated uncertainty bounds, which cover the whole domain of combined gain uncertainty and time delay uncertainty. However, the control design studies in this research (Chapter 6) show very good results using simple first-order weights for the input uncertainty description.

118

5 A Structured Uncertainty Model

Model uncertainty due to flow dynamics could be represented by a multiplicative output uncertainty, as well. This approach has the disadvantage that the uncertainty bounds can no longer be approximated by time delays.

5.5 Measurement uncertainty An additional source of uncertainty are the temperature measurements. The dynamic behavior of a temperature sensor is well approximated by a first-order lag: KT G T ( s ) = ------------------1 + sT T

(5.13)

The time constant T T of this transfer function depends on the temperature measurement position. While the time constant will usually be clearly below 1 minute if the sensor is placed in the liquid phase, we have to expect time constants up to 10 minutes if it is placed in the vapor phase. In the case of the industrial distillation column under investigation, a position in the liquid phase cannot be guaranteed because of the small head on the plates. Therefore we have to consider time constants up to 10 minutes for the control design. The gain K T of the sensor model G T depends on the sensor calibration and on the heat loss to the environment. The sensor can easily be calibrated with high accuracy. However, the dynamic effects of the heat loss due to variations of the environment temperature are difficult to estimate. They concern mainly the low-frequency range and cause a slow bias variation of the temperature measurements. This effect is comparable to variations of the setpoints for the control system. The stability of the control system is not affected if large bias variations are avoided. They would lead to product compositions which are very distinct from those at the design operating points. A good thermal isolation of the temperature sensors is thus recommendable. Because the μ-analysis and μ-synthesis guarantees the performance for the worst case, it is proposed to include the model for the temperature

5.6 Specification of the controller performance

119

sensors with a gain K T = 1 and a time constant T T = 10 min into the column model. Thereby further uncertainty blocks can be avoided. It will be easily recognized later, that shorter time constants TT will not endanger the closed-loop system’s stability due to the large output uncertainties w y specified for the controller design in the upper i frequency range.

5.6 Specification of the controller performance The uncertainty model discussed above is structured. A controller design or a robust performance analysis requires the framework of the structured singular value μ, which expects the disturbance inputs d (feed composition, flow rate), the reference inputs r, and the control error to be in a H ∞ -norm bounded set. This is illustrated by the frequency shaped plant in Figure 5.6.

p We(s)

d r

e

Wd(s)

+

K(s)

-

GΔ(s) T24 T10,T44

Figure 5.6: Performance specification for the uncertain plant

The uncertain plant G Δ ( jω ) describes the nonlinear behavior of the distillation column. The performance objective is defined as making the weighted control error e to be in the set

120

5 A Structured Uncertainty Model

sup

W e ( jω )e ( jω )

2

≤1

∀ω ∈ R +

(5.14)

d ≤1 r 2

The following H ∞ bound is equivalent to this specification: T

d →p r

( jω )

≤1

(5.15)



W e ( s ) is a (usually diagonal) matrix of transfer functions which shapes the maximum allowed amplitude of the transfer function from [ d, r ] T to e. If W e is large in a certain frequency range, only a small control error is allowed there. The matrix W d ( s ) shapes the frequency content of the disturbances and setpoint changes. In the case of our distillation column, variations of the feed composition and feed flow rate will affect the medium and lower frequency range. First order lags shape the frequency content of these two disturbances quite well. Because measurement noise enters the control loop at the same position as the reference inputs, the corresponding weights are chosen to be constant. The weighting functions chosen are discussed in the following chapter.

5.7 Summary The complete uncertainty model is shown in Figure 5.7. It consists of the input uncertainty (5.1), the model uncertainties (5.8) and (5.11), and the performance specifications. Simple dynamic models of the temperature sensors are included in the column models. This relatively complex uncertainty model has the advantage that the entire operating range of the distillation column is covered. The large conservatism of an unstructured uncertainty description is avoided. Therefore, with design procedures based on the μ-synthesis or μ-optimization, we can expect high controller performance for the entire operating range.

r

d

Wd

+



e

We

p

K

Wu

Δu +

+

ρi

~u

GR 1/2

1/2

-

+ ξ o δGΙ

Figure 5.7: Complete uncertainty model

ρo

GI +

+

ξi +

Wy

1

0

0

T24

0

0 δy 3

0 δy 0 2

δy

T10, T44

ηo

+

+

ηi

5.7 Summary 121

122

5 A Structured Uncertainty Model

The input uncertainty bounds w u i are easily shaped. Only a few reflections are necessary about the steady-state error and the frequency where a 100% error is to be expected. However, the output uncertainties w y j are more difficult to shape. During the controller design procedure it is often necessary to adjust them iteratively until nonlinear simulations show a satisfactory closed-loop dynamics. This problem is discussed further in Chapter 6.

5.8 References [5.1] Lundström, P., S. Skogestad, and Z.-Q. Wang: “Uncertainty Weight Selection for H-Infinity and Mu-Control Methods,” Proc. 30th Conference on Decision and Control, Brighton, U. K. (1991) [5.2] McDonald, K. A.: “Characterization of Distillation Nonlinearity for Control System Design and Analysis,” The Shell Process Control Workshop, ed. D. M. Prett and M. Morari, Butterworth, Boston, 279-290 (1987) [5.3] Postlethwaite, I., J.-L. Lin and D.-W. Gu: “Robust Control of a High Purity Distillation Column Using μ-K Iteration,” Proc. 30th Conference on Decision and Control, Brighton, U. K (1991) [5.4] Skogestad, S., M. Morari, and J. C. Doyle: “Robust Control of IllConditioned Plants: High-Purity Distillation,” IEEE Trans. Automatic Control, 33, 12, 1092-1105 (1988) [5.5] Skogestad, S., and P. Lundström: “Mu-Optimal LV-Control of Distillation Columns,” Comp. Chem. Eng., 14, 4/5, 401-413 (1990)

6.1 Introduction

123

Chapter 6 μ-Optimal Controller Design

6.1 Introduction While the synthesis and analysis of controllers using the structured singular value μ (SSV) has attracted considerable attention among aerospace and electrical engineers (e.g., [6.8], [6.9]), it has been less commonly considered by process control engineers. One reason for that might be the lack of adequate structured uncertainty models for chemical processes. The uncertainty model discussed in the previous chapter forms a suitable basis for a μ-optimal controller design. Since this uncertainty model covers the dynamic behavior of the industrial distillation column for the entire operating range, the resulting controllers guarantee stability and performance for all operating points. This chapter presents the results of a μ-optimal controller design for the LV control structure of the distillation column. After a summary of the most useful aspects of the SSV, the design of state-space controllers by μ-synthesis is demonstrated. Because the implementation of state-space controllers in a distributed control system is a troublesome project, the design of controllers with fixed and easy-to-implement structures (PID control structures) is considered in a special section. A comparison of the controller’s performances in the time-domain terminates this chapter.

124

6 μ-Optimal Controller Design

6.2 The structured singular value The uncertainty model approximating the nonlinear dynamic behavior of the industrial distillation column (see Chapter 5) includes several simultaneous uncertainty blocks ( δ i, Δ j ), thus representing a structured uncertainty model. Most of the well-known robust control design methods (e.g., H∞, LQG/LTR) are based on unstructured uncertainty descriptions. The application of these methods on such uncertainty models often introduces unnecessary conservatism in controller design, because these methods combine all the uncertainties in one large, fully occupied uncertainty block. Thus the special structure of the uncertainties is neglected. This conservatism can be avoided by the use of the structured singular value μ, which was introduced in 1982 by J. C. Doyle ([6.5], [6.6]). The structured singular value μ so far has seldom been discussed in textbooks. Therefore the most important facts about μ are summarized within the following three sections. The discussion is restricted to complex uncertainty blocks. Results for mixed real/complex uncertainties can be found in [6.15]. The references [6.4], [6.12], and [6.14] contain additional informations. 6.2.1 Representation of structured uncertainties The definition of the structured singular value presumes that the uncertainty model for a plant is rearranged into a special form, as shown in Figure 6.1. The plant P consists of the process models and the weighting functions. It has three sets of inputs and outputs: The first set of inputs and outputs is highly important. Within the uncertainty model, this set represent the output and input signals of the uncertainty blocks. In our case the inputs to the uncertainty blocks are the signals ρ o, ξ o and η o . The corresponding outputs are the signals ρ i, ξ i , and η i . The second set of inputs consists of all external signals (disturbances d, reference inputs r), while the third set of inputs consists of all manipulated inputs u. The corresponding set of outputs contains the outputs p,

6.2 The structured singular value

125

Δ

ρi

ρo

Uncertainties

ξi = z

v =

ηi

ξo ηo

P

d r

p

Plant u

y

K Controller Figure 6.1: Standard representation of an uncertain plant. The definition of the vectors z and v is related to Figure 5.7

subject to any performance measure (e.g., the weighted control error), and the measured plant outputs y, respectively. If the uncertainty model is structured (i.e., it contains more than one uncertainty block) the matrix Δ is a block diagonal matrix with all uncertainty blocks on the main diagonal. In case of the uncertainty model for the distillation column considered in this thesis, the following block structure holds: Δ = diag ( Δ u, δ G I 3, δ y , δ y , δ y Δ u ∈ C 2 × 2, δ G ∈ C, δ y ∈ C ) 1

2

i

3

or, alternatively, Δ = diag ( Δ u, δ G I 3, δ y , δ y , δ y Δ u ∈ C 2 × 2, δ G ∈ R, δ y ∈ C ) 1

with

2

i

3

Δj



≤ 1, δ j



(6.1)

≤1

For an unstructured uncertainty model, the matrix Δ is a fully occupied matrix without predefined structure.

6 μ-Optimal Controller Design

126

The rearrangement of an uncertainty model into the standard form is always possible. The MATLABTM μ-Analysis and Synthesis Toolbox [6.1] as well as the Robust Control Toolbox [6.3] provide efficient tools for that purpose. 6.2.2 Definition of the structured singular value Let X be the set of all Δ matrices with a given, fixed block-diagonal structure: ⎧ m × mj ⎫ X = ⎨ diag ( δ 1 I r , …, δ s I r , Δ 1, …, Δ f ) δ i ∈ C, Δ j ∈ C j ⎬ 1 s ⎩ ⎭

(6.2)

The structured singular value [6.7] of the Matrix M ∈ C m × m with m = ∑ r i + ∑ m j (Fig. 6.2) is defined by ⎧ 1 ⎪ --------------------------------------------------------------------------------------------⎪ min { σ max ( Δ ) ( det ( I + MΔ )=0 ) } ⎪ μΔ ( M ) = ⎨ Δ ∈ X ⎪ ⎪ ⎪ 0 if no ( Δ ∈ X ) solves det ( I + MΔ )=0 ⎩

(6.3)

Hence 1 ⁄ μ ( M ) is the size of the smallest matrix Δ which moves a pole of the system shown in Figure 6.2 onto the imaginary axis. In the case of a nominally stable system M, 1 ⁄ μ ( M ) is the size of the smallest destabilizing matrix Δ . In case of a nominally unstable system M,

Δ M Figure 6.2: M-Δ feedback connection

6.2 The structured singular value

however, μ ( M ) misleading.

127

is not defined, and the numerical results are

Some important properties of μ [6.7] Let D be the set of diagonal scaling matrices: D = { diag D 1 , …, D S, d S + 1 I m , …, d S + F I m 1

Di ∈ C

ri × ri

,

D i =D i*

F

(6.4)

> 0, d S + j ∈ R , d S + j > 0 }

and let U be the set of block-diagonal unitary matrices ⎧ ⎫ k ×k U = ⎨ diag ( U 1, U 2, …, U n ) U i ∈ C i i, U i* U i =I k ⎬ i ⎩ ⎭

(6.5)

With these definitions the following properties of μ hold: ρ ( M ) ≤ μ ( M ) ≤ σ max ( M )

(6.6)

μ ( DMD –1 ) = μ ( M )

(6.7)

max ρ ( UM ) ≤ μ ( M ) ≤ inf σ max ( DMD – 1 )

(6.8)

U∈U

D∈D

Property (6.6) reflects the advantage of the structured singular value μ: In the presence of structured uncertainty, usually the inequality holds. Therefore, μ is smaller than the maximum singular value. The invariance of μ to diagonal scaling is indicated by property (6.7), which is essential for the approximate calculation of the structured singular values as well as for the DK-iteration for μ-synthesis. No direct way has been found yet to calculate μ exactly. All algorithms for the numerical computation calculate upper and lower bounds according to property (6.8). Both bounds represent an optimization problem. The optimization problem for the upper bound is convex. For simple block structures with 2S + F ≤ 3 the upper bound is guaranteed to be equal to μ Δ ( M ) .

6 μ-Optimal Controller Design

128

The optimization problem for the lower bound is not convex and its calculation may converge to local maxima. Nevertheless numerical experience indicates that usually the difference between upper and lower bounds is within 5%, and almost always within 15% ([6.12], [6.14]). 6.2.3 Robustness of stability and performance Before we start to discuss robust stability and performance within the framework of the structured singular value, we have to join the plant P and the known controller K of the standard configuration in order to close the control loop (Fig. 6.3). This is easily done by a linear fractional transformation [6.13]: M ( P, K ) = F l ( P, K ) = P 11 + P 12 K ( I – P 22 K ) – 1 P 12

(6.9)

The resulting plant M has two sets of inputs and outputs:

M 11 M 12 v = p M 21 M 22

z (6.10)

d r

Δ z

d r

v

M(P,K)

p

Figure 6.3: Representation of uncertain control system with controller K and plant P combined into the system M

6.2 The structured singular value

129

The input sets are (1) the outputs from the uncertainty block Δ , and (2) the disturbance and reference inputs. The outputs, in turn, are the inputs to the uncertainty block Δ and the set of performance measures p. Theorem 6.1: Robust stability Let BX be the set of all block diagonal matrices with a particular structure and with infinity-norm-bounded submatrices: BX = { diag ( δ 1 I r , …, δ s I r , Δ 1, …, Δ f ) s

1

δ i ∈ C, Δ j ∈ C

mj × mj

⎫ , δ i ∞ ≤ 1, Δ j ∞ ≤ 1 ⎬ ⎭

(6.11)

The system shown in Figure 6.3 remains stable for all Δ ∈ BX if and only if sup μ Δ ( M 11 ) < 1 ω

(6.12)

Proof: see [6.6] Theorem 6.1 allows the stability analysis of control systems with structured uncertainties. If we plot the upper and lower bounds of μ ( M 11 ) for enough frequency points in the frequency range of interest and find that the maximum value of μ is smaller than one, the control system is stable for the uncertainties specified with the assumption δ i ∞ ≤ 1 , Δ j ∞ ≤ 1 . If μ ( M ) exceeds one for any frequency, the control system is not guaranteed to be stable. However, for smaller uncertainties with δ i ∞ ≤ 1 ⁄ ( sup μ Δ ) and Δ j ∞ ≤ 1 ⁄ ( sup μ Δ ) stability is guaranteed. ω

ω

Theorem 6.2: Robust performance The performance of the control system is robustly achieved if and only if sup μ ˜ ( M ) < 1 ω

Proof: see [6.6]

Δ

with

Δ = diag [ Δ, Δ P ]

(6.13)

6 μ-Optimal Controller Design

130

For the application of theorem 6.2 we have to add one uncertainty block Δ P to the uncertainty structure Δ (Fig. 6.4). Imagine that the performance specification of the control system is met for all allowed disturbance matrices Δ in the set BX. In this case the output p is bounded by p



≤ 1 for all inputs [ d T, r T ] T with [ d T, r T ] T



≤ 1 . If we close the

loop from p to [ d T, r T ] T by introduction of the block Δ P with Δ P the system will be stable. But if any block Δ P with Δ P





≤ 1,

≤ 1 destabilizes

the loop p → [ d T, r T ] T , the specified performance cannot be achieved for all possible plants within the specified set. Therefore, in the framework of the SSV, the robust performance problem is handled like a stability problem. A test for robust performance will be similar to a test for robust stability. Because the test for robust performance includes robust stability, it usually will be sufficient.

z

Δ 0 0 ΔP

v

d r

M(P,K)

p

Figure 6.4: The robust performance setup for the SSV-framework

6.3 The design model The design of μ-optimal controllers for the industrial distillation column is based on the uncertainty model developed in Chapter 5 (see Figure 5.7 on page 121). For this model, the weighting functions for the

6.3 The design model

131

• input uncertainties, • output uncertainties, • reference inputs, disturbance inputs, and • controller performance are to be selected. All weighting functions are chosen as diagonal matrices: W d ( s ) = diag [ w x ( s ), w F ( s ), w r ( s ), w r ( s ) ] F

10

44

W u ( s ) = diag [ w u ( s ), w u ( s ) ] L

(6.14)

(6.15)

V

W y ( s ) = diag [ w y ( s ), w y ( s ), w y ( s ) ]

(6.16)

W e ( s ) = diag [ w e ( s ), w e ( s ) ]

(6.17)

10

10

44

24

44

The selection of the weighting functions is primarily done on the basis of physical considerations: Feed disturbances: The variations of the feed composition and the feed flow rate will affect the lower frequency range. The frequency contents of these disturbances are modelled by first-order lags. Typical weights chosen here for the control design are 1 w x ( s ) = K x ---------------------- with K x = 0.1 mol/mol, T x = 180 min (6.18) F F1 + T F F x s F

1 w F ( s ) = K F ------------------- with K F = 6 mol/min, T F = 120 min 1 + TF s

(6.19)

Reference inputs: The reference inputs r can be used to model setpoint changes as well as measurement noise. They are chosen as constant weights, representing setpoint changes and measurement noise of ±0.2 °C:

6 μ-Optimal Controller Design

132

w r ( s ) = w r ( s ) = 0.2 10

44

(6.20)

Input uncertainties: An uncertainty of 10% is assumed for both manipulated variables within a wide frequency range. For higher frequencies the uncertainty is expected to be much higher. An uncertainty of more than 100% is assumed for ω > 0.5 rad/minute: 1 + 20s w u ( s ) = w u ( s ) = 0.1 ----------------------L V 1 + 0.02s

(6.21)

Output uncertainties: The resolution of the temperature measurements is limited. The maximum deviation of the pressure-compensated temperatures from their setpoints is usually significantly smaller than 1 °C. A measurement error for ΔT of 10% seems to be reasonable. In the higher frequency range the output uncertainty is affected by model mismatch. An assumption of a 100% error for ω ≈ 1 ⁄ 16 rad/minute has shown good results in controller design. Adjusting of this 100% crossover is one of the possibilities to influence the high-frequency behavior of the resulting controller. Typical uncertainty weights are 1 + 167s w y ( s ) = w y ( s ) = w y ( s ) = 0.1 ----------------------10 44 24 1 + 1.67s

(6.22)

Performance weights: The performance weights “punish” the control error in the frequency domain. These weights have been selected as first-order lags with a large steady-state gain, which forces nearly integrating behavior of the controller. The cut-off frequency of these weights is a matter of optimization: If the frequency is too high, robust performance cannot be achieved. On the other hand a cut-off frequency specification significantly lower than the maximum attainable frequency may lead to an unsatisfactory controller design. This holds especially for the μK-iteration which is discussed later. A typical performance specification, which allows a 0.01 °C steady-state offset, is given by 1 w e ( s ) = w e ( s ) = 100 --------------------------10 44 1 + 27800s

(6.23)

6.4 Controller design with μ-synthesis

133

All weights above are illustrated by Figure 6.5 Uncertainty weights

10

10

10 Magnitude

Magnitude

10

Input and performance weights

2

1

wy wu

0

-1

10 -4 10

10

10

2

we wF

0

wr wx

-2

F

-4

-2

0

10 -6 10

2

10 10 10 Frequency (rad/min)

-4

-2

10 10 Frequency (rad/min)

10

0

Figure 6.5: Weights for μ-synthesis

6.4 Controller design with μ-synthesis The objective of the μ-synthesis is the calculation of a stabilizing controller K without restriction on the controller order and its structure, which minimizes the SSV for all frequencies: K = arg

inf K stabilizing

μ ˜ ( M ( P, K ) ) Δ

(6.24)



As it is not possible to calculate the SSV exactly, the design task (6.24) is usually approximated by the upper bound for μ (6.8) K = arg

inf

inf

K stabilizing D ∈ D

σ max ( DM ( P, K )D – 1 )



(6.25)

The aim of the μ-synthesis is perfectly reached, if the maximum value of μ ˜ for the closed-loop system (Figure 6.4) is below one. Δ

6 μ-Optimal Controller Design

134

6.4.1 Synthesis algorithms The μ-synthesis is not a trivial task. Yet no algorithms have been developed which allow a one-step solution of the μ-synthesis problem (6.24). The known algorithms require the repeated calculation of an H ∞ problem, alternating with a scaling of the plant. These algorithms cannot guarantee convergence of the iteration. DK-Iteration The synthesis problem (6.25) is a simultaneous optimization problem of the frequency-dependent scaling matrices D and the controller K. Because no direct solutions exist, Doyle [6.7] proposes an iterative approach: If we keep the diagonal scaling matrices D constant, the miniinf σ max ( DMD –1 ) ∞ forms the convex H ∞ problem mization of D∈D

K = arg inf DF l ( P, K )D – 1 K



(6.26)

If we fix the controller K, equation (6.26) represents a convex optimization problem for the diagonal scaling matrices D. These scaling matrices are optimized by a μ-analysis of the closed-loop system: μ ˜ [ F l ( P, K ) ] ≈ inf σ max ( DMD –1 ) Δ

(6.27)

D∈D

The frequency-dependent scaling matrices D are approximated with ˜ ( s ) . Alternating the H controller stable, rational transfer functions D ∞ synthesis and the optimal scaling, convergence is achieved for most μsynthesis problems after several iterations. The iteration procedure is illustrated in Figure 6.6. The DK-iteration is finished either if the solution does not show any further improvement or if μ < 1 . However, convergence cannot be guaranteed: Both of the single optimization problems are convex, but not the overall optimization problem (6.25). The optimized scaling matrices D are an optimal solution for the local optimization problem (6.27), but they are not optimal for the global optimization problem (6.25). Therefore, the DK-iteration may converge to local minima.

6.4 Controller design with μ-synthesis

135

K 0 = arg inf F l ( P, K ) ∞ K

μ ˜ [ F l ( P, K 0 ) ] = inf σ max ( DMD – 1 ) D∈D Δ

˜ (s) Fit D(s) with stable, min. phase transfer functions D

˜ ( s )F ( P, K )D ˜ –1 ( s ) K 1 = arg inf D l ∞ K

μ ˜ [ F l ( P, K 1 ) ] = inf Δ

D∈D

σ max ( DMD –1 )

.............

Figure 6.6: DK-Iteration

μK-Iteration A new algorithm for μ-synthesis has been proposed by Lin et al. [6.11]. Instead of fitting the scaling matrices D, this algorithm is based on a fit of the frequency-dependent SSV with a stable rational transfer function. At each iteration step, the plant is premultiplied with a diagonal matrix of the μ-approximating transfer function. Thus the peaks of μ within the frequency range of interest are punished, and the algorithm tries to flatten the μ-curve. The convergence of this algorithm is not proved. The authors present “a reasoned argument for believing that the sequence will converge” [6.11]. A scheme of the μK-algorithm is

6 μ-Optimal Controller Design

136

K 0 = arg inf F l ( P, K ) ∞ K

μ ( jω ) = μ ˜ [ F l ( P, K 0 ) ] Δ

Fit

μ ( jω ) μ˜ 0 ( jω ) ≈ ----------------------μ ( jω ) ∞

, μ˜ 0 ( s ) stable, miminum phase

K 1 = arg inf μ˜ 0 ( s )F l ( P, K ) ∞ K

μ ( jω ) = μ ˜ [ F l ( P, K 1 ) ] Δ

Fit

μ ( jω ) μ˜ 1 ( jω ) ≈ ----------------------μ ( jω ) ∞

K 2 = arg inf μ˜ 1 ( s )μ˜ 0 ( s )F l ( P, K ) ∞ K

.......... Figure 6.7: μK-Iteration

given in Figure 6.7. This iteration scheme has certain disadvantages: It usually converges more slowly than the DK-iteration, and the convergence properties are strongly dependent on the fit of the μ-curve. Even if convergence is achieved, μ is not minimized for all frequencies.

6.4 Controller design with μ-synthesis

137

6.4.2 Applying the DK-Iteration The application of the DK-iteration to our design problem was not successful because convergence of the algorithm is not attainable. Most likely, the main problem is the fit of the fully occupied 2×2 or 3×3 (including or excluding the measurement of T24, respectively) block of the D-scaling matrices. This block results from the repeated scalar uncertainty block δ G I . If we fit each position of this block with a scalar, stable, and minimum-phase transfer function, a minimum-phase behavior for the resulting MIMO system is guaranteed. However, unnecessary conservatism is introduced thereby, since a minimum phase behavior is only required for the MIMO system, but not for the single scalar transfer functions. This problem remains to be solved. 6.4.3 Applying the μK-Iteration The application of the μK-Iteration does yield convergence. However, it is necessary to slightly modify the algorithm. The premultiplication of the μ-curve-fitting transfer function μ i ( s ) increases the order of the design model at each iteration step. This easily leads to models with more than 200 states. Unacceptable calculation times and numerical problems result therefrom. This problem is avoided by an order reduction step after the augmentation of the plant and before the H ∞ design. The order reduction method utilized is a balanced truncated realization [6.1]. Experiences A typical course of the iteration is shown by Figure 6.8. In the frequency range where the performance specification μ ˜ ( jω ) < 1 is not achieved, Δ the upper bound of the SSV is forced down at each iteration step. After six steps, the solution is reached for which no further improvement is possible. If we look at the frequency range between 10 –2 and 10 –1 rad/ min, we discover that the first controller exhibits much better robust performance than the final design. This results from the “flattening”

6 μ-Optimal Controller Design

138

Structured singular value (RP)

1 1.5

2 3 4 5

1

6

0.5

0 -6 10

-5

10

-4

10

-3

10

-2

10

-1

10

10

0

Frequency (rad/min) Figure 6.8: Convergence of the μK-Iteration

behavior of the μK-Iteration, which leads to any solution of the design task, but not necessarily to the optimal one. This behavior of the iteration scheme may lead to strange results. Even if the robust performance criterion is achieved, the simulation of the closed-loop behavior may exhibit an insufficiently damped oscillation. It has to be emphasized here that such an oscillation is consistent with the performance specification. This performance of the closed-loop system is specified in the frequency-domain rather than in the time-domain! Slightly increasing the performance requirement or the uncertainty specifications usually removes this problem. Another problem of the μK-iteration is the small convergence area: The performance specification W e ( s ) (see (6.17)) has to be close to the maximum achievable performance, otherwise the iteration does not converge. As a last drawback the long computation times have to be

6.4 Controller design with μ-synthesis

139

mentioned. The CPU time for a design usually exceeds 2h on a SUN SPARC 2 workstation! Analytically or numerically linearized models? In chapter 4 two different types of linear models have been developed. The main differences between these model types consist of the low frequency gains and the representation of flow dynamics. In the case of the analytically linearized models, the relative uncertainty in the low frequency range due to variation of the steady-state operating points is essentially smaller (see Chapter 5). For both types of models, state-space controllers have been designed. In order to achieve an acceptable controller design (oscillation free) with analytically linearized models, the low-frequency gains of the disturbance weights (6.18) and (6.19) must be approximately doubled. With respect to the higher frequency range the fact of unmodelled flow dynamics does not dominate the shape of the output uncertainty weights w y for the tray temperature T10 and T44. Both weights may be kept equal for both types of linear models. However, the uncertainty weight for the temperature measurement in the middle of the column T24 has to be increased for the analytically linearized models due to the unmodelled flow dynamics. In accordance with these adaptations of the weighting functions, the resulting state-space controllers yield nearly identical closed-loop behavior with a small advantage from using the numerically linearized models. Complex or real uncertainty block δ G Within the structured uncertainty models, the uncertainty δ G may be chosen as a complex or real uncertainty. The choice as a real uncertainty reduces the uncertainty for the entire frequency range. However, the resulting closed-loop behavior exhibits insufficiently damped oscillations. To avoid this problem, the performance weights and the disturbance weights must be increased. The state-space controller designed with the modified weights are not superior to the design obtained with a complex uncertainty block δ G .

6 μ-Optimal Controller Design

140

μK-Iteration results for three temperature measurements The input vector of the controller may consist either of the pressurecompensated tray temperatures T10 and T44, or of all three temperatures. The temperature T24 is close to the feed tray and its response to feed flow disturbances is faster than that of the other two temperatures. Therefore, an improved controller design should result from this additional temperature measurement. Several synthesis attempts have shown that it is possible to increase the performance specification by circa 30% up to 1 w e = 100 --------------------------i 1 + 20580s

(6.28)

Structured singular value

After the convergence of the μK-Iteration, the final controllers were reduced to an order 20 using a balanced truncated realization of the control system [6.1]. The μ-plot for the reduced-order controller (Figure 6.9) using the uncertainty, input, and performance weights (6.14)–(6.23) demonstrate the excellent robustness properties of this controller. 1 RP

0.5

0 -5 10

Figure 6.9: Robust performance and stability for the μ-optimal state-space controller (controller inputs: T10, T44, T24)

RS

-3

-1

10 10 Frequency (rad/min)

1

10

An analysis of the nominal closed-loop system (with plant model GN) is shown in Figure 6.10 and Figure 6.11. The singular values of the transfer function from the reference signals r to the controlled outputs (Fig. 6.10 a) y indicate a good set-point tracking for the frequency range

6.4 Controller design with μ-synthesis

Tr → y

1

10

0

10

-1

10

-2

10 -5 10

Td → y

1

Magnitude

10

Magnitude

141

0

10

-1

10

-2

-3

10

10

-1

10 -5 10

1

Frequency (rad/min)

10

-3

-1

10

1

10

Frequency (rad/min)

a)

10

b)

Figure 6.10: Singular values for the nominal closed-loop system with the μ-optimal state-space controller (controller inputs: T10, T44, T24) a) Transfer function from reference signals to controlled output signals b) Transfer functions from disturbance signals to controlled output signals Dash-dotted line: T F → y , solid line: T x → y F

10

Magnitude

10

10

10

1

0

-1

-2

-3

10 -5 10

-4

10

-3

10

-2

10

-1

10

0

10

1

10

Frequency (rad/min) Figure 6.11: Singular values of the sensitivity function at u for the nominal closedloop system with the μ-optimal state-space controller (controller inputs: T10, T44, T24)

142

6 μ-Optimal Controller Design

of interest. The singular values of the individual transfer functions from the two disturbance inputs to the controlled outputs (Fig. 6.10 b) show a maximum of the sensitivity to these disturbances in the midfrequency range. While in the high-frequency range the sensitivity is smaller due to the low-pass characteristics of the plant, the large controller gains cause an effective compensation in the low-frequency range. The plot of the sensitivity at u (Fig. 6.11) S u ( s ) = [ I + K ( s )G ( s ) ] –1

(6.29)

confirms the good robustness properties in the common unstructured uncertainty representation. The maximum value (≈1.6) guarantees a stability phase margin of at least 35° [6.4]. The simulation of step responses using the rigorous dynamic model described in Chapter 3 demonstrates the closed-loop behavior in the time-domain. Two disturbances are simulated: An increase of the feed composition from 0.8 to 0.9 mol/mol, and an increase in the feed flow by 3.6 mol/min. Figure 6.12 shows the top and product impurities as well as the control errors for these disturbances and for maximum and minimum feed flow rates. To estimate the sensitivity to errors in the manipulated variables, a 10% error of the controller outputs ΔL and ΔV for the same test bench has been simulated. The results are represented by the thin lines in Figure 6.12. The steady-state offsets of the product compositions are caused by controlling pressure-compensated temperatures on trays 10 and 44 instead of the product compositions. The simulation results confirm the good robustness properties, especially the low sensitivity to errors in the manipulated variables. At both operating points, the overshoot of the control error is small. For the compensation of the first disturbance — an increase in feed composition — reflux and boilup must be reduced. The second disturbance is an increase in feed flow rate, which has to be compensated by an increase in reflux and boilup (see Figure 6.13). The large difference between reflux and boilup even at steady-state has various reasons. First, the reflux and feed are subcooled and a partial condensation of the

6.4 Controller design with μ-synthesis

Ft=0=20 mol/min F(t=0)=20 mol/min

0.015

0.010

0.005 0

0.4

Top composition Bottom composition

10

20 30 Time (h) [h]

0.020 Composition (mol/mol) Composition

Composition Composition (mol/mol)

0.020

143

40

0.010

0.4

Top composition Bottom composition

10

20 30 Time [h] (h)

40

Ft=0=46 mol/min F(t=0)=46 mol/min

0.2 Temperature (K) Temperature

0.2 Temperature (K) Temperature

0.015

0.005 0

Ft=0=20 mol/min F(t=0)=20 mol/min

0.0 -0.2 -0.4 -0.6 0

Ft=0=46 mol/min F(t=0)=46 mol/min

Control error T-44

20 30 Time (h) [h]

-0.2 -0.4

Control error T-10

10

0.0

40

-0.6 0

Control error T-10 Control error T-44

10

20 30 Time [h] (h)

40

Figure 6.12: Simulation results with μ-optimal state space controller (controller inputs: T10, T44, T24) for an increase in feed composition (0.8 → 0.9 mol/mol) at t=0 h and an increase of feed flow rate (+ 3.6 mol/min) at t=20 h Upper plots: Product composition Lower plots: Control error L, V equal to controller output ΔL with +10% error, ΔV with –10% error

6 μ-Optimal Controller Design

144

70

Ft=0=20 mol/min F(t=0)=20 mol/min

Ft=0=46 mol/min F(t=0)=46 mol/min 140

Flow (mol/min) Flow rate [mol/min]

Flow (mol/min) Flow rate [mol/min]

60

50 Boilup Reflux

40

120 Boilup

100 Reflux

30

80

20 0

60 0

10

20

Time Time (h) [h]

30

40

10

20

30

40

Time Time (h) [h]

Figure 6.13: Simulation results with μ-optimal state-space controller (controller inputs: T10, T44, T24) for an increase in feed composition (0.8 → 0.9 mol/mol) at t=0 h and an increase of feed flow rate (+ 3.6 mol/min) at t=20 h L, V equal to controller output ΔL with +10% error, ΔV with –10% error

vapor flow thus increases the liquid flow rates below the corresponding trays. Second, the major part of the feed leaves the column as top product. If we compare the plots for the minimum and maximum feed flow rate, we recognize an essentially slower rejection of the feed composition disturbance at the maximum feed flow rate. A distinct improvement of the performance at maximum feed flow rate is not possible using a linear time-invariant feedback controller. Higher controller gains would improve the disturbance compensation at this operating point, but simultaneously destabilize the control loop at the minimum feed flow rate. A closer look at these figures demonstrates that especially at high feed flow rates the controller response is more sluggish for changes in

6.4 Controller design with μ-synthesis

145

feed composition but not for disturbances in the feed flow rate. This fact is explained by the course of the manipulated variables in Figure 6.13. An increase in feed composition at minimum feed flow rate forces the controller to reduce the reflux and the boilup by ≈11 mol/min, while at maximum feed flow both flow rates have to be reduced by ≈30 mol/min! Since in practice a step change of feed composition is improbable, the rejection of feed flow variations has much higher significance. μK-Iteration results for two temperature measurements When only the pressure-compensated temperature on trays 10 and 44 are used as the controller inputs, it becomes extraordinarily difficult to achieve convergence of the μK-Iteration and an oscillation-free closedloop dynamics. For design purpose, the same weights as in the previous design for all three temperature measurements have been used.

Structured singular value

The final controller was reduced to order 20 by a balanced truncated realization of the control system [6.1]. The μ-plots of the reduced order controller (Figure 6.14) shows worse robustness properties of this controller in the higher frequency range (compared to the controller with three measured temperatures as input), only just matching the robustness and stability criteria.

1 RP

0.5

0 -5 10

Figure 6.14: Robust performance and stability for the μ-optimal state-space controller (controller inputs: T10, T44)

RS

-3

-1

10 10 Frequency (rad/min)

1

10

An analysis of the sensitivity functions (Figure 6.15) exhibits a small maximum sensitivity at the control error e, but an evidently reduced

6 μ-Optimal Controller Design

146

Sensitivity at e 1

10

0

Magnitude

10

-1

10

-2

10

-3

10 -5 10

10

-4

-3

10

-2

10

-1

10

0

10

1

10

Frequency (rad/min) Sensitivity at u

1

10

0

Magnitude

10

-1

10

-2

10

-3

10 -5 10

10

-4

-3

10

-2

10

-1

10

0

10

1

10

Frequency (rad/min) Figure 6.15: Singular values of the sensitivity functions at e (upper plot) and at u (lower plot) for the nominal closed-loop system with the μ-optimal state-space controller (controller inputs: T10, T44)

6.4 Controller design with μ-synthesis

147

stability margin at u. This illustrates the direct relationship between the SSV and the common unstructured robustness measures. Nevertheless the simulation results in Figure 6.16 demonstrate a high controller performance, paired with a larger sensitivity to input errors. If we compare the controllers with 3 temperatures and 2 temperatures as input, we must state that the “control qualities” in the time-domain are very similar. As mentioned before, the intuition of a control engineer is to expect a better performance for more measurements due to a faster state-estimation. This is obviously not the case! It will be possible to give an explanation for this result in the further course of this chapter.

0.020

Composition (mol/mol) Composition

Composition (mol/mol) Composition

0.020

Ft=0=20 mol/min F(t=0)=20 mol/min

0.015

0.010

Ft=0=46 mol/min F(t=0)=46 mol/min

0.015

0.010

Top composition

Top composition

Bottom composition

0.005 0

10

20 30 Time (h) [h]

Bottom composition

40

0.005 0

10

20 30 Time (h) [h]

40

Figure 6.16: Simulation results with μ-optimal state-space controller (controller inputs: T10, T44) for an increase in feed composition (0.8 → 0.9 mol/mol) at t=0 h and an increase of feed flow rate (+ 3.6 mol/min) at t=20 h L, V equal to controller output ΔL with +10% error, ΔV with –10% error

6 μ-Optimal Controller Design

148

6.5 Design of controllers with fixed structure In the process industry PID or advanced PID control structures are very common. Therefore, the implementation of a controller design and its acceptance are substantially improved if the design is based on PID control structures. The corresponding design objective is the μ-optimal tuning of simple control elements (such as PID controllers, first-order lags) within a fixed control structure: K = arg

inf K stabilizing K with fixed structure

μ ˜ (M) Δ



(6.30)

The solution of this design objective is extremely difficult. Because no synthesis methods exist, (6.30) must be solved by a parameter optimization approach. During this optimization the SSV has to be calculated repeatedly for a number of frequency points. However, the maximum of the SSV may be very sensitive to the number of frequency points calculated. In order to simplify the numerical treatment, the design objective can be approximated by a summation of the cube of the SSV for all k frequency points: k

θ = arg inf θ

∑ μΔ3˜ { Fl [ P, K ( θ ) ] }

(6.31)

i=1

Summing the cube, large values of the SSV have much more weight and the design objective becomes closer to (6.30). The calculation of the SSV presumes nominally stable control loops. Within μ-synthesis, the controllers are calculated by solving an H ∞ problem, which always guarantees nominal stability. However, during a parameter optimization, nominally unstable control loops may be generated. Therefore, the design objective (6.31) must be supplemented with the boundary condition for nominal stability: Re [ λ max ( F l { G i, K } ) ] < 0

(6.32)

6.5 Design of controllers with fixed structure

149

A second boundary condition is the robust stability criterion, which should be fulfilled for the final parametrized controller: μ Δ { F l [ P, K ( θ ) ] }



<1

(6.33)

This constrained parameter optimization problem is solved by sequential quadratic programming [6.10]. In contrast to the μ-synthesis methods, this approach has shown a high reliability, at the price of even higher computation times. However, the excellent results justify the effort. 6.5.1 Diagonal PI(D) control structures The diagonal PI(D) control structure (Figure 6.17) is the simplest and most frequently used composition control structure for distillation columns. Due to the high interaction between the two control loops, this structure is difficult to tune, and the response to setpoint changes or disturbances is known to be very sluggish.

– r10

PID1

+ r44

PID2

+



L

V

Distillation Column with inventory control

T10

T44

Figure 6.17: Diagonal PID control structure

The design model The design model for this optimization is the same uncertainty model as that used for the μ-synthesis, excluding the temperature measurement

6 μ-Optimal Controller Design

150

on tray 24. The weighting functions are the same transfer functions as discussed in section 6.3. Results for PI control The matrix transfer function of the diagonal PI control structure is given by

L(s) = V(s)

1 + TI 1 s KR 1 --------------------TI 1 s 0

0

e 10 ( s )

1 + TI 2 s e 44 ( s ) KR 2 --------------------TI 2 s

(6.34)

Table 6.1 summarizes the results of the parameter optimization for the analytically and numerically linearized column models, as well as for a complex and mixed real/complex μ-analysis. A comparison of the different optimization results shows quite similar parameters for a complex μ-analysis and a mixed real/complex μ-analysis. However, a significant difference exists between the numerically and the analytically linearized models: Using the analytically linearized models, the time constants TI 1 are much smaller and the corresponding low-frequency gains are much higher. The reason for that are the smaller low-frequency gains of these linear models. An underestimation of low-frequency uncertainty results therefrom. Simulations with these controller designs show a faster, but insufficiently damped closed-loop behavior. Of course, with an increase in the output uncertainty of the Table 6.1: Results for the diagonal PI control structure Model lineari- μ-analysis KR1 zation (mol/min/°C) Numerical

Analytical

TI1 (min)

KR2 (mol/min/°C)

TI2 (min)

Complex

–14.09

137

2.49

34

Mixed R/C

–11.27

141

3.15

52

Complex

–11.52

49

6.92

56

Mixed R/C

–14.58

60

6.36

41

6.5 Design of controllers with fixed structure

151

design model, the design can be improved. This leads to results comparable to those obtained with the numerically linearized models. These experiences corresponds to those of the μ-synthesis. Consequently, the design with the analytical linearized models is not further discussed. Due to the extremely large computation times using mixed real/complex μ-analysis and the very similar optimized tuning, the further discussion will focus on the optimization with complex μanalysis and numerically linearized models. The upper bounds for robust stability and performance (numerically linearized models, complex μ-analysis) using the diagonal PI control structure are shown in Figure 6.18. While stability is guaranteed for the specified uncertainties and for the entire frequency range, the performance specification is not met in the lower frequency range. However,

Structured singular value

robust performance is achieved within the upper frequency range.

2 1.5

RP

Figure 6.18: Robust performance and stability for diagonal PI control

1 0.5 0 -5 10

RS

-3

10

-1

10

1

10

Frequency (rad/min)

The transfer functions from the reference and disturbance inputs to the temperature measurements for the nominal closed-loop systems (Figure 6.19) shows a high condition number for the tracking behavior within the most important frequency range. This means a high sensitivity of the tracking behavior to the direction of the reference inputs.

6 μ-Optimal Controller Design

152

1

Tr → y

10

0

Magnitude

Magnitude

10

10

-1

10

-2

10 -5 10

Td → y

1

0

10

-1

10

-2

-3

10

-1

10

10 -5 10

-3

-1

10

10

Frequency (rad/min)

Frequency (rad/min)

a)

b)

1

10

Figure 6.19: Singular values for the nominal closed-loop with diagonal PI-controller a) Transfer function from reference to output signals b) Transfer functions from disturbance to output signals Dash-dotted line: T F → y , solid line: T x → y F

These conclusions are confirmed by the results of the nonlinear simulations (Figure 6.20). They demonstrate the sluggish disturbance rejection of the optimally tuned diagonal PI control. However, as expected from the robust performance plot, the maximum control error is sufficiently small. Another positive result is the small sensitivity to input uncertainty. Results for PID control The use of real PID control instead of PI control gives additional degrees of freedom for the controller design. Since true differential behavior is not realizable, the parameters for PID controllers with a first order lag in series are optimized (real PID controllers). The following transfer function for the controllers holds:

6.5 Design of controllers with fixed structure

0.015

0.010

0.005 0

0.4

Top composition Bottom composition

10

20 30 Time (h) [h] Time

0.020

Composition (mol/mol) Composition

Composition Composition (mol/mol)

0.020

Ft=0=20 mol/min F(t=0)=20 mol/min

0.010 Top composition Bottom composition

10

20 30 Time [h] (h)

40

Ft=0=46 mol/min F(t=0)=46 mol/min

0.2

Temperature (K) Temperature

Temperature (K) Temperature

0.015

0.4

0.2 0.0 -0.2

-0.6 0

Ft=0=46 mol/min F(t=0)=46 mol/min

0.005 0

40

Ft=0=20 mol/min F(t=0)=20 mol/min

-0.4

153

Control

error T-10

Control

error T-44

10

20 30 Time [h] (h)

0.0 -0.2 -0.4

40

-0.6 0

Control error T-10 Control error T-44

10

20 30 Time [h] (h)

40

Figure 6.20: Simulation results with diagonal PI control for an increase in feed composition (0.8 → 0.9 mol/mol) at t=0 h and an increase of feed flow rate (+ 3.6 mol/min) at t=20 h Upper plots: Product composition Lower plots: Control error L, V equal to controller output ΔL with +10% error, ΔV with –10% error

6 μ-Optimal Controller Design

154

e 10 ( s ) L ( s ) = G K1 0 V(s) 0 G K2 e 44 ( s )

(6.35)

1 + TI i s + TI i TD i s 2 G Ki ( s ) = KR i ---------------------------------------------------TI i s ( 1 + sTL i )

(6.36)

with

The optimal tuning results show unacceptably large controller gains in the high-frequency range. A high amplification of the measurement noise can be avoided by various methods: • “Punishment” of high frequency controller output by additional weighting functions • Limitation of high frequency gains by additional boundary conditions In order to keep the uncertainty model invariant, the differential behavior of the controller was limited by a minimum bound of 2 min for the time constants TL of the first-order lags. The resulting tuning constants are given in Table 6.2. Table 6.2: Results for the diagonal real PID control structure Controller

KR mol/min/°C

TI (min)

TD (min)

TL (min)

PID 1

–15.97

101

7.41

2.00

PID 2

4.40

39.0

15.2

7.16

Results achieved with numerically linearized model and complex μ-analysis

The μ-plots (Figure 6.21) for the diagonal PID control structure show an improvement of the robust performance. However, the design objective of robust performance

Structured Singular Value

6.5 Design of controllers with fixed structure

155

1.5 RP Figure 6.21: Robust performance and stability for diagonal real PID control

1

0.5

0 -5 10

RS

-3

-1

10 10 Frequency (rad/min)

1

10

μ ˜ { F l [ P, K ( θ ) ] } Δ



<1

(6.37)

is by far not reached. The simulation results in Figure 6.22 illustrate the same sluggish behavior as was obtained with the PI controllers. The improvement is a slight reduction of the settling time. An analysis of the controller’s singular values shows large highfrequency gains despite a limitation of the minimum filter time constant TL (Figure 6.23). This makes a first-order filter for the reference inputs

Magnitude

10

10

2

Figure 6.23: Singular values of the diagonal PID controller

0

-2

10 -5 10

-3

-1

10 10 Frequency (rad/min)

1

10

necessary. A reduction of these high-frequency controller gains would be

6 μ-Optimal Controller Design

156

Ft=0=20 mol/min F(t=0)=20 mol/min

0.020

Composition (mol/mol) Composition

Composition (mol/mol) Composition

0.020

0.015

0.010

Ft=0=46 mol/min F(t=0)=46 mol/min

0.015

0.010

Top composition

Top composition

Bottom composition

0.005 0

10 20 30 Time (h) [h] Time

Bottom composition

40

0.005 0

10

20 30 Time [h] Time (h)

40

Figure 6.22: Simulation results with diagonal PID control for an increase in feed composition (0.8 → 0.9 mol/mol) at t=0 h and an increase of feed flow rate (+ 3.6 mol/min) at t=20 h L, V equal to controller output ΔL with +10% error, ΔV with –10% error

possible. However, decreasing high-frequency limits annihilate the improvements achieved over the diagonal PI control structure. Summarizing the results for the diagonal PI(D) control structure, we can conclude that this control structure is absolutely not suited for a high performance. 6.5.2 PI(D) control structures with two-way decoupling The major disadvantage of the diagonal PI(D) control structures is the neglect of the interactions between the two control loops. These interactions can be partially cancelled by use of decoupling techniques.

6.5 Design of controllers with fixed structure

157

A simple controller structure with decoupling is shown in Figure 6.24.



r10 +

+

Distillation Column

C2 C1

r44

T10

L

PID1

PID2

V

with inventory control

T44



Figure 6.24: PID control structure with static decoupling

The decoupling elements can be static (static decoupling) or dynamic (dynamic decoupling). The tuning of the decoupling control structure for a distillation column is difficult. Often decouplers are based on an inversion of the plant’s transfer function G ( s ) . The resulting closed-loop behavior is very sensitive to input uncertainty and decoupler errors. Summarizing the research results, Skogestad in [6.17] concludes that (two-way) decouplers should never be used for high-purity distillation columns with the LV-configuration. On the other hand one-way decoupling seems to be less sensitive to input uncertainty and should be preferred [6.18]. Results for static two-way decoupling The simplest decoupling structure is static decoupling. Here the two decoupling elements C 1 and C 2 are constant factors. The results for this structure are obtained with the same weighting functions and with the same uncertainty model as used for the diagonal PI(D) control structure.

6 μ-Optimal Controller Design

158

Table 6.3 summarizes the μ-optimal parameters for PI and real PID Table 6.3: μ-optimal parameters for PI(D) control with static decoupling Controller or KR decoupler No. (mol/min/°C)

TI (min)

TD (min)

TL (min)

C (–)

1

–5.21

22.8





–0.0240

2

3.71

46.8





1.11

1

–13.1

51.6

7.83

8.43

–0.217

2

4.56

62.1

5.11

3.07

1.03

Results achieved with numerically linearized model and complex μ-analysis

control with static decoupling. The high-frequency gains of the PID controller are small enough that no boundary conditions concerning this criterion were necessary. The results for the decouplers are somewhat surprising. They indicate that the optimal decoupling is very close to a one-way decoupling! Let us examine this control structure in detail: The decoupler parameter C 2 is close to one. Therefore any variation of the output of the top composition controller causes a simultaneous increase or decrease of reflux and boilup by almost the same magnitude. Thus this controller shapes the composition profile within the column by an adaptation of the separation. The other decoupler parameter C 1 is small. Consequently the output of the bottom composition controller has a small effect on the reflux. This controller moves the composition profile within the column. In light of this interpretation, the limited advantage of an additional temperature measurement in the middle of the distillation column is easily explained. Since no setpoint is available for such a temperature measurement, an improved feedback may be calculated neither for the composition profile nor for the composition profile’s position.

6.5 Design of controllers with fixed structure

159

1

0.5

0 -5 10

RP

RS

-3

-1

10 10 Frequency (rad/min)

1

10

Structured singular value

Structured singular value

This special behavior of the control system has its significant advantages for the closed-loop behavior. The μ-plots in Figure 6.25 demon1

RP

0.5

0 -5 10

RS

-3

-1

10 10 Frequency (rad/min)

10

1

Figure 6.25: Robust performance and stability for PI control (left) and real PID control (right) with static decoupling

strate the superior robust performance of the decoupling control structures. While the optimal tuning for PI control shows a peak of the robust performance plot within the medium frequency range, the PID control structure shows nearly flat and significant smaller structured singular values. Using a decoupling control structure, the additional degrees of freedom in the controller design allow a significantly better controller performance, especially in the important mid-frequency range. Because of the better performance using PID controllers, the further discussion focuses on that control structure. The singular value plots of the loop transfers from the reference and disturbance signals to the output signals (Figure 6.26) illustrate the better controller performance as well. The condition numbers of T r → y are much smaller than those of the diagonal PI(D) control structure, and the tracking behavior is significantly improved. The simulation results (Figure 6.27) confirm the fundamentally improved controller performance. The sluggish behavior has vanished, and the maximum control errors are comparable to those obtained with

6 μ-Optimal Controller Design

160

1

Tr → y

1

10

Magnitude

Magnitude

10

0

10

-1

10

-2

10 -5 10

Td → y

0

10

-1

10

-2

-3

10 -5 10

-1

10

10

-3

-1

10

10

Frequency (rad/min)

Frequency (rad/min)

a)

b)

1

10

Figure 6.26: Singular values for the nominal closed-loop system for PID control with static decoupling a) Transfer function from reference to output signals b) Transfer functions from disturbance to output signals Dash-dotted line: T F → y , solid line: T x → y F

the diagonal PI(D) control structures. While the sensitivity to input uncertainty has increased, it is still small. PID Control with dynamic decoupling Using lead-lag transfer functions for the decoupler elements C 1 and C 2 1 + TDs C i ( s ) = KC --------------------1 + TLs

(6.38)

a dynamic decoupling structure is realized. The additional degrees of freedom allow a further improvement of the control design. The resulting optimal tuning constants are listed in Table 6.4. The simulation results exhibits a performance which is insignificantly worse than that of the μ-optimal state-space controllers (Figure 6.28). However, the more difficult initialization of a control structure with dynamic decoupling in a distributed control system is a disadvantage.

6.5 Design of controllers with fixed structure

0.015

0.010

0.005 0

0.020

Composition (mol/mol) Composition

Composition (mol/mol) Composition

0.020

Ft=0=20 mol/min F(t=0)=20 mol/min

10

161

Ft=0=46 mol/min F(t=0)=46 mol/min

0.015

0.010

Top composition

Top composition

Bottom composition

Bottom composition

20 30 Time [h] (h)

40

0.005 0

10

20 30 Time (h) [h]

40

Figure 6.27: Simulation results for PID control with static decoupling for an increase in feed composition (0.8 → 0.9 mol/mol) at t=0 h and an increase of feed flow rate (+ 3.6 mol/min) at t=20 h L, V equal controller output ΔL with +10% error, ΔV with –10% error

Table 6.4: μ-optimal parameters for PID control with dynamic decoupling Controller

KR (mol/min/°C)

KC (–)

TI (min)

TD (min)

TL (min)

PID 1

–22.2



80.2

19.6

44.8

PID 2

5.68



59.4

12.6

24.7

C1



–0.138



117

7.42

C2



1.07



53.0

71.43

Results achieved with numerically linearized model and complex μ-analysis

6 μ-Optimal Controller Design

162

0.015

0.010

0.005 0

0.020

Composition (mol/mol) Composition

Composition (mol/mol) Composition

0.020

Ft=0=20 mol/min F(t=0)=20 mol/min

10

Ft=0=46 mol/min F(t=0)=46 mol/min

0.015

0.010

Top composition

Top composition

Bottom composition

Bottom composition

20 30 Time (h) [h] Time

40

0.005 0

10

20 30 Time (h) [h] Time

40

Figure 6.28: Simulation results for PID control with dynamic decoupling for an increase in feed composition (0.8 → 0.9 mol/mol) at t=0 h and an increase of feed flow rate (+ 3.6 mol/min) at t=20 h L, V equal controller output ΔL with +10% error, ΔV with –10% error

6.5.3 PID control structures with one-way decoupling The results for two-way decoupling have shown optimal results for decoupling structures which are close to one-way decoupling. In this section the optimal tuning results for one-way decoupling are discussed. This control structure is particularly easy to implement in a distributed control system and simple to initialize. In order to keep the decoupler as simple as possible, the discussion is limited to static one-way decoupling. Two different decoupler structures are possible if we set either C 1 or C 2 of the control structure shown in Figure 6.24 to zero. While the results

6.5 Design of controllers with fixed structure

163

Table 6.5: μ-optimal parameters for PID control with static one-way decoupling Controller or KR decoupler No. (mol/min/°C)

TI (min)

TD (min)

TL (min)

C (–)

1

–10.5

45.7

2.18

5.01

0

2

5.35

67.4

13.4

13.9

1.05

Structured singular value

Results achieved with numerically linearized model and complex μ-analysis

1 RP

0.5

0 -5 10

Figure 6.29: Robust performance and stability for real PID-control with one-way decoupling

RS

-3

-1

10 10 Frequency (rad/min)

1

10

for the two-way decoupling lead us to expect a good performance for the first case ( C 1 = 0 ), no inference is possible for the second case ( C 2 = 0 ). In fact, the optimization results show insufficient performance for the second case ( C 2 = 0 ). Therefore a reversal of the decoupling control structure with shaping of the composition profile by the bottom composition controller and moving the composition profiles position by the top composition controller does not lead to results comparable to those obtained with the other decoupling structure. The tuning parameters for the controller with C 1 = 0 can be found in Table 6.5. The corresponding μ-curves (Figure 6.29) let us expect a

6 μ-Optimal Controller Design

164

0.015

0.010

0.005 0

0.020

Composition (mol/mol) Composition

Composition (mol/mol) Composition

0.020

Ft=0=20 mol/min F(t=0)=20 mol/min

10

Ft=0=46 mol/min F(t=0)=46 mol/min

0.015

0.010

Top composition

Top composition

Bottom composition

Bottom composition

20 30 Time [h] (h)

40

0.005 0

10

20 30 Time [h] (h)

40

Figure 6.30: Simulation results for PID control with static one-way decoupling for an increase in feed composition (0.8 → 0.9 mol/mol) at t=0 h and an increase of feed flow rate (+ 3.6 mol/min) at t=20 h L, V equal controller output ΔL with +10% error, ΔV with –10% error

performance somewhere between that of the PI control with static twoway decoupling and that of the real PID control with static decoupling. The simulation results in Figure 6.30 support this interpretation. Therefore this controller represents a structure which is simple and easily implemented in a distributed control system, distinguished by a sufficiently high controller performance.

6.6 Summary

165

6.6 Summary The comparison of the state-space controllers obtained by μ-synthesis with PID control structures obtained by μ-optimization leads to surprising results. The frequently heard opinion that state-space controllers are much superior to PID control structures apparently is not true for this distillation column. The PID control structures with decoupling exhibit nearly the same performance as that achieved with state-space controller of a higher order, provided that the PID control structures are optimally tuned. The visual results of the μ-curves and simulation plots shall be supported by numerical measures. For purposes of comparison, the integral square of the control errors t=40h



ISE =

2 ( t ) + e 2 ( t ) ]dt , [ e 10 44

(6.39)

0

and the integral of the time-multiplied absolute control errors t=40h

ITAE =



[ e 10 ( t ) + e 44 ( t ) ] t dt

(6.40)

0

have been calculated and summed up for both operating points and all controllers. While ISE punishes especially large control errors, the ITAE performance measure has a higher importance for the process industry because it punishes any undesirably sluggish disturbance rejection. Both criteria, relative to the result for the state-space controller using 3 temperature measurements, can be found in Table 6.6. The last two columns in this table state the maximum absolute value of the SSV (RP) and the value of the optimization criterion k

f(θ) =

∑ μΔ˜ { Fl [ P, K ( θ ) ] } 3

i=1

(6.41)

6 μ-Optimal Controller Design

166

Table 6.6: Comparison of controllers in time-domain Relative ISE

Relative ITAE

Max. μ ˜

State-space controller, 3 temp. measurements

1.00

1.00

0.85

1.00

State-space controller, 2 temp. measurements

0.82

1.05

1.04

1.54

Diagonal PI control

3.13

2.89

2.14

9.03

Diagonal PID control

2.08

1.87

1.53

3.98

PI control with static twoway decoupling

2.42

1.74

1.13

1.55

PID control with static two-way decoupling

1.44

1.23

0.91

1.19

PID control with dynamic two-way decoupling

1.18

1.12

0.88

1.12

PID control with static one-way decoup. (C1=0)

1.99

1.51

0.97

1.34

Control structure

Δ

Relative

∑ μ Δ3˜ ( jω )

relative to the value for the state-space controller using 3 temperature measurements. The high correlation of the ITAE and the optimization criterion are obvious. The single exception is the state-space controller using 2 temperature measurements, which may be caused by the convergence problems mentioned before. This table effectively illustrates the high performance achieved with simple and easily realized PID-control structures. The μ-optimization approach has proved to be an efficient tool for the optimal design of controllers with fixed structure.

6.7 References

167

6.7 References [6.1] Balas, G. J., J. C. Doyle, K. Glover, A. Packard, and R. Smith: μAnalysis and Synthesis Toolbox, MUSYN Inc., Minneapolis MN, and The MathWorks, Inc., Natick, MA (1991) [6.2] Balas, G. J., A. K. Packard, and J. T. Harduvel: “Application of μSynthesis Techniques to Momentum Management and Attitude Control of the Space Station,” Proc. 1991 AIAA Guidance, Navigation and Control Conference, New Orleans, LA (1991) [6.3] Chiang, R. Y., M. G. Safonov: Robust Control Toolbox User’s Guide, The Mathworks Inc., Natick, MA (1992) [6.4] Dailey, R. L.: “Lecture Notes for the Workshop on H∞ and μ Methods for Robust Control,” IEEE Conference on Decision and Control, Brighton (1991) [6.5] Doyle, J. C.: “Analysis of Feedback Systems with Structured Uncertainties,” IEE Proc., 129, Pt. D., No. 6, 242-250 (1982) [6.6] Doyle, J. C.: “Performance and Robustness Analysis for Structured Uncertainty,” Proc. of the 21st Conference on Decision and Control, (1982) [6.7] Doyle, J. C.: “Structured Uncertainty in Control System Design,” Proc. of the 24th Conference on Decision and Control, Ft. Lauderdale, FL (1985) [6.8] Doyle, J., K. Lenz, and A. Packard: “Design Examples Using μSynthesis: Space Shuttle Lateral Axis FCS During Reentry,” NATO ASI Series F: Computer and Systems Science, 34, 128-154 (1987)

6 μ-Optimal Controller Design

168

[6.9] Enns, D. F.: “Rocket Stabilization as a Structured Singular Value Synthesis Design Example,” Control Systems, 11, 4, 67-73 (1991) [6.10] Grace, A.: Optimization Toolbox — User’s Guide, The MathWorks, Inc., Natick, MA (1990) [6.11] Lin, J.-L., I. Postlethwaite, and D.-W. Gu: “μ-K Iteration: A New Algorithm for μ-synthesis,” Automatica, 29, 219-224 (1993) [6.12] Maciejowski, J. M.: Multivariable Feedback Design, AddisonWesley Publishing Company, Wokingham, England (1989) [6.13] McFarlane, D. C., and K. Glover: “Robust Controller Design Using Normalized Coprime Factor Plant Descriptions,” Lecture Notes in Control and Informations Science, 138, Springer-Verlag, Berlin (1990) [6.14] Packard, A., J. Doyle, and G. Balas: “Linear Multivariable Robust Control With a μ Perspective,” Trans. of the ASME, 115, 426-438 (1993) [6.15] Packard, A., and J. Doyle: “The Complex Structured Singular Value,” Automatica, 29 1, 71-109 (1993) [6.16] Shinskey, F. G., “Distillation control for Productivity and Energy Conservation,” 2nd ed., McGraw-Hill, New York, 194-203 (1984) [6.17] Skogestad, S.: “Dynamics and Control of Distillation Columns – A Critical Survey,” Preprints of the 3rd IFAC Symposium on Dynamics and Control of Chemical Reactors, Distillation Column and Batch Processes, April 26-29, 1992, College Park, MD, 1-25 (1992)

6.7 References

169

[6.18] Skogestad, S., and M. Morari: “Implications of Large RGA Elements on Control Performance,” Ind. Eng. Chem. Res., 26, 23232330 (1987) [6.19] Skogestad, S., and P. Lundström: “MU-Optimal LV-Control of Distillation Column,” Comp. Chem. Eng., 14, 4/5, 401-413 (1990)

170

6 μ-Optimal Controller Design

7.1 Introduction

171

Chapter 7 Controller Design for Unstructured Uncertainty — A Comparison

7.1 Introduction A controller design for the entire operating range of the distillation column (see Chapter 6) requires a structured uncertainty model incorporating two linear models, and a huge computational effort. Naturally, the question arises what controller performance and robustness properties can be achieved if we use simpler design methods, based on just one plant model for the nominal operating point (Model GN) and classical design methods or simple unstructured uncertainty bounds. A few of these simpler methods are discussed in this chapter. They are applied in a straightforward manner, and the design results are not guaranteed to represent the optimum achievable controller performance. However, the results give an impression of the limits and inherent problems of the application of design methods based on simpler uncertainty concepts, and they allow a comparison with the μ-optimal results presented in the previous chapter. The weighting functions of

172

7 Controller Design for Unstructured Uncertainty — A Comparison

the structured uncertainty model used for the μ-analysis are the same as those used in the previous chapter.

7.2 Diagonal PI-control A diagonal PI-control scheme seems to be most frequently used in conventionally controlled distillation columns. Usually these PI controllers are tuned on-line. Due to the large time constants of the composition dynamics, we cannot expect this on-line tuning approach to lead to a controller performance close to the optimum. The attempt to use tuning rules such as Ziegler-Nichols for the individual SISO loops often results in an unstable MIMO closed-loop system, because these tuning rules do not take the interaction between the two control loops into account. While the following two simple and model based tuning methods make use of the classical design methods, they try to pay attention to the loop interactions. Both methods lead to a nominally stable controller design. However, sufficient stability margins for the closed-loop system at all possible operating points cannot be guaranteed. 7.2.1 The BLT method The Biggest Log Modulus Tuning was proposed by Luyben in 1986 ([7.5], [7.6]). This method is a multivariable extension of the classical Nyquist stability criterion. The closed-loop system (Tr→y) with a square nominal model G ( s ) = G u → y ( s ) and a diagonal PI control law K(s) is given by y ( s ) = [ I + G ( s )K ( s ) ] – 1 G ( s )K ( s )r ( s )

(7.1)

The characteristic equation of the multivariable system is the scalar equation det ( I + G ( s )K ( s ) ) = 0

(7.2)

7.2 Diagonal PI-control

173

If we plot (7.2) as a function of frequency, the number of right half-plane zeros of the closed-loop characteristic equation are determined. In order to make this multivariable plot like the SISO scalar Nyquist plot, Luyben introduces a new function W(s): W ( s ) = – 1 + det ( I + G ( s )K ( s ) )

(7.3)

The closer this function approaches the (-1,0) point in the Nyquist plot, the closer the MIMO system is to closed-loop instability. The design objective is defined as W ( jω ) L lm = 20 log ------------------------1 + W ( jω )

≤ 2p ∀ω ∈ R +

(7.4)

where p is the number of inputs/outputs of G(s). The proposed tuning procedure starts with independent Ziegler-Nichols settings for PIcontrollers of the individual control loops. In a second step these settings are detuned by a factor F K ZN K i = -------------i F

TI i = F TI ZN

(7.5)

i

in order to achieve the design objective (7.4). Results of the BLT tuning The tuning results for the nominal model G N ( s ) of the distillation process are listed in Table 7.1. A detuning factor F of 3.82 was necessary to achieve the design objective (7.4). The proportional gain KR of the top composition controller is too large for satisfactory setpoint tracking and

Table 7.1: Tuning constants with BLT-method Controller

KR (mol/min/°C)

TI (min)

PI 1

–47.1

95.2

PI 2

6.74

171.8

174

7 Controller Design for Unstructured Uncertainty — A Comparison

Structured singular value

measurement noise attenuation. A plot of the structured singular values (with the same uncertainty and performance weights as used in the previous chapter) illustrates the insufficient robust stability and robust performance of this composition control design (Fig. 7.1). However, any further detuning would reduce the low and highfrequency gains of the bottom composition controller to an absolutely insufficient level. 3 RP

2

Figure 7.1: μ-plots for a diagonal PI-control law tuned with BLTmethod

1 RS

0 -5 10

-3

-1

10 10 Frequency (rad/min)

1

10

7.2.2 Sequential loop closing The idea of the sequential loop closing was introduced by Mayne ([7.8], [7.9]). First, a SISO controller is designed for one pair of input and output variables. When this design has been completed, the corresponding control loop is closed and the next pair of input and output variables is chosen. Thus the interaction between the control loops is taken into account. This design procedure is illustrated in Figure 7.2. It is an advantage of this method that each single loop can be designed using classical methods. However, this method has some severe drawbacks: First, the selection of the first one or two input/output pairs may have a deleterious effect on the behavior of the remaining loops [7.7]. There exists little help for this sequence problem. Second, this method cannot guarantee robustness for the entire operating range. Especially if the plant G(s) is not diagonal dominant, that means the condition

7.2 Diagonal PI-control

175

G3(s)

– K 1 (s)



G1(s)

y1

u1

y2

u2

G(s) .....

u3

y3 .....

K 2 (s)

G2(s)

Figure 7.2: Sequential loop closing

G ii ( jω ) > G ij ( jω ) ∀ω ∈ R +

(7.6)

is not satisfied, we have to expect robustness problems. Design Results The sequential loop closing idea has been applied to the composition control problem represented by the nominal model G N ( s ) . For each SISO loop, a phase margin of at least 60 degrees and for both controllers a maximum high frequency gain of 18 mol/min/°C has been required. The results of both possible design sequences and with a minimal integral absolute error (IAE) for the rejection of feed composition and feed flow disturbances (with respect to the linear model) are summarized in Table 7.2. An analysis of the robustness to unstructured peturbations shows maximum values for the sensitivities of Se=2.6, and Su=2.1 for the

176

7 Controller Design for Unstructured Uncertainty — A Comparison

Table 7.2: Results of the sequential loop closing Design sequence

KR 1 (mol/min/ °C)

TI 1 (min)

KR 2 (mol/min/ °C)

TI 2 (min)

Top → Bottom

–18.0

101.9

10.09

55.8

Bottom → Top

–18.0

52.6

8.78

214.5

3

2

1

RP RS

0 -5 10

10

-3

-1

10

1

10

Structured singular value

Structured singular value

Top → Bottom design sequence, and of Se=2.4 and Su=1.9 for the sequence Bottom → Top. These stability margins are insufficient. The results of the analysis using the structured uncertainty model are illustrated by the μ-plots in Figure 7.3. Both controller designs can neither guarantee robust performance nor robust stability. 3

RP 2

1

RS 0 -5 10

-3

-1

10

10

Frequency (rad/min)

Frequency (rad/min)

a)

b)

10

1

Figure 7.3: μ-plots for the sequential loop closing designs a) Top → Bottom design sequence b) Bottom → Top design sequence

7.2.3 Optimized robust diagonal PI-control The objective of this controller design is a maximization of the disturbance rejection capabilities with the boundary conditions of sufficient stability margins. As a measure of the disturbance rejection capabilities the IAE as defined by

7.2 Diagonal PI-control

177

T End

IAE =



[ e 10 ( t ) + e 44 ( t ) ]dt

(7.7)

0

is a suitable measure. It is calculated for step responses to feed composition and feed flow rate of the closed-loop system. If we tune both PIcontrollers in order to minimize the IAE-criterion, the robustness properties of the closed-loop system form boundary conditions for the minimum achievable IAE. Stability bounds in terms of the sensitivity at the plant input and output are well established. If we require a phase margin of at least 35 degrees (which is relatively small), the following sensitivity bounds hold S e ( jω ) = [ I + G ( jω )K ( jω ) ] – 1 < 1.7

∀ω ∈ R +

(7.8)

S u ( jω ) = [ I + K ( jω )G ( jω ) ] –1 < 1.7

∀ω ∈ R +

(7.9)

The optimal parameters which minimize the IAE criterion are found either by trial and error or by a constrained parameter optimization. Results The results for this design approach are given in Table 7.3. The corresponding μ-plots (Fig. 7.4) illustrate the improved robust stability properties compared to the previous two methods. While design guarantees robust stability, the robust performance is substantially worse than the

Table 7.3: Tuning constants with optimizing method Controller

KR (mol/min/°C)

TI (min)

PI 1

–5.10

600.0

PI 2

4.92

86.2

Structured singular value

178

7 Controller Design for Unstructured Uncertainty — A Comparison

8 6 Figure 7.4: μ-plots for diagonal PI-controller designed by optimizing method

RP

4 2 RS

0 -5 10

-3

-1

10 10 Frequency (rad/min)

1

10

μ-optimal design of a diagonal PI-controller design (see Figure 6.18, page 151). An analysis of the controller behavior in the time domain (Figure 7.5) shows extremely sluggish disturbance rejections.

7.3 PI-control with decoupling The basic idea of decoupling is a reduction of the loop interactions. If we increase the diagonal dominance of the system, the design task takes on more the characteristics of a multiloop SISO design problem. However, as emphasized already in the previous chapter, a reduction of the loop interactions does not automatically imply better control. Due to an increased sensitivity to model and input errors, the maximum performance of a controller exhibiting sufficient stability margins may be strongly reduced even compared to that of a diagonal PI controller. In the simplest case, as discussed in this section, the plant behavior is altered by a pre- or postmultiplied constant “compensating” or “decoupling” matrix. Different approaches for the selection of these interaction reducing matrices are proposed: Davison [7.3] recommends a steady-state decoupling of the process. For the “decoupled” process G * ( s ) holds

7.3 PI-control with decoupling

Ft=0=20 mol/min F(t=0)=20 mol/min

0.015

Top composition Bottom composition

0.010

0.005 0

10

20 30 [h] Time (h)

0.020

Composition (mol/mol) Composition

Composition (mol/mol) Composition

0.020

179

40

Ft=0=46 mol/min F(t=0)=46 mol/min

0.015

Top composition Bottom composition

0.010

0.005 0

10

20 30 Time (h) [h]

40

Figure 7.5: Simulation results with diagonal PI controller (designed by optimizing method) for an increase in feed composition (0.8 → 0.9 mol/mol) at t=0 h and an increase of feed flow rate (+ 3.6 mol/min) at t=20 h Upper plots: Product composition Lower plots: Control error L, V equal to controller output ΔL with +10% error, ΔV with –10% error

G * u → y ( s ) = G u → y ( s )G –1 u → y ( 0 ) or

(7.10)

G * u → y ( s ) = G – 1 u → y ( 0 )G u → y ( s ) With a state space representation of the process, the decoupling matrix is calculated according to G – 1 u → y ( 0 ) = ( CA –1 B ) –1

(7.11)

180

7 Controller Design for Unstructured Uncertainty — A Comparison

The choice of a premultiplication or postmultiplication of this interaction reducing matrix is another degree of freedom for the controller design. Mayne [7.9] proposes a reduction of the high-frequency interactions of the plant. The corresponding decoupling matrix is calculated by G – 1 u → y ( j∞ ) = ( CB ) – 1

(7.12)

As before, the choice of a pre- or postcompensation has to be decided during the controller design. Ryskamp [7.11] suggests a decoupling scheme which is based on the idea of a composition profile control: The difference in the temperature deviations should be used to set the reflux ratio, and the sum of the two temperature deviations should be used to set the reboiler heat duty. This scheme is called “implicit decoupling.” Another interesting approach, based on a singular value decomposition (SVD) of the process at steady state, is presented by Brambilla et al. [7.1]. Let the SVD of the steady-state transfer matrix of the process G u → y ( 0 ) be G u → y ( 0 ) = UΣV T ,

(7.13)

where U and V are unitary matrices and Σ is a diagonal matrix containing the singular values Σ = diag ( σ 1, σ 2 ) . A plant-inverting ˆ (at plant input) according to this SVD is the matrix compensator D ˆ = VΣ –1 U T D

(7.14)

In order to avoid a high sensitivity to input errors due to the perfect decoupling at steady state, Brambilla et al. [7.1] introduce a matrix F F = αI + ( 1 – α )Σ

(7.15)

˜ as and define a new compensation matrix D ˜ = VFΣ – 1 U T D

(7.16)

7.3 PI-control with decoupling

181

The single parameter α with α = 0…1 allows a continuous shift between a plant-inverting compensator ( α = 1 ) and a compensator which does not remove the effect of the directionality of the process ( α = 0 ). The tuning parameter α has to be chosen on the basis of (1) the magnitude of the assumed errors in the model, (2) the sensitivity of the ˜ ), and (3) the process to the model errors (Relative Gain Analysis of D required performance in terms of reduction of interactions and direc˜ ). tionality (Relative Gain Analysis of G u– 1→ y ( 0 )D Design results The four proposed compensation matrices are summarized in Table 7.4. In order to calculate “optimal” controllers, the optimization approach described in section 7.2.3 has been applied to the different compensated plants. However, it was not possible to achieve any acceptable controllers using the proposed compensation matrices, except for the SVDbased compensator. This SVD-based compensator is distinguished by almost the same one-way decoupling structure as we obtained as a result of the μ-optimal decoupling (see Chapter 6). Table 7.4: Compensator matrices Type of compensator Position of compensator Decoupling at

Compensator matrix

Plant input or plant output

– 0.636 0.168 – 0.728 0.195

ω = ∞

Plant input or plant output

0.380 – 0.295 0.875 – 0.193

Implicit decoupling

Plant output

–1 –1 –1 1

SVD-based compensation ( α = 0.8 )

Plant input

0.901 0.082 0.955 0.391

ω = 0

Decoupling at

The parameters of the IAE-optimal PI controllers (with respect to an additional boundary condition for the proportional gains KR i < 18 mol/

Structured singular value

182

7 Controller Design for Unstructured Uncertainty — A Comparison

RP

1 Figure 7.6: μ-plots for SVD-based compensation with optimally tuned diagonal PI control.

RS

0.5

0 -5 10

-3

-1

10 10 Frequency (rad/min)

1

10

min/°C) are given in Table 7.5. The μ-plots for this controller design Table 7.5: Optimal PI tuning constants for plant with SVD-based compensation Controller

KR (mol/min/°C)

TI (min)

PI 1

–18.0

47.4

PI 2

18.0

116.0

(Figure 7.6) demonstrate good robust performance and robust stability. However, the simulation results (Figure 7.7) show an insufficiently damped oscillation at higher frequencies for the minimum feed flow rate. The damping of these oscillations is significantly better for a +10% error in the change of the reflux L and a –10% error in the change of boilup V. These unwanted oscillations are allowed by the performance specification in the frequency domain! They require a detuning of the controllers’ proportional gains which on the other hand, reduces the controller performance.

7.3 PI-control with decoupling

Ft=0=20 mol/min F(t=0)=20 mol/min

0.020

Composition (mol/mol) Composition

Composition (mol/mol) Composition

0.020

183

0.015

0.010

Ft=0=46 mol/min F(t=0)=46 mol/min

0.015

0.010

Top composition

Top composition

Bottom composition

0.005 0

10

20 30 Time [h] Time (h)

Bottom composition

40

0.005 0

10

20 30 Time [h] Time (h)

40

Figure 7.7: Simulation results with SVD based compensator and diagonal PI control for an increase in feed composition (0.8 → 0.9 mol/mol) at t=0 h and an increase of feed flow rate (+ 3.6 mol/min) at t=20 h Upper plots: Product composition Lower plots: Control error L, V equal to controller output ΔL with +10% error, ΔV with –10% error

184

7 Controller Design for Unstructured Uncertainty — A Comparison

7.4 H∞ optimal design The H∞-norm minimizing design ([7.2], [7.4], [7.7], [7.10]) of multivariable controllers have proved to be a powerful method for robust, modelbased controllers. H∞ Design specification The closed-loop system with the plant G(s) and the controller K(s), augmented with the weighting functions W d ( s ) , W e ( s ) , W u ( s ) , and W y ( s ) is outlined in Figure 7.8. This scheme is often called S/KS/Tweighting scheme. The matrix W d ( s ) is a diagonal matrix of transfer functions and represents the frequency content of the feed composition, feed flow rate, and reference input signals. The selection of these input weights is discussed in section 6.3. The same weighting functions are applied here.

We(s)

ze(s)

Wu(s)

zu(s)

Wy(s)

zy(s)

d(s) r(s)

Wd(s)

G(s)

+

e(s)

K(s)

y(s)

u(s)



Figure 7.8: Augmented closed-loop system with weighting functions for the H∞ design

H∞ optimal design

185

All other weighting functions are chosen as diagonal frequency-dependent weights because the performance and robustness properties are equal for all channels: W e ( s ) = diag [ w e ( s ), w e ( s ) ]

(7.17)

W u ( s ) = diag [ w u ( s ), w u ( s ) ]

(7.18)

W y ( s ) = diag [ w y ( s ), w y ( s ) ]

(7.19)

The performance of the closed-loop system is specified in terms of the sensitivity function by the weighting function W e ( s ) . A first-order lag with a static gain of 100 has been specified to achieve a nearly integrating behavior. The bandwidth of the closed-loop system is limited by the weighting function W y ( s ) , which punishes the transfer function T [ d T, r T ] T → y from the disturbance and reference signals to the plant outputs. A firstorder lead-lag transfer function is suitable for this task. A weighting of the plant inputs allows a frequency-dependent limitation of the control energy and helps to achieve sufficient stability margins for the sensitivity function at u. As done with W y ( s ) , a first-order lead-lag transfer function has been selected. The poles and zeros of the weighting functions were adjusted until the sensitivity functions at e and at u of the closed-loop system had attained approximately the same peak values as the μ-optimal controller design (with 2 temperature measurements), a high performance, and T

d →z r

( jω )

≈1

(7.20)



were achieved. The best weighting functions are given by 1 w e ( s ) = 100 -----------------------1 + 9520s

(7.21)

186

7 Controller Design for Unstructured Uncertainty — A Comparison

1 + 520s w u ( s ) = 0.5 --------------------1 + 13s

(7.22)

1 + 1500s w y ( s ) = 0.1 -----------------------1 + 1.5s

(7.23)

Design results

Structured singular value

Despite the fact that the singular values of sensitivity functions for the H∞- design (Figure 7.10) and for the μ-synthesis (Figure 6.15) are nearly identical, the μ-analysis shows significant differences. The μ-plots of the H∞ design (Figure 7.9) show much higher peak values in the low and mid-frequency ranges. The simulation results (Figure 7.11) allow a conclusion with respect to the larger structured singular values: The sensitivity of the closed-loop performance to errors in the manipulated variables is large. A reduction of this sensitivity to plant input errors was not possible using the common S/KS/T weighting scheme.

1.5

1

RP Figure 7.9: μ-plots for H∞ optimal controller

RS

0.5

0 -5 10

-3

-1

10 10 Frequency (rad/min)

1

10

H∞ optimal design

187

Sensitivity at e

1

10

0

Magnitude

10

-1

10

-2

10

-3

10 -5 10

-4

10

-3

10

-2

10

10

-1

0

10

10

1

Frequency (rad/min) Sensitivity at u

1

10

0

Magnitude

10

-1

10

-2

10

-3

10 -5 10

-4

10

-3

10

-2

10

10

-1

0

10

10

1

Frequency (rad/min) Figure 7.10: Singular values of the sensitivity functions at e (upper plot) and at u (lower plot) for the nominal closed-loop system with the H∞ controller

188

7 Controller Design for Unstructured Uncertainty — A Comparison

0.015

0.010

0.005 0

0.020

Composition (mol/mol) Composition

Composition (mol/mol) Composition

0.020

Ft=0=20 mol/min F(t=0)=20 mol/min

10

Ft=0=46 mol/min F(t=0)=46 mol/min

0.015

0.010

Top composition

Top composition

Bottom composition

Bottom composition

20 Time [h] (h)

30

40

0.005 0

10

20

30

40

Time Time (h) [h]

Figure 7.11: Simulation results with the H∞-controller for an increase in feed composition (0.8 → 0.9 mol/mol) at t=0 h and an increase of feed flow rate (+ 3.6 mol/min) at t=20 h Upper plots: Product composition Lower plots: Control error L, V equal to controller output ΔL with +10% error, ΔV with –10% error

7.5 Summary

189

7.5 Summary The application of design methods for unstructured uncertainty to the composition (or temperature) control problem shows that it is extraordinarily difficult to obtain performances which are comparable to those of the μ-optimal controllers. Despite the high effort for a robust tuning of the PI-control structures, it was not possible to achieve any satisfactory result. Better results were obtained using the H ∞ -minimization approach. The resulting state-space controller guarantees stability for the entire operating range and the singular values of the sensitivity functions (Se, Su) are nearly identical to those of the μ-optimal state-space controller. Nevertheless, the high sensitivity to input uncertainty demonstrates the limits of simple unstructured uncertainty bounds. Even a robust controller design based on an unstructured uncertainty model tends to be very sensitive to input uncertainty at operating points different from the design point. The advantages of a μ-optimal controller design as presented in Chapter 6 are obvious.

7.6 References [7.1] Brambilla, A., and L. D’Elia: “Multivariable Controller for Distillation Column in the Presence of Strong Directionality and Model Errors,” Ind. Eng. Chem. Res., 31, 536-543 (1992) [7.2] Dailey, R. L.: “Lecture Notes for the Workshop on H∞ and μ Methods for Robust Control,” 1991 IEEE Conference on Decision and Control, Brighton, December 9-10 (1991) [7.3] Davison, E. J.: “Multivariable tuning regulators: The feedforward and robust control of general servomechanism problems,” IEEE Trans. Aut. Control, AC-21, 35-47 (1976) [7.4] Glover, K., and J. C. Doyle: “A State Space Approach to H∞ Optimal Control,” Lecture Notes in Control and Information Sciences, 135, 179-218, Springer-Verlag, Berlin (1989)

190

7 Controller Design for Unstructured Uncertainty — A Comparison

[7.5] Luyben, W. L.: “Simple Method for Tuning SISO Controllers in Multivariable Systems,” Ind. Eng. Chem. Process Des. Dev., 25, 654-660 (1986) [7.6] Luyben, W. L.: Process Modeling, Simulation, and Control for Chemical Engineers, 2nd ed., McGraw-Hill, New York (1990) [7.7] Maciejowski, J. M.: Multivariable Feedback Design, AddisonWesley Publishing Company, Wokingham (1989) [7.8] Mayne, D. Q.: “The design of linear multivariable systems,” Automatica, 9, 201-207 (1973) [7.9] Mayne, D. Q.: “Sequential design of linear multivariable systems,” Proc. IEE., 126, 6, 568-572 (1979) [7.10] Raisch, J., L. Lang, und E.-D. Gilles: “H∞-Reglerentwurf für Zwei- und Dreistoffdestillationsprozesse”, at, 41, 6, 215-224 (1993) [7.11] Ryskamp, C. J.: “Explicit vs. implicit decoupling in distillation control,” Chemical Process Control II, American Institute of Chemical Engineers, New York, 361-375 (1982)

8.1 Introduction

191

Chapter 8 Feedforward Controller Design

8.1 Introduction It is a drawback of feedback control that a corrective action necessitates a deviation of the controlled variables from their setpoints. This disadvantage can be overcome by the use of feedforward control. A major and probably the most frequent disturbance of a distillation column is a change in the feed flow rate. Because the feed flow rate is always measured, it can be used as a controller input. An appropriately designed feedforward controller takes most of the necessary corrective action before the product compositions and the controlled tray temperatures change. However, because of model errors and other unmeasured disturbances a feedforward controller alone will never be able to yield perfect control so that feedback control will still be needed. Within this chapter, the design of linear time-invariant feedforward controllers for our distillation column is discussed. The proposed design methods take into account the wide operating range of the distillation column and the unmeasured feed composition.

192

8 Feedforward Controller Design

8.2 The design problem 8.2.1 The design objective The objective of feedforward control is a reduction of the control error in presence of feed flow rate disturbances. The main problem is the nonlinear behavior of distillation columns. The perfect control action for a rejection of a feed flow disturbance depends on the actual and measured feed flow rate and the unmeasured feed composition. A controller design for one operating point may be unsatisfactory at any others. Consequently, it is impossible to design a perfect linear timeinvariant feedforward controller for the entire operating range of a distillation column. Hence the design objective is a feedforward controller which improves the compensation of feed flow disturbances for the largest possible part of the operating range, but never makes it worse. A perfect solution of this design objective would be an enormous task. Fortunately, the ideas discussed in Chapter 5 lead to very good results: If we design the feedforward controller simultaneously for the models G R ( s ) (representing minimum feed flow rate and maximum feed composition) and G I ( s ) (representing maximum feed flow rate and minimum feed composition), we obtain a design which improves the compensation of feed flow disturbances for the entire operating range. 8.2.2 One-step or two-step design? The design of feedforward controllers is feasible either in a one-step design, simultaneously with the feedback controller, or as a second step for the closed-loop system (Fig. 8.1) [8.3]. The design of a feedforward controller for the open-loop system is not recommended because the feedback controller shifts the poles and, consequently, affects the dynamics of the system. A μ-optimal one-step design using the uncertainty structure presented in chapter 5 is tempting. For that purpose, the uncertainty structure is slightly modified by the additional input to the controller, i.e., the

8.2 The design problem

193

a) One-step design

Δ

xF F r

p

P wF(s)

KF(s) K(s)

b) Two-step design Δ Step 1: Feedback design d r

P

p

K(s)

Δ

Step 2: Feedforward design F r

P* wF(s)

KF(s)

+

p

+

K(s)

Figure 8.1: Design of feedforward controllers a) Simultaneous design with feedback controller b) Design as a second step for the closed loop system. The weighted plant P* may be a simpler uncertainty structure than the plant P.

194

8 Feedforward Controller Design

weighted feed flow signal. However, this approach has certain drawbacks: • Convergence is unattainable using the μK-Iteration in our case • Using the μ-optimization approach, the one-step design needs significantly more computation time than the two-step design. (The computing cost is proportional to {number of parameters}k with k>2.) • For acceptable results, the weighting function for the feed flow signals must be modified: Small improvements in the compensation of feed flow disturbances cause a dominance of the reference and the feed composition inputs with regard to the performance specification. Very small gains in the feedforward part result therefrom. Consequently, the discussion is focused on the design of feedforward controllers for the closed-loop system, i.e., as a second design step. Since feedforward control does not affect any stability properties of the closedloop system, the design is relatively simple. It is discussed by means of two examples.

8.3 H∞-minimization The H∞-minimization [8.4] is well suited for a feedforward controller design. Before we use the numerical tools available (e.g., [8.1], [8.2]), we have to build up a closed-loop plant with a previously designed feedback controller K(s). As an example, the μ-optimal state-space controller using all 3 temperature measurements is selected (see section 6.4.3). If we wish to improve the compensation of feed flow disturbances for the plant models G R ( s ) as well as for G I ( s ) , we have to close the feedback loops for both models separately, define the desired performance, and limit the high-frequency output of the feedforward controller K F ( s ) . The design plant is outlined in Figure 8.2. The performance weight W e ( s ) is a diagonal matrix of the transfer functions w e ( s )

H∞-minimization

195

Wu ( s )

F u

KF ( s )

uF

GR ( s ) –

+ F

K (s)

+ We ( s )

F u

+ +

zu

ze

GI ( s ) – K (s)

Figure 8.2: The augmented plant for a design of the feedforward controller KF(s) by H∞-minimization

W e ( s ) = diag [ w e ( s ), w e ( s ), w e ( s ), w e ( s ) ]

(8.1)

It demands the same performance for both column models and both controlled temperatures. The transfer function w e ( s ) is chosen as a first-order lag with a high static gain. The pole of w e ( s ) is adjusted until T F → z ∞ ≈ 1 is achieved. The final transfer function becomes 1 w e ( s ) = 100 -----------------------1 + 2380s

(8.2)

If we do not specify any high-frequency limits of the feedforward controller output, we obtain a controller with large high-frequency gains. This is undesirable because measurement noise and short-time feed flow fluctuations cause unnecessarily, large control actions. Using a diagonal transfer function matrix W u ( s ) for the feedforward controller output u F according to

196

8 Feedforward Controller Design

W u ( s ) = diag [ w u ( s ), w u ( s ) ] F

(8.3)

F

with the lead-lag transfer function 1 + 104s w u ( s ) = 0.5 --------------------F 1 + 2.5s

(8.4)

a controller behavior similar to a first order lag is obtained. The singular values of the controller and the transfer functions from the disturbance inputs to the control error (for the nominal model) are shown by Figure 8.3. If we compare Figure 8.3 b with Figure 6.10 b, we recognize the significant improvement of the feed flow disturbance compensation (dash-dotted lines). KF

1

0

10

-1

10

-2

10 -5 10

Td → y

1

10

Magnitude

Magnitude

10

0

10

-1

10

-2

-3

-1

10

10

1

10

10 -5 10

10

-3

-1

10

Frequency (rad/min)

Frequency (rad/min)

a)

b)

1

10

Figure 8.3: a) Singular values of the feedforward controller b) Singular values of the transfer function from the disturbance inputs d to the controlled output signals y for the nominal model GN with feedforward and feedback control. Solid line: T x → y , dash-dotted line: T F → y F

Nonlinear simulations confirm these expectations (Figure 8.4). In the interest of consistency, the same disturbances are simulated as in all previous chapters. Of course, the response to the step changes in the feed composition remains identical to the one shown in Figure 6.12.

H∞-minimization

Ft=0=20 mol/min F(t=0)=20 mol/min

0.010

0.4 Temperature Temperature(K)

Composition Composition (mol/mol)

0.015

0.005 0

Top composition Bottom composition

10

20 30 Time [h] (h)

0.020

40

Ft=0=20 mol/min F(t=0)=20 mol/min

0.2 0.0 -0.2 -0.4 -0.6 0

Control error T-44

20 30 Time [h] (h)

0.015

0.010

0.4

Top composition Bottom composition

10

40

20 30 Time [h] (h)

40

Ft=0=46 mol/min F(t=0)=46 mol/min

0.2 0.0 -0.2 -0.4

Control error T-10

10

Ft=0=46 mol/min F(t=0)=46 mol/min

0.005 0

Temperature Temperature(K)

Composition Composition (mol/mol)

0.020

197

-0.6 0

Control error T-10 Control error T-44

10

20 30 Time [h] (h)

40

Figure 8.4: Simulation results with μ-optimal state space controller (controller inputs: T10, T44, T24) and feedforward controller for an increase in feed composition (0.8 → 0.9 mol/mol) at t=0 h and an increase of feed flow rate (+ 3.6 mol/ min) at t=20 h Upper plots: Product composition Lower plots: Control error L, V equal to controller output ΔL with +10% error, ΔV with –10% error

198

8 Feedforward Controller Design

However, the maximum control error during the compensation of the feed flow disturbance is approximately halved, and for the maximum feed flow rate it is reduced even more.

8.4 Optimization approach The implementation of state-space controllers in a distributed control system is difficult. Of course, this holds for feedforward controllers as well. Most desirable are feedforward controllers with a simple and easily implementable structure. The singular values of the H∞ norm-minimizing state-space controller suggest a feedforward controller structure with a first-order lag and different gains for the outputs to the reflux L and the boilup V:

KF ( s ) =

KR L KR V

1 ---------------1 + Ts

(8.5)

The parameters of this simple control structure are computed by a constraint parameter optimization [8.5]. The objective may be of different kind: One possibility is the minimization of the H∞ norm of the transfer function T F → z for the plant shown in Figure 8.2. This design objective has certain disadvantages, however: • Due to the few degrees of freedom resulting from using this simple controller structure, it is not possible to obtain a controller which is close to design specifications for a wide frequency range. • The H∞-norm minimizing parameters strongly depend on the allowed maximum for γ , with γ = T F → z ∞ . If we allow γ > 5 controller designs with large enough gains ( KR L , KR V ) are obtained. But for a performance specification allowing γ ≈ 1 , we attain small controller gains and the improvement of the disturbance compensation is insufficient

8.4 Optimization approach

199

Most of the feed flow disturbances entering this distillation column are step changes. Consequently, we are able to define an appropriate design objective in the time domain. It is the minimum absolute control error for a step change in the feed flow rate. The design objective becomes [ T, KR L, KR V ] = arg

inf

[ T, KR L, KR V ]

E

(8.6)

with Te

E =

∫ { e10R ( t )

+ e 44 ( t ) + e 10 ( t ) + e 44 ( t ) }dt . R

I

(8.7)

I

0

The performance measure E is calculated for a step response to the plant input F, employing the plant illustrated by Figure 8.5. F u



+ F

KF ( s )

uF

GR ( s )

K (s)

+

e 10

R

e 44

R

F u

+ +

GI ( s ) – K (s)

e 10

I

e 44

I

Figure 8.5: Plant structure for the optimization of feedforward controller parameters

If we select the μ-optimal PID-controller with one-way decoupling as the feedback controller K (see section 6.5.3), and limit the time constant T

200

8 Feedforward Controller Design

by a lower bound of 5 minutes, the following simple optimal feedforward controller results: 1 K F ( s ) = 1.5 -------------------2.6 1 + 5.0s

(8.8)

The singular values of the feedforward controller are shown in Figure 8.6 a. In Figure 8.6 b we find the singular values of the transfer functions T d → y for nominal closed loop system with this feedforward controller. It demonstrates the low sensitivity of the feedback and feedforward controlled distillation column to variations of the feed flow rate. KF

1

10

Magnitude

Magnitude

10

0

10

-1

10

-2

10 -5 10

10

10

1

Td → y

0

-1

-2

-3

-1

10

10

10

1

10 -5 10

Frequency (rad/min)

10

-3

-1

10

1

10

Frequency (rad/min)

a)

b)

Figure 8.6: a) Singular values of the feedforward controller with fixed structure b) Singular values of the transfer functions for the nominal closed loop system from the disturbances inputs d to the controlled output signals y (Feedback and feedforward control). Solid line: T x → y , dash-dotted line: T F → y F

The simulation results (Figure 8.7) demonstrate that the maximum deviation of the product compositions for a step change of the feed flow rate is very small. A comparison with the simulation results for the same feedback controller but without feedforward control in Figure 6.30 confirms the substantial improvements by this simple feedforward controller.

8.5 Summary

Ft=0=20 mol/min F(t=0)=20 mol/min

0.015

0.010

0.005 0

0.020

Composition (mol/mol) Composition

Composition (mol/mol) Composition

0.020

201

10

Ft=0=46 mol/min F(t=0)=46 mol/min

0.015

0.010

Top composition

Top composition

Bottom composition

Bottom composition

20 30 Time (h) [h] Time

40

0.005 0

10

20 30 Time [h] (h)

40

Figure 8.7: Simulation results with μ-optimal PID controller with one-way decoupling and a simple feedforward controller for an increase in feed composition (0.8 → 0.9 mol/mol) at t=0 h and an increase of feed flow rate (+ 3.6 mol/min) at t=20 h L, V equal to controller output ΔL with +10% error, ΔV with –10% error

8.5 Summary The compensation of feed flow disturbances can be improved by using feedforward controllers. H∞-norm minimization and the minimization of the control errors in the time domain (for feedforward controllers with fixed structure) are efficient design methods. Frequency domain as well as time-domain results demonstrate the pleasing improvements which are obtained even by a feedforward controller of order one. A comparison of the ISE and ITAE criteria (see section 6.6) in Table 8.1 demonstrates

202

8 Feedforward Controller Design

improvements up to 50%! As mentioned previously, the maximum structured singular value μ is not a good performance measure if we include the feedforward control in the structured uncertainty model. Table 8.1: Comparison of controllers in time-domain Relative ISE

Relative ITAE

Max μ ˜

State-space controller, 3 temp. measurements

1.0

1.0

0.85

State-space controller, 3 temp. measurements and feed forward control

0.63

0.51

0.86

PID control with static one-way decoupling (C 1=0)

1.99

1.51

0.97

PID control with static one-way decoupling (C 1=0) and simple feedforward control

1.13

0.87

1.05

Control structure

Δ

8.6 References [8.1] Balas, G. J., J. C. Doyle, K. Glover, A. Packard, and R. Smith: μAnalysis and Synthesis Toolbox, MUSYN Inc., Minneapolis MN, and The MathWorks, Inc., Natick, MA (1991) [8.2] Chiang, R. Y., and M. G. Safonov: Robust Control Toolbox — User’s Guide, The MathWorks, Inc., Cochituate Place, Natick, MA (1992) [8.3] Christen, U., M. F. Weilenmann, and H. P. Geering: “Design of H2 and H∞ Controllers with Two Degrees of Freedom,” Proc. of the 1994 American Control Conference, Baltimore, MA (1994)

8.6 References

203

[8.4] Glover, K., and J. C. Doyle: “A State Space Approach to H∞ Optimal Control,” Lecture Notes in Control and Information Science, 135, 179-218, Springer-Verlag, Berlin (1989) [8.5] Grace, A.: Optimization Toolbox — User’s Guide, The Mathworks, Inc., Natick, MA (1990)

204

8 Feedforward Controller Design

9.1 Introduction

205

Chapter 9 Practical Experiences

9.1 Introduction In simulations the performance of controllers is tested in a sterile environment. Lacking measurement noise, operator actions, and varying environmental conditions, the results of these simulations represent a well established basis for a comparison of different controller designs. However, only the implementation of a controller in the real plant proves its performance. While in the literature a great number of design methods has been proposed and the resulting controllers have been tested by simulations, only very few results of an implementation at a real industrial distillation column have been reported. This chapter complements the simulation results presented in previous chapters with the results of a controller implementation in the distributed control system (DCS) which is coupled with this distillation column. The first section describes the implementation including the handling of constraints. Further sections discuss the use of pressure compensated temperatures, the controller performance observed, and economic aspects. A short summary concludes the chapter.

206

9 Practical Experiences

9.2 Controller implementation In the research field the objective of any control design is a high controller performance. A control design implemented in an industrial environment must consider many additional aspects. A few of them are listed below. Simple implementation: As mentioned previously, state-space controllers are difficult to implement in a DCS. Therefore the control scheme should be based on fixed low order structures, e.g. on PID control or on advanced PID control structures. Robustness: The control design must guarantee stability for the entire operating range of the column, including time variations due to corrosion of trays, transmitter drifts, etc. Easy to initialize: The switch from manual to automatic control must be simple and easy to understand. Operators often are semiskilled workmen who cannot and should not be expected to have an engineering background. A complex initialization procedure of a control scheme unnecessarily increases the risk of errors and requires an intensive operator training. Handling of constraints: Constraints are necessary to prevent the column from flooding, weeping, overpressure, overtemperature, etc. Often it is sufficient to limit reflux and reboiler heat duty. Performance: Despite the requirements listed above, the performance of the control scheme should still be high. Comparing the different control schemes proposed within this thesis, the PID control structure with one-way decoupling including the simple feedforward controller evolves as the best compromise among all these requirements. The control scheme is simple to initialize1, robust to plant 1. Initialization of the control scheme (see Figure 9.1): First, the output of the top composition controller in manual mode is adjusted to achieve rR ≅ Lactual. Then the top composition controller is switched to automatic mode. Second, the output of the bottom composition controller in manual mode is adjusted to achieve rQ ≅ Qactual. After that the controller is switched to automatic mode.

9.2 Controller implementation

207

uncertainty, it allows a simple handling of constraints, and it exhibits a high performance in simulations. This control scheme has been implemented in the DCS installed at the plant considered here (i.e., an Eckardt PLS 80E). The controller inputs are estimated tray compositions xˆ i , which for the operators have proved to be easier to understand than pressure compensated temperatures. The proportional gains of the controllers are easily converted for these controller inputs. A scheme of the implementation is shown in Fig. 9.1. The handling of constraints is realized by using the anti-windup facility of the standard PID controller blocks within the DCS. The following ideas have been realized: • If the setpoint for the reflux controller r R becomes smaller than its minimum limit R min , the top composition is allowed to rise above the setpoint (⇒ top composition purer than required), and the top composition controller must be prevented from windup. • If the setpoint for the reflux controller r R exceeds its maximum limit R max , the top composition is allowed to decrease (⇒ top composition less pure than required), and the top composition controller again must be prevented from windup. However, if at all possible, this case should be avoided. • Equivalent constraints hold for the bottom composition controller. This policy establishes individual constraints for the top composition control loop as well as for the bottom composition control loop. Since we have to include the feedforward control and the one-way decoupling, the outputs of the composition controllers are limited by the following four signals entering the anti-windup facility of the PID controllers: R FB1, max = R max – R FF

(9.1)

R FB1, min = R min – R FF

(9.2)

F

T 44

p 51

p0

T 10

43 p 0 + ------ ( p 51 – p 0 ) 50

9 p 0 + ------ ( p 51 – p 0 ) 50 Composition estimation

Composition estimation

Feedforward controller

LAG

xˆ 44

xˆ 10

KQF/RF

R FF

+ Θ3 p + Θ4 p2

Θ 1 + Θ 2 ( T + T Corr )

pˆ 44

pˆ 10

+ θ3 p

θ 1 + θ 2 ( T + T Corr )

Q FF

R min – R FF

R max – R FF

Q min – Q FF – Q FB1

Q max – Q FF – Q FB1

Q FB2, max

Q FB2

R FB1, min

Q FB2, min

PID 2

R FB1

R FB1, max

Feedback controllers with constraints

PID 1

rQ

rR

PID Q

PID R

Q

% Valve

R

% Valve

Figure 9.1: Controller implementation

+

Decoupling

Q FB1

KQ/R

+

208 9 Practical Experiences

9.3 Composition estimators

209

Q FB2, max = Q max – Q FF – Q FB1

(9.3)

Q FB2, min = Q min – Q FF – Q FB1

(9.4)

These individual constraints make unnecessary the configuration of variable structure control in the DCS. However, the maximum constraint of the reflux may lead to a top product quality significantly below the product specification, which is much more undesirable than a deterioration of the bottom product quality. Fortunately, simulation as well as practical experiences have shown that the reboiler heat duty exceeds its maximum limit Q max first. In this case, the behavior of the control scheme is identical to that of a single composition control scheme with reboiler heat duty set at maximum Q max and top composition controlled by reflux flow rate. If the reflux as well as the reboiler heat duty reach their minimum constraints, both products become purer than desired.

9.3 Composition estimators While the implementation of the controllers did not cause any particular problems, the correct parametrization of the composition estimators was very troublesome. In a first step the parameters of the estimators were calculated by regression of {Tpx} data (see Chapter 2). However, the correlation of the estimated compositions on tray 10 and 44 with the product compositions analyzed once a day proved to be unsatisfactory. Hence operating data were recorded for two weeks. Since the feed composition was almost constant, it was possible to compare these measurement data with tray compositions calculated by steady-state simulations. Minimizing the errors between the estimated and the calculated tray compositions, a correction of the estimated tray pressures by 20% was necessary to correct the estimates. Since pressure sensors on tray 10 and 44 are not installed, the pressures on these trays are calculated by a linear interpolation between top pressure and bottom pressure (see Fig. 9.1). The error in the pressure compensation

210

9 Practical Experiences

might have been caused by this interpolation. Other error sources could have been incorrect {Tpx} data or pressure measurements. Once the parameters of the estimators had been adjusted, these simple estimators worked fairly satisfactorily. Nevertheless, the compensation of the pressure variations’ influence on the tray temperatures is the limiting factor for the overall performance of the control scheme. This will be shown in more detail in the following section. Of course, the effort for the parametrization of the estimators is fairly high and the performance of the control scheme is limited by them. In view of these two points, the installation of on-line gas chromatographs could be preferable. However, in our case the light component polymerizes at temperatures exceeding a certain level. Since polymerization plugs a gas chromatograph in a short time, the use of pressure compensated temperatures or estimated tray compositions as controller inputs is indispensable.

9.4 Controller performance The controller performance observed matched the simulation results quite well. Figures 9.2 and 9.3 depict the recorded deviations of the estimated tray compositions from their setpoints in the presence of several feed flow disturbances and at two different feed flow rates. The large measurement noise of the estimated tray compositions is caused by the noisy pressure readings in the column bottom. Using a first-order low-pass element in series with the bottom pressure measurement, the noise could be significantly reduced. Unfortunately, at the time of the installation of the control scheme, the capacity of the DCS was exhausted. Even for the configuration of this simple element, there was no space left. As soon as new capacity is available, the bottom pressure measurement will be filtered. During the recording of these operating data, the setpoints were kept constant. In Figure 9.2 the feed flow rate was increased by 40 l/h in four steps. The feed flow rate at t=0 h was 260 l/h (49 mol/min), while the feed composition was approximately 0.85 mol/mol. Although the feed flow

9.4 Controller performance

211

Feed

Deviation of flow rate (l/h)

20 10 0

-10 0

10

20

30

40

50

60

70

Deviation of estimated composition (mol-%)

Tray 10

5

-0.5

0

-5 0

0.5 10

20

30

40

50

60

70

Deviation of compensated temperature (°C)

Time (h)

Deviation of estimated composition (mol-%)

Tray 44

6 -0.5 0 0.5 -6 0

10

20

30

40

50

60

70

Deviation of compensated temperature (°C)

Time (h)

Time (h) Figure 9.3: Recorded operating data with installed PID control scheme including one-way decoupling and feedforward control. Top: Deviation of feed flow rate from 170 l/h (32 mol/min) Middle: Deviations of estimated tray composition and of pressure compensated temperature from setpoint on tray 10 Bottom: Deviations of estimated tray composition and of pressure compensated temperature on tray 44

212

9 Practical Experiences

Deviation of flow rate (l/h)

Feed

40 20 0 0

20

40

60

80

100

5

-0.5

0

-5 0

20

40

Time (h)

60

80

0.5 100

Tray 44

6

-0.5 0 0.5 -6 0

20

40

60

80

100

Deviation of compensated temperature (°C)

Deviation of estimated composition (mol-%)

Deviation of estimated composition (mol-%)

Tray 10

Deviation of compensated temperature (°C)

Time (h)

Time (h) Figure 9.2: Recorded operating data with installed PID control scheme including one-way decoupling and feedforward control. Top: Deviation of feed flow rate from 260 l/h (49 mol/min) Middle: Deviations of estimated tray composition and of pressure compensated temperature from setpoint on tray 10 Bottom: Deviations of estimated tray composition and of pressure compensated temperature on tray 44

9.4 Controller performance

213

rate was out of the design range, the reflux and boilup remained within the range covered by the controller design. In Figure 9.3, the feed flow rate at t=0 h was 170 l/h (32 mol/min) and it was increased only once by 10 l/h. The control errors in presence of these feed flow disturbances remain extraordinary small. In fact, it is almost impossible to separate the control error from the measurement noise and the effect of all other unknown disturbances. This proves the high performance of this simple advanced PID control scheme. The advantages of the controller implementation are demonstrated best by a comparison of the product compositions analyzed once a day before and after the installation. At the beginning of this project, top and bottom composition were controlled manually. The results are shown on the left-hand sides of Figure 9.4 and 9.5. Obviously, the average product compositions are found far from their setpoints, and the variations of the product compositions are very large. The right-hand sides of Figure 9.4 and 9.5 show the analysis results beginning after the adjustment of the composition estimators. Clearly, the variations of the product compositions are much smaller and the average product compositions are close to the desired results. However, despite the high performance of the control scheme as illustrated by Figure 9.2 and 9.3, significant variations of the product compositions can still be detected. Please remember that pressure measurements of tray 10 and 44 are lacking. Therefore the influence of the large pressure variations (bottom pressure: 120-190 mbar) to the tray temperatures cannot be perfectly compensated and an adjustment of the controller setpoints depending on the feed flow rate is necessary. Since the results presented are achieved with almost constant setpoints, the results will improve even further as the operators gain more extensive experience with the setpoints.

214

9 Practical Experiences

Manual operation

Controlled

0.3 Top composition (mol/mol)

Top composition (mol/mol)

0.3 0.25 0.2 0.15 0.1 0.05 0 0

35 Days

0.25 0.2 0.15 0.1 0.05 0 0

70

Manual operation

0.25 Bottom composition (mol/mol)

Bottom composition (mol/mol)

44

Controlled

0.25 0.2 0.15 0.1 0.05 0 0

22 Days

35 Days

70

0.2 0.15 0.1 0.05 0 0

Figure 9.4: Analysis data of top and bottom product Top: Top composition 1-xD Bottom: Bottom composition xB Dashed line: Average composition

22 Days

44

9.4 Controller performance

215

Controlled

0.6

0.5

0.5 Relative frequency

Relative frequency

Manual operation

0.6

0.4 0.3 0.2 0.1 0

0.4 0.3 0.2 0.1 0

0 0.1 0.2 0.3 Top composition (mol/mol)

0 0.1 0.2 0.3 Top composition (mol/mol)

Controlled

0.6

0.5

0.5

0.4 0.3 0.2

Relative frequency

Relative frequency

Manual operation

0.6

0.4 0.3 0.2

0.1

0.1

0

0

0 0.1 0.2 0.3 Bottom composition (mol/mol)

0 0.1 0.2 0.3 Bottom composition (mol/mol)

Figure 9.5: Histograms of analyzed product qualities Top: Top composition 1-xD Bottom: Bottom composition xB

216

9 Practical Experiences

9.5 Economic aspects The management decision for or against the installation of a control system depends primarily on the economic feasibility and to some degree on ecological improvements. In our case the installation of the control scheme yields the following most important improvements: • More uniform product qualities ⇒ Less overpurification necessary ⇒ Energy savings (which is an ecological advantage, too) • Reduced mean of light component in bottom ⇒ More top product with a market value of > 250000 $/a • Increased maximum column load ⇒ The installation of an additional column can be avoided These pay-offs are complemented by side effects, for example a deeper understanding of column dynamics by the operating staff, which as a consequence achieved a better operation of other columns in the same plant. These financial benefits must be weighed against the investment costs. Hardware and software expenses exclusively for this project total approximately 50000 $. It is not unreasonable to estimate the necessary engineering effort for a similar project to be less than half a man year. Therefore the economic benefits are on a very positive side.

9.6 Summary The results of the implementation of the PID control structure with oneway decoupling and feedforward control on the real plant confirm the high performance of this simple control scheme indicated by simulations. The main problem of the implementation was, except for overcoming high psychological resistances, the correct parametrization of the composition estimators. A solution of this problem never would have been possible without an extensive comparison of simulation and operating data. Nevertheless, the use of pressure compensated tempera-

9.6 Summary

217

tures or estimated tray compositions remains the limiting factor of the overall performance of the implemented control scheme. The economic advantages achieved by this simple control scheme exceed the financial effort by far.

218

9 Practical Experiences

10.1 Introduction

219

Chapter 10 Conclusions and Recommendations

10.1 Introduction This thesis treats all the necessary steps for a composition control design for an industrial binary distillation column. Each of these steps produced new insights into various aspects of the control design. Since a chronological discussion of these steps would lead to a thematic confusion, they are summarized in the four sections • Controller synthesis • State-space or PID control? • How many temperature measurements? • Column models This thesis does not presume to present a final solution to all distillation control problems. The ideas presented come up against many gaps in research, limits of distributed control systems, and problems of cooperation between industry and university. In the last section the most important aspects of these topics are discussed.

220

10 Conclusions and Recommendations

10.2 Controller synthesis This thesis discusses the design of robust controllers for the dual composition control problem of an industrial binary distillation column. Distillation columns are usually operated over a wide range of feed compositions and feed flow rates. Consequently, a controller must guarantee stability and a high performance not only at a single operating point, but for the entire operating range of the distillation column. The common robust controller design methods are based on unstructured uncertainty models, for example a multiplicative uncertainty at plant output. An estimate of the corresponding uncertainty bounds has shown that these bounds are too large to allow any controller design. Nevertheless a solution of the design problem is possible. It is based on a structured uncertainty model which to a large extent avoids the unnecessary conservatism of an unstructured uncertainty description. This model treats the nonlinear column behavior as several simultaneous uncertainties and quite well describes the column dynamics for all operating points within the predefined operating range. Utilizing this uncertainty model, a feedback controller synthesis requires the framework of the structured singular value μ. The appropriate design methods are the μK-Iteration for the synthesis of statespace controllers and a constraint parameter optimization for the synthesis of controllers with fixed structure. These methods lead to feedback controllers which are distinguished by a high controller performance and guaranteed stability within the entire operating range, paired with a low sensitivity to errors in the manipulated variables. A drawback of this design approach is the high effort for uncertainty modelling and computation of the controllers. In principle, comparable results could be obtained and the computational effort could be significantly reduced by using design methods based on arbitrary small unstructured uncertainty bounds. However, it has been shown that these common design methods are not well suited for ill-conditioned plants such as high-purity distillation columns.

10.3 State-space or PID control?

221

The ideas of the feedback controller synthesis can be extended to the feedforward control design. A simultaneous controller design for the closed-loop models at maximum and minimum column load using H ∞ minimization (for state-space feedforward controllers) or optimization in the time-domain (for feedforward controllers with fixed structure) yield controllers, which greatly improve the compensation of feed flow disturbances. The theoretical and simulation results are confirmed by the results of the practical implementation of a simple PID control structure with oneway decoupling and a simple feedforward control scheme. The very satisfactory controller performance achieved without any expensive online composition analyzers leads to high economic and ecologic benefits which justify the effort of the control design and implementation.

10.3 State-space or PID control? A comparison of the different state-space controllers with optimally tuned advanced PID control structures has demonstrated an unexpected result: • The performance of μ-optimally tuned advanced PID control structures is only insignificantly worse than the performance of high-order state-space controllers This statement is of great significance for industrial practice. It holds for the feedback as well as for the feedforward control design. The implementation of advanced PID control structures in a distributed control system requires much less effort than that of state-space controllers and increases the acceptance of the control design by the operators. It must be emphasized, however, that the high performance of the PID control structure is achieved with unconventional controller settings. The optimal tuning of PID control structures with decoupling for this distillation column caused an additional insight. The optimal controller performance is achieved with an implicit decoupling scheme where in essence

222

10 Conclusions and Recommendations

• the bottom composition is controlled by moving the composition profile, and • the top composition is controlled by intensifying or weakening the S-form of the composition profile. Since the position and shape of the composition profile at steady-state depends essentially on the actual and unmeasured feed composition, it is difficult to make any inference from a composition or temperature measurement in the column middle to the manipulated variables. Very similar considerations hold for the relative performance of the state-space controllers. For the same reason, the estimation of the composition profile by the inherent observer has no advantages. The better performance results only from the higher degree of freedom in the controller design, which allows a higher performance in the low- and mid-frequency range without destabilizing the closed loop system in the high-frequency range.

10.4 How many temperature measurements? A comparison of a control design including a temperature measurement in the middle of the column with a design excluding this measurement leads to the following statement • Additional temperature or composition measurements in the middle of the distillation column have no significant influence on the maximum controller performance. The reason for the very limited advantage of additional temperature measurements for the control design is their unknown setpoint, which depends on the actual, unmeasured feed composition. The high performance of the control design can be achieved with just two pressurecompensated temperatures or two estimated tray compositions. Dispensing with additional measurements reduces the installation costs of the control system and increases its economic viability. However, if regression models are used to estimate the product compositions

10.5 Column models

223

based on temperature and flow measurements, additional temperature measurements are of great advantage.

10.5 Column models All results of this thesis are based directly or indirectly on models of the distillation column. Especially the model-based adjustment of the composition estimators clearly proved that such process models are absolutely necessary. However, the control design may be based on linear models that include or exclude flow dynamics. Within the structured uncertainty model, a multiplicative uncertainty is included for each measured tray temperature, whose uncertainty bounds exceed 100% for frequencies above 1/16 rad/min. Since the flow dynamics affect the high-frequency range, the following statement is justified: • Including or excluding flow dynamics in the linear models is insignificant for the controller design. This has an impact on the design effort. If a controller design can be based on an analytical linearization of a simple model for the composition dynamics at particular steady states, a rigorous dynamic model is not absolutely necessary. The steady states of a column may be calculated with common flowsheeting programs such as ASPEN PLUSTM or PROCESSTM and the controllers designed can be tested using a simplified nonlinear model without flow dynamics.

10.6 Recommendations 10.6.1 Academic research Multicomponent distillation: The results of this thesis are based on the example of a single binary distillation column. While the adaptation of these results to other binary columns is expected to be straightforward,

224

10 Conclusions and Recommendations

the uncertainty modelling of multicomponent distillation columns needs additional research. μ-synthesis: The robustness analysis of controllers using the structured singular value μ has shown to be a reliable and outstanding tool. However, the convergence properties of the corresponding algorithms for μ-synthesis (DK-iteration, μK-iteration) are insufficient. More robust algorithms are absolutely necessary. Decentralized control: Generally, the design of robust controllers with simple structures is at an early stage of development. In the case of this distillation column, it was relatively easy to propose potential control structures and to solve the design objective with a constrained parameter optimization (μ-optimization). However, dealing with many more control loops simultaneously, the problem of the loop pairing is still not solved. For example, the high performance of the controllers in this thesis has been obtained using the LV control configuration. Since common methods for control structure selection (single loop pairing of controlled variables and manipulated inputs) try to minimize the interactions between the individual control loops, certainly these methods favor other control configurations. Therefore methods for the selection of control structures are necessary, which include simple multivariable control schemes. Similar arguments hold for the controller tuning. The current methods for the tuning of multiloop SISO control schemes are known to be either very conservative or else to lack robustness. Better methods would be very desirable. 10.6.2 Decentralized control systems Today a control engineer in the research field is familiar with modern and flexible software tools such as MATLABTM or MATRIXXTM. His first contact with a decentralized control system (DCS), even with a modern one, arouses feelings of working in the analog computing era. The replacement of the old consoles with a computer seems to be the only idea for the development of the DCS. The inherent possibilities for a faster, more flexible, and simplified controller implementation are not exhausted yet.

10.6 Recommendations

225

10.6.3 Cooperation industry—university Often the industry complains of the inadequate cooperation between university and industry. Some typical problems of such a cooperation were encountered during the course of this project. The main problem is the divergence between the interest of the partners in the project. University researchers are interested in deeper insights into basic problems and their solution, while the process industry wants a rapid solution of the actual problem. Additionally the contact persons in industry are chronically overworked with everyday problems, thus unable to spend enough time to concern themselves with such a project. This leads to an insufficient flow of communication. Consequently, both partners speak different languages: the university researcher does not understand the industrial needs, while the industrial counterpart does not understand the mathematical methods. Therefore it is of high importance that • the aims and responsibilities of both partners in the project are spelled as clearly as possible • at least one control engineer of the industrial partner actively follows the progress of the project If these two points could be kept in mind, many problems between industry and university could be avoided.

226

10 Conclusions and Recommendations

227

Curriculum vitae

Name

Hans-Eugen Musch

Date of birth

June 19, 1965

Place of birth

Freiburg im Breisgau, Germany

Nationality

German

1971-1975

Primary school

1975-1984

Humanistic gymnasium Kolleg St. Sebastian at Stegen near Freiburg

1984

Abitur

1984-1985

Military service

1985-1989

Chemical engineering studies at the ETH Zurich

1989

Masters degree in Chemical Engineering (“Diplom”)

Since 1990

Research assistant at the Measurement and Control Laboratory, ETH Zurich

Related Documents


More Documents from "BHARGAV RAMI REDDY"

June 2020 15
May 2020 1
May 2020 0
May 2020 0