LICENTIATE T H E S I S
Department of Computer Science and Electrical Engineering
Robust Methods for Control Structure
Selection in Paper Making Processes
Miguel Castaño Arranz
ISSN: 1402-1757 ISBN 978-91-7439-188-6
Luleå University of Technology 2010
Miguel
Castaño
Ar
ranz
Rob
ust
Methods
for
Contr
ol
Str
uctur
e
Selection
in
Paper
Making
Pr
ocesses
structure selection in paper making
processes
Miguel Casta˜
no Arranz
Dept. of Computer Science and Electrical Engineering
Lule˚
a University of Technology
Lule˚
a, Sweden
Supervisor:
Printed by Universitetstryckeriet, Luleå 2010 ISSN: 1402-1757
ISBN 978-91-7439-188-6 Luleå 2010
Process industries have to operate in a very competitive and globalized environment, requiring efficient and sustainable production processes. As a result, production targets need to be translated into control objectives which are usually formulated as perfor-mance specifications of the process, i.e. tracking of references or rejection of process disturbances. This is often a hard and difficult task which involves assumptions and sim-plifications because of the process complexity. Complexity arises often due to the large scale character of a process, i.e. a pulp and paper can host thousands of control loops. A critical step in the design of these loops is the choice of the structure of the control, which means that controllers need to be placed between sensors and actuators.
Current methods for control structure selection include the Interaction Measures (IMs). The IMs help the designer to select a subset of the most significant input-output channels, which will form a reduced model on which the control design will be based. The IMs are traditionally evaluated using a nominal model of the process. However, all process models are affected by uncertainties as simplifications and approximations are unavoidable during modeling. Thus, the validity of the control structure suggested by the IMs cannot be assessed by only analyzing the nominal model. The first part of this thesis focuses in analyzing the sensitivity of the IMs to model uncertainties in order to determine a robust control structure which is feasible for all the uncertainty set.
It also becomes clear that, control structure selection requires extensive knowledge about how the multiple process variables are interconnected. The second part of this thesis focuses on creating IMs which can help the control designers to understand the propagation of effects in the process, and express this propagation in directed graphs for an intuitive understanding of the process which will help to design a feasible con-trol structure. These methods have been inspired by coherence analysis used in brain connectivity.
Neurons and neural populations interact with each other in different brain processes related to events as perception, or cognition. Electroencephalography (EEG) is a measure of electrical activity in the brain which is acquired from sensors positioned on the surface of the head, each of the electrodes collects the aggregated voltage of a neuron population. Analyzing the flow of information between populations of neurons allows to understand the communication between different parts of the brain in different brain processes. In a very similar way, analyzing the flow of information between variables in an industrial process will provide designers with the required information to understand the behavior of the plant.
Introduction
1
1 – The complexity of structures 3
2 – Control Structure Design 7
2.1 Procedure for control design of large scale systems . . . 7
2.2 Loop interaction . . . 8
2.3 Control structure selection using Interaction Measures . . . 10
3 – Contribution 11 4 – Thesis Structure 13 References 15
Part I. Robust Analysis of Process Interactions
17
I.1 – Introduction 19 I.2 – Preliminaries 21 2.1 Representing linear systems . . . 212.2 System gramians . . . 21
2.3 Norms on systems . . . 22
2.4 Index Arrays . . . 24
2.5 Representing model uncertainty . . . 25
I.3 – Bounds on Interaction Measures from parametric models with uncertainty description 27 3.1 Model uncertainty in the Nyquist diagram and connection to the bounds of tr(PjQi) . . . 27
3.2 Model uncertainty in the Bode diagram and connection to the bounds of the H2norm . . . 31
I.4 – Confidence bounds on Interaction Measures from estimated nonparametric models in the frequency domain 35 4.1 Estimation of the FRF of the system . . . 36
4.2 Robust estimation of Interaction Measures . . . 41 vii
I.5 – Paper A. Confidence bounds on tr(PQ) from estimated
nonpara-metric models in the time domain 47
5.1 Introduction . . . 49
5.2 Preliminaries . . . 50
5.3 Biased estimation of tr(PjQi) . . . 52
5.4 Unbiased Estimation of tr(PjQi) . . . 53
5.5 Confidence intervals on the estimation of tr(PjQi) . . . 55
5.6 Case study: Bark boiler . . . 56
5.7 Conclusions . . . 60
5.A Implementing the algorithm for the calculation of the confidence bounds 60 I.6 – Creating bounds for Index Arrays 65 I.7 – Conclusions and future work 69 7.1 Conclusions . . . 69
7.2 Future work . . . 70
References 71
Part II. New methods for interaction analysis and
visualiza-tion of complex processes
75
II.1 – Introduction 77 II.2 – Paper B. New methods for interaction analysis of complex processes using weighted graphs 79 2.1 Introduction . . . 812.2 Signal Flow Graphs Representing complex processes . . . 83
2.3 Scaling Independet Representation . . . 84
2.4 The Quaruple Tank Process . . . 85
2.5 Norms and normalizations for quantifying the significance of the process interconnections . . . 89
2.6 Methods for structural and functional analysis formulated from a visual-ization representation. . . 91
2.7 Interpretation . . . 99
2.8 Scaling Sensitivity . . . 102
2.9 Closed Loop Analysis . . . 102
2.10 Conclusions . . . 107
2.A Sensitivity of the Hankel norm to time delays . . . 108
II.3 – Case study: Stock Preparation Plant 111 3.1 Implementation of the stock preparation plant in ProMoVis . . . 113
3.2 Analysis of the stock preparation plant with ProMoVis . . . 113
II.4 – Conclusions and future work 115
After finishing my period as Erasmus student in LTU, I went back to my home country Spain to finalize my studies on industrial engineering. A that moment my perspectives of how to continue with my professional life were quite confusing and uncertain. Dr Wolfgang Birk encouraged me to continue my formation with Ph.D. studies and offered me a position in LTU. My first deep and sincere acknowledgments are then to Wolfgang Birk for leading me into my real vocation as a researcher, and for for his advice, support and understanding. Thank you.
Many other are those who directly or indirectly contributed to the results of this thesis. My first acknowledgements regarding contribution on results are to my collaborator Dr.
Bj¨orn Halvarsson for the fruitful discussions and the results derived.
I would like to thank to those that with Wolfgang Birk and me, created a major impact or contribution on the current prototype of the software tool ProMoVis, with special mention to Dr. Johan Karlsson, Tilda Nordin, and Dr. Andreas Johansson.
Thanks also to Pablo Fern´andez de Dios for the first implementation of the interface
using Simulink.
The collaboration of the industrial partners of the SCOPE consortium is gratefully ac-knowledged, with special mention of the participants in MeSTA. Thanks to SCA Obbola AB for the experiments performed on the stock preparation plant and the bark boiler, which were of great importance for the illustration of some of the methods included in this thesis.
I would like to express my gratitude to Professor Johan Schoukens and the depart-ment ELEC in Vrije Universiteit Brussel (VUB) for the school Measuring, Modelling and Simulation of (Non)linear Dynamic Systems.
Regarding funding, the contributions provided by Swedish Governmental Agency for Innovation Systems (VINNOVA) and the Hjalmar Lundbohm Research Center (HLRC) are gratefully acknowledged.
In my work environment, I would like to thank to the numerous colleges and friends from the Department of Computer Science and Electrical Engineering at LTU for this friendly environment and full of stimulating discussions. Special thanks to my office mate and friend Soheil Salehpour for bearing with my daily eccentricities.
From the personal perspective, I want to say thank you to my parents and to Maite
Nicol´as for their faithful support without which I could have not done this work.
Lule˚a, November 2010
Miguel Casta˜no Arranz
The work conducted in this thesis has been done in the MeSTA project (Methods for Structural Analysis in Pulping).
The project MeSTA is run by ProcessIT Innovations under the SCOPE industry programme for the pulp and paper industry.
As a result of the project, a prototype of a software tool for Process Modeling and Visualization (ProMoVis) has been created.
ProMoVis is being used in the paper mills of SCA Obbola AB and Billerud Karlsborg AB since 2009. ProMoVis can visualize the mill processes and the control system, and it has been proved to be useful in obtaining and transferring process knowledge, and in helping the control engineers to generate feasible control structures.
ProMoVis is still to be further developed and is meant to include the methods for robust control structure selection and for analysis and visualization of complex processes created and developed in MeSTA until the end of the project in December 2011.
The complexity of structures
The control engineer of a complex industry process such as a pulp and paper mill faces similar challenges to those of a builder of a Gothic cathedral in medieval times.
A crucial competence of a master of works in medieval times was mastering the art of controlling the balance of forces in a delicate and complex structure (see Fig. 1.1). For this purpose, the architect had at his disposal several elements and strategies such as arches, vaults, timber frame roofs, columns and counterforts [1]. A key step on the design was to make structural decisions on where to place this elements for an optimal distribution of forces. In the same way, a control engineer has at his disposal different control elements differentiated as control strategies such as PID, MPC, H∞, Feed-forward actions, LQG, . . . , and this elements have to be structurally placed for obtaining an optimal distribution of control actions in a largely complex system which is the pulp and paper mill (see Fig. 1.2).
Figure 1.1: Structure of the Gothic cathedral Notre Dame of Rhemis (source [1]).
4 The complexity of structures
Figure 1.2: The Kraft process (source Metso). Simplification of a possible structure of a pulp and paper mill.
By structurally placing and designing these elements, the builder aimed for certain objectives like maximizing the open areas at the facade for placing windows and doors, minimizing the number and with of columns in the naves for increasing the number of followers which can fit in and the visibility they experience, maximizing the height and slenderness of the building, whilst minimizing the cost of the resources. In addition, the construction needs to be stable and robust against disturbances such as wind (continuous disturbance) and possible attacks with siege weapons (instantaneous disturbance). The control engineer of a modern pulp and paper mill also need to choose where to structurally place the controllers in a complex structure, and design them to fulfill certain control objectives, like maximizing production and certain quality parameters whilst minimizing the consumption. This has to be done while keeping a stable process which is also robust to certain continuous or instantaneous process disturbances, and even failure situations. In both examples, the structural decisions have to be taken in a largely complex system in which local actions or decisions can derive in unexpected consequences and failures at other parts of this interconnected structure. Therefore, even if it is possible to introduce hierarchies for decomposing, analyzing and structurally designing such a complex system, a holistic perspective has to be always considered.
The knowledge of a master of works was much more a know-how than a theoretical knowledge. In order to create a stable and robust structure, the master had to take risks and make experiments, always bearing in mind former mistakes and accidents, like collapsing of vaults [1]. The advances in elemental theory research on physics and the wide range of computer simulating tools allows the modern architects to take a much more theoretical approach, in contrast with the profoundly empirical methods used by the medieval builder. For the actual framework of process control, and even if there exists
a wide selection of methods for control structure selection, being the first one dating form 1966 ([2]), these methods are usually published in research papers and rarely form part of educational books or courses. However, there has been a large research contribution in the field of control structure selection during the last decade, and this topic starts to be addressed as sections in educational books about multivariable control ([3]), or even having full dedicated books ([4]). This large delay between research and education, and finally industry application ([5]) implies that, control engineers in current paper mills still address a strongly empirical approach to control structure design, basing their decisions in know-how or common sense principles and experience.
This is the motivation for one of the objectives of the MeSTA project, which is to bring these methods for control structure selection into a software tool which will help the control engineers in paper mills and other complex industrial plants to take robust decisions in the control structures to be used. A prototype of this software tool has already been created and named ProMoVis ([6]) and is still evolving with late results on control structure selection.
Far from being only a software development project, MeSTA is also involved in re-searching methods for robust control structure design. Robustness is needed to deal with the uncertainty and inaccuracy which are inherent properties of any model. Current methods for control structure design are based only on a nominal model of the analyzed system, but there is clearly a need of creating robust methods which would validate a decision considering the existing model uncertainty. The robustness of the methods for control structure selection is an issue which has recently received increasing attentions, like the computation of the RGA ([7]), the HIIA ([8]) or the Gershgorin bands ([9]) for uncertain process models. However, a full consolidated method for selecting a robust control structure is still sought.
Continuing with the used analogy, it was soon natural in this world of medieval complexity, that it was necessary to have good tools for communication between the patron, the master of works, and other involved parts. To understand each other, plans, diagrams and working models were used on the building site like the one in Fig. 1.3. Those visualization techniques evolved and became more sophisticated, supporting a good understanding and communication and deriving in a more efficient work.
Figure 1.3: Distribution of forces in a vault (source [10]).
Figure 1.4: Importance of the actuators for each measure-ment in the refining section of a stock preparation plant.
6 The complexity of structures
Visualization and communication tools and techniques are also needed in a modern pulp and paper mill. Diagrams and flow sheets are used by control and process engineers for communicating process knowledge to other workers such as operators or to decision
boards. It is an objective of MeSTA to explore and evolve visualization techniques
and combine them with mathematical tools in order to create diagrams which allow to understand and communicate the process behavior (see Fig. 1.4).
An important added value of such a tool is the ability to maintain and transfer process knowledge. In both cases of this analogy, if the responsible of evolving and maintaining the complex system is substituted by a new employee, the new incorporate is likely to be lost in this world of complexity which includes solutions and structures which have rarely been documented, and whose sense was in the mind of the previous designer. ProMoVis, within the same software tool, allows to create and edit process diagrams, to define process models and analyze them for selecting an optimal control structure, and to define implemented control solutions such as soft sensors and controllers and use this information for the analysis of interactions in the closed loop. It is clear that, having at disposal well maintained process information in ProMoVis will largely reduce the integration period of a new control engineer in a complex industry plant such as a pulp and paper mill.
Control Structure Design
Control structure design is a step in the procedure of the design of a control system in which actuators and sensors are grouped in sets of variables with strong interconnections. Each of the actuator/sensor sets will determine a subsystem, and independent controllers will be designed for each of the subsystems. A proper structure with well tuned controllers should derive in low loop interactions which would otherwise derive in a performance degradation.
2.1
Procedure for control design of large scale
sys-tems
The design of a control system for a multivariable process usually involves the following steps [11, 12]:
- Identify production targets.
- Select a subset of sensors and actuators which will be used for control. - Determine a model for the system.
- Select a control structure by selecting a subset of the most significant actuator-measurement interconnections on which control will be based.
- Translate production targets into performance specifications for the controller. - Synthesis of the controller parameters using an adequate control strategy (i.e. PID,
MPC, LQG, . . . ).
- Controller implementation.
The work in this thesis focuses in the design step of the control structure selection. The objective is to select the simplest control structure which derives in low loop inter-action.
8 Control Structure Design
2.2
Loop interaction
Interaction analysis plays an important role in large scale control systems. The design of a control structure which is not considering important process interconnections will most likely derive in a large loop interaction having a performance degradation as direct consequence.
Figure 2.1: Interacting subsystem formed by the pulp recirculation conducts in a Stock Prepa-ration Plant.
Consider the simple case in which in which the recirculated flow of pulp from two parallel refiners in a stock preparation plant of a paper mill is to be controlled (see Fig. 2.1). The recirculated flow from each of the refiners is measured by a flow indicator
(F I1 and F I2), being these flows controlled by two independent SISO controllers (C11
and C22) for tracking the references (R1and R2) which are set by and operator. Each of
the controllers decide the control action on a valve (V1 and V2 respectively) in order to
maintain the desired flow references. A linearization of the process can be described as a multivariable transfer function of the form:
F I1 F I2 = G11 G12 G21 G22 | {z } G V1 V2
where Gij is the transfer function from Vj to F Ii. The control action on the actuators
can be described as a function of the reference and the measurements as:
V1 V2 = C11 0 0 C22 | {z } C R1− F I1 R2− F I2
where the matrix C represents the controller. The structure here resulting from a di-agonal controller matrix C is known as decentralized control. In a decentralized control structure, the process actuators and controlled variables are grouped in pairs, and SISO controllers are synthesized for each of the pairs. This structure is usually the preferred one in process industry due to its simplicity of design and maintenance. However, a de-centralized control structure assumes a free of interaction process, which means that the off-diagonal terms of the transfer function matrix G are 0 (see Fig. 2.2). In our example this is clearly a simplification of the real case, and the existence of off-diagonal elements in G will involve a performance degradation due to loop interaction which will depend on their significance. If the off-diagonal elements of G are very small in comparison with the diagonal ones, and for a well tuned decentralized controller, then the process behaves almost like SISO independent loops ([13]).
For a better understanding of the nature and importance of loop interaction, two cases will be differentiated ([14]) using the example of the flow recirculation.
- One way interaction. This is the case when, for example, G21 6= 0 and G12 = 0 (see
Fig. 2.3). It is clear that a reference change in R2 will not perturb the other loop,
but a reference change in R1will require a control action in V1 which will perturb
F I2 and make the controller C22 react and actuate con V2 in order to reject the
interaction disturbance. A possible solution for this would be to add a feed forward action in the second loop which will compensate for the interaction disturbance.
This will be reflected in a nonzero element C21in the controller matrix C.
- Two ways interaction. Both off-diagonal terms are significant (see Fig. 2.1).
Consid-ering a change of the reference R1, the existing path G21 will, as in the previous
example, create a perturbation in the other control loop which will require an actu-ation on V2, but in addition, this will involve another perturbactu-ation on the original
loop through the existing path in G12. This loop interaction is likely to create large
oscillations in the control loops, and can even affect the closed loop stability. In this example a full MIMO controller would be usually required, which is reflected in a controller matrix C with no zero elements.
Figure 2.2: Example of a free of inter-action subsystem.
Figure 2.3: Example of a subsystem presenting one way interaction.
10 Control Structure Design
The complexity of the interaction effects grows with the number of process variables. In this example only the recirculation system of a stock preparation plant is considered but this is only a subsystem which is interacting with the rest of the stock preparation plant, which is in turn another subsystem of superior hierarchies.
The complete plant of a large scale system is usually decomposed in reduced subsys-tems with low reciprocal interaction. The resulting structure of a plantwide controller is sought to be a block diagonal, with each of the blocks being the designed controller of a subsystem in the decomposition.
2.3
Control structure selection using Interaction
Mea-sures
Control structure selection is the act of selecting the nonzero elements in the controller matrix C. This structure is preferred to be as sparse as possible whilst presenting low loop interaction. The maximum sparsity is obtained with a decentralized controller, which is likely to derive in large loop interaction in the case of a process with large off-diagonal elements. The maximum complexity will be a full MIMO controller with
no zero elements in the controller matrix C. This structure will have the potential
of synthesizing a controller with the minimum achievable loop interaction, whilst the design, implementation and maintenance of such a controller would present a large cost and complexity.
Interaction measures (IM) help the designer to select a control structure by selecting a subset of the most important input-output interconnections which will for a reduced model on which control will be based.
A large amount of literature has been published on IMs since Bristol introduced the Relative Gain Array (RGA) in 1966 ([2]) as an indicator for choosing input-output pairings in decentralized control structures which is based on the DC gain of the process. The RGA has been largely used since then, and many authors have proposed different variants to deal with some of its limitations, like the dynamic RGA (DRGA) which is an extension of the original concept of the RGA in the frequency domain and is able to consider process dynamics [15]. The RGA can provide with feasible suggestions for choosing a decentralized pairing, and in combination of the Niederlinski Index [16], can verify the stability and integrity of the chosen pairing.
However, the RGA gives no information on how to increase the controller complexity in the case of having a process where a decentralized control structure will derive in large loop interaction. The Index Arrays (IA) where created to overcome with this limitation and are able to suggest sparse control structures. IAs exploit the concept of overall used dynamics of a system. This means that the index array gives information on the amount of the dynamics described by the reduced model which is used in the controller design. Or in other words, if a large amount of the dynamics are considered during the controller design, the there is a large probability that the control structure is also viable.
Contribution
During 2008 and 2010 the author of this thesis contributed in 6 international publi-cations. All the contributions are within the field of control structure selection, but two areas can be differentiated.
Robust analysis of process interactions
Paper 1. M. Casta˜no and W. Birk, “A new approach to the dynamic RGA analysis
of uncertain systems,” IEEE International Conference on Computer-Aided Control Systems (CACSD), September 2008.
Paper 2. B. Halvarsson, M. Casta˜no, and W. Birk, “Uncertainty bounds for
gramian-based interaction measures,” The 14th WSEAS International Conference on Sys-tems, July 2010.
Paper 3. M. Casta˜no, W. Birk and B. Halvarsson, “Unbiased estimation of
gramian-based process interactions with confidence bounds,” Submitted to IFAC World Congress 2011.
Paper 1 describes a sensitivity analysis of the RGA to model uncertainty. Its contents are not included in this thesis, since it is a formalization of the results produced by the author in his Master Thesis [17], which can be consulted for a more detailed description and variety of examples.
Paper 2 includes a pre-study in which a previously existing inequality for creating bounds on the squared Hilbert-Schmidt norm of a SISO subsystem [8] is tightened and compared with a new proposed method. The main contribution of the author was to
create and develop this new method in collaboration with Bj¨orn Halvarsson and Wolfgang
Birk. The considered examples in the paper were favorable to the new method, which has been evolved in this thesis and included in chapters I.3 and II.1.
Paper 3 has been included in this thesis as Paper A (Chapter I.5). The contribution of Wolfgang Birk with the models from the Bark Boiler, supervision, advice, reviewing
and writing efforts is gratefully acknowledged. The contribution of Bj¨orn Halvarsson in
12 Contribution
the creation of the method described and in the reviewing of the manuscript is gratefully acknowledged. Special thanks to SCA Obbola AB for the experiments performed on the bark boiler, and from which the models used in the paper were derived.
New method for structural and functional analysis of complex processes
Paper 4. M. Casta˜no and W. Birk, “New methods for structural and functional
anal-ysis of complex processes,” International Conference on Computer-Aided Control Systems (CACSD), July 2009.
Paper 5. W. Birk, A. Johansson, M. Casta˜no, S. R¨onnb¨ack, T. Nordin, and N.-O.
Ekholm, “Interactive modeling and visualization of complex processes in pulp and paper making,” Control Systems 2010, September 2010.
Paper 6. M. Casta˜no and W. Birk, “New methods for interaction analysis of complex
processes using weighted graphs,” Submitted to Journal of Process Control. Paper 6 is a journal extension of Paper 4. Paper 6 has been included in this thesis as Paper B (Chapter II.2). The contribution of Wolfgang Birk with supervision, advice, reviewing and writing efforts is gratefully acknowledged.
In Paper 5, the contribution of the author was creating a case study in which the stock preparation plant in SCA Obbola AB is analyzed. This contribution has been extracted from the paper and included as Chapter II.3 of this thesis. The review and advice of the rest of the authors in creating this case study is fully acknowledged. Special thanks to SCA Obbola AB for the information provided on the control structure of the stock preparation plant, and for the experiments performed from which models were created and analyzed in this case study.
Thesis Structure
This thesis Introduction includes in Chapter 1 the existing background in the design of control structures in process industry, which is the motivation for the conducted re-search. A description of the problem of control structure selection is given in Chapter 2. The contribution of this thesis to the filed of control structure selection is described in Chapter 3.
The results reported in this thesis were divided in two parts with the following con-tents:
Part I. Robust Analysis of Process Interactions.
The robustness of two previously existing IMs for control structure selection is analyzed. New methods for either determining bounds on the considered IMs from uncertain models or for creating confidence estimations on the IMs are introduced. These methods allow to take robust decisions on control structure selection. Part II. New methods for interaction analysis of complex processes using
weighted graphs.
The introduced methods give visual and intuitive information on the process struc-tural and functional properties. These methods can be used for acquiring and transferring process knowledge, and have direct applications in the design of con-trol structures.
Each of the parts and the introduction has their own reference list at their end. Two papers have been inserted as chapters inside the parts, conserving their own and independent references, which are listed at the end of the paper.
The numbering of equations, figures, chapters and sections resets for each of the parts and the introduction.
Examples are numbered only within the included papers, and they maintain an inde-pendent numbering.
There is a list of acronyms at the end of the thesis.
[1] M. Henry-Claude, L. Stefanon, and Y. Zaballos, Principles and elements of medieval
church architecture in Westerdn Europe. Fragile, 2007.
[2] E. Bristol, “On a new measure of interaction for multi-variable process control,” IEEE Transactions on Automatic Control, vol. AC-11, no. 1, pp. 133 – 134, 1966.
[3] S. Skogestad and I. Postlethwaite, Multivariable Feedback Control, 2nd ed. John
Wiley and Sons, 2005.
[4] A. Khaki-Sedigh and B. Moaveni, Control Configuration Selection for Multivariable
Plants. Springer, 2009.
[5] “Future perspectives of distributed control systems software for the process indus-try,” The fourth Saudi engineering conference, vol. III, November 1995.
[6] W. Birk, A. Johansson, M. Casta˜no, S. R¨onnb¨ack, T. Nordin, and N.-O. Ekholm,
“Interactive modeling and visualization of complex processes in pulp and paper making,” Control Systems 2010, September 2010.
[7] D. E. S. Dan Chen, “Relative gain array analysis for uncertain process models,” AIChE Journal, vol. 48, no. 2, pp. 302–310, 2002.
[8] B. Moaveni and A. Khaki Sedigh, “Input-output pairing analysis for uncertain mul-tivariable processes,” Journal of Process Control, vol. 18, no. 6, pp. 527 – 532, 2008. [9] D. Chen and D. Seborg, “Robust Nyquist array analysis based on uncertainty de-scriptions from system identification,” Automatica, vol. 38, no. 3, pp. 467 – 475, 2002.
[10] K. D. Alexander, R. Mark, and J. F. Abel, “The structural behavior of medieval ribbed vaulting,” Journal of the Society of Architectural Historians, vol. 36, no. 4, pp. pp. 241–251, 1977.
[11] M. Morari, Y. Arkun, and G. Stephanopoulos, “Studies in the synthesis of control structures for chemical processes. I. Formulation of the problem. process decompo-sition and the classification of the control schemes.” AIChE J., vol. 26, no. 2, pp. 220–232, 1980.
16 References
[12] M. van de Wal and B. de Jager, “A review of methods for input/output selection,” Automatica, vol. 37, no. 4, pp. 487 – 510, 2001.
[13] M. E. Salgado and A. Conley, “Mimo interaction measure and controller structure selection,” International Journal of Control, vol. 77, no. 4, pp. 367–383, 2004. [14] W. Birk, “Towards incremental control structure optimisation for process control,”
dec. 2007, pp. 1832 –1837.
[15] M. Witcher and T. McAvoy, “Interacting control systems: Steady-state and dynamic measurement of interaction,” ISA Transactions, vol. 16, no. 3, pp. 35–41, 1977. [16] A. Niederlinski, “A heuristic approach to the design of linear multivariable
interact-ing control systems,” Automatica, vol. 7, pp. 691–701, 1971.
[17] M. Casta˜no, Sensitivity of variable pairing in multivariable processcontrol to model
Robust Analysis of Process
Interactions
Introduction
A process model with a description of the model uncertainty can be understood as a set of models which includes the nominal model. Therefore, the value of the Interaction Measures (IMs) may differ for different models in the set. Including model uncertainties in the interaction analysis would make possible to take robust decisions in the control structure selection.
Two interaction measures are considered in this thesis: the Σ2 [1] and the PM [2].
These are Index Arrays (IAs) created from the H2norm and the squared Hilbert-Schmidt
(HS) norm respectively. The objective of this part is to create methods for computing the possible bounds on these norms, which can later be converted in the bounds of their corresponding IA as described in Chapter I.6.
Chapter I.3 describes methods for calculating the bounds on the considered norms for uncertain parametric process models affected by multiplicative uncertainty.
However, creating parametric models is usually a complicated and time consuming task. This complexity increases as the number of process variables increases, and esti-mating nonparametric models is usually a more general and simpler approach which can be addressed without the need of choosing the model order and structure.
Chapter I.4 introduces a methodology for obtaining bounds on the estimation of
the H∈ and the squared HS norms of each of the input-output channels of a linear
multivariable system. The robust estimation of the considered norms is performed in the frequency domain from a nonparametric model of the Frequency Response Function (FRF).
A large advantage of this method is derived from existing modeling approaches for the identification of nonlinear systems. Traditional methods for control structure selection require linear models for their use. In the case of facing nonlinear systems, IMs are usually applied to a linearized process model around a working point, without any quantification of the degree of nonlinearity of the systems which would validate the decision of analyzing
it within the linear framework. By using periodic excitation signals with a tailored
spectra, the degree of contribution of the process nonlinearities can be quantified [3], and the best linear approximation of the system can be estimated as a nonparametric
20 Introduction
model of the FRF. It is suggested in this thesis to use such an excitation design in case of having to perform an analysis of interactions of a nonlinear process with unknown degree of nonlinearity. If, after quantifying the contribution of the process nonlinearities it is found that the system can be considered as weakly nonlinear, then the obtained
best linear approximation can be used to compute a robust estimation of the H∈and the
squared HS norms which will allow to take robust decisions in control structure selection. If the process is found to be strongly nonlinear, then any of the existing methods for control structure selection for nonlinear systems can be applied.
Possible extensions of the RGA concept to nonlinear systems have been proposed. The extension proposed in [4] assumes that the process outputs can be expressed as an invertible function of the process inputs. Another possible extension was introduced in [5], requiring a non-linear parametric model of the process, and using feedback linearization to compute the DC-Gains of the non-linear process. However, this measurements, as extensions of the RGA, are only able to suggest decentralized control structures. A large work in the interaction analysis of nonlinear systems still remains to be done.
Chapter I.5 includes an alternative to the estimation of the squared HS norm in the frequency domain introduced in Chapter I.4. The method created confidence bounds on the squared HS norm based on estimation of a nonparametric model of the system impulse response, which is obtained by a linear regression on process data in the time domain obtained with white noise input excitation. This estimation method is simpler and general.
Using any of the approaches to obtain a robust estimation of the considered norms in a multivariable process will allow to make a robust selection of a subset including the most significant input-output interconnections on which the controller will be based. In the case of requiring more sophisticated models for the controller synthesis, the modeling effort will only be placed on this set of significant interconnections.
Preliminaries
2.1
Representing linear systems
A time-invariant, linear continuous dynamic system can be represented in a state space representation of the form:
˙x = Ax + Bu (2.1a)
y = Cx + Du (2.1b)
The process inputs, outputs and internal states are collected in the vectors u ∈ Rn,
y ∈ Rm and x ∈ Rp respectively, and the system is represented by the quadruple
(A, B, C, D) with A ∈ Rp×p, B ∈ Rp×n, C ∈ Rm×p and D ∈ Rm×n.
A multivariable continuous time linear process can also be represented by its transfer function matrix G(s), which for a proper system, can be obtained from the state space representation as:
G(s) = C(sI − A)−1B + D
where s is the Laplace variable, and the size of the transfer function matrix is m × n.
Each of the elements Gij(s) of the matrix G(s) is the elementary subsystem which
con-nects the jthinput with the ithoutput. When the system is represented in the frequency
domain, the Frequency Response Function (FRF) is denoted as G(jω).
2.2
System gramians
Given a continuous time system in state space form as in Equation (2.1), the controllabil-ity (P) and observabilcontrollabil-ity (Q) gramians are obtained by solving the following continuous-time Lyapunov equations [6]:
AP + P AT+ BBT = 0
22 Preliminaries
ATQ + QA + CTC = 0
The controllability gramian quantifies how hard it is to control a system state from the system inputs, and the observability gramian quantifies how hard it is to observe a system state from the system outputs. Some properties of the gramians are [7]:
- The minimal energy required to steer the states of the system from 0 to xf is
x∗ fP
−1x f.
- The maximal energy of the output obtained by observing a system with initial state x0, is x∗0Qx0
- The states which are difficult to reach are in the span of the eigenvectors of P which correspond to small eigenvalues, and the states which are difficult to observe are in the span of the eigenvectors of Q which correspond to small eigenvalues.
Therefore, inspecting the eigenvalues of P quantifies how hard it is to control the system states from the system inputs, and inspecting the values of Q quantifies how hard it is to observe a system state from the system outputs.
Gramians can be used to assess the process interactions and select viable control structures, giving birth to several interaction measures. An early contribution in the field of gramian based IMs dates from 1996, where A. Khaki-Sedigh and A.Shahmansourian [8] presented an input-output pairing algorithm based on the cross-gramian matrix [9]. The algorithm was using an indicator which was later classified as an IM named Dynamical Input-Output Pairing Matrix (DIOPM) [10]. However, the term gramian based IMs became popular after the publication in 2000 of the work by A. Conley and M. Salgado [11], where the Participation Matrix (PM) was first introduced. Other gramian based
IMs were later introduced. In 2002 B. Wittenmark and M. Salgado introduced the
Hankel Interaction Index Array (HIIA) [12], which is the same indicator as the DIOPM, but reinvented using the largest singular value of the product P Q instead of the largest singular value of the squared cross-gramian matrix to quantify the channel significance, being both values equal. In 2003 W. Birk and A. Medvedev [1] suggested the use of other gramian based norms for the quantification of process dynamics, and introduced the IM
denoted as Σ2.
2.3
Norms on systems
In the work reported in this thesis, the following norms are used for quantifying the importance of the process dynamics.
The squared Hilbert-Schmidt norm
The squared Hilbert-Schmidt (HS) norm is equal to the sum of the squared singular values of the Hankel operator (Hankel Singular Values) [13]. The Hankel operator maps
the past system inputs into the future system outputs, and the Hankel Singular Values (HSVs) have important implications in quantifying the system dynamics [7].
The squared nonzero HSVs are equal to the eigenvalues of tr(PQ) and therefore, the squared HS norm can be computed as the trace of the product P Q:
||Gij||2
HS= tr(PjQi)
It has been mentioned that the gramians contain important information on the sys-tem dynamics. However, gramians depend on the state realization. To overcome this dependance, the product P Q can be used to quantify the system dynamics. P Q collects important information about the controllability and observability of the system [11], be-sides it is a positive-semidefinite matrix whose eigenvalues don’t depend on the state realization.
Relationship of tr(PQ) with the frequency domain. The Nyquist diagram is a
well known tool for the analysis of dynamic systems. It was shown in [13] that tr(P Q) is equal to the area enclosed by the oriented Nyquist diagram divided by π, with oriented meaning clockwise direction:
tr(PjQi) = π−1Ac(Γij(ω)), (2.3)
where Ac(.) means area enclosed by, andΓij(ω) = Gij(jω)with ω ∈ Revaluated from
−∞ to ∞ in the continuous-time case.
Relationship of tr(PQ) with the time domain. In this thesis, the following
rela-tionship will be used for sampled systems to compute the If a discrete impulse response for each of the channels of a multivariable system is know, tr(PjQi) can be computed as:
tr(PjQi) = N X
k=0
k (hij(k))2
were hij(k) are the coefficients of the impulse response, k is the order of the coefficients,
and N is the length of the impulse response.
The H2 norm
For any of the Gij(s) SISO systems, if Gij(s) is stable and strictly proper, ||Gij(s)||2can
be expressed as ||Gij(s)||2= s 1 2π Z ∞ −∞ |Gij(jω)|2dw (2.4)
The H2norm has a strong connection with the controllability and observability
grami-ans. In the case that the process is given by state space description (A, B, C, 0), the H2
norm of each of the elementary SISO subsystem (A, Bj, Ci) can be computed as
||Gij(s)||2=
q
CiPjCT
24 Preliminaries
where Pjis the controllability gramian of the SISO subsystem, and Cithe i-th row of C.
It has been shown [14] that the squared H2norm of each elemental SISO subsystem
can be interpreted as the coupling in terms of the energy transmission rate between the past inputs and the current output.
2.4
Index Arrays
In Section 2.3 the H2and the squared HS norms have been introduced. Both of them can
be used to quantify the importance of each of the input-output channels in a multivariable system. The Index Arrays were defined as an array including the norms of each of the elementary subsystems divided by the total sum of the norms of all the elementary subsystems.
Note that, due to the normalization, all the elements of given index array add up to 1. Therefore the larger elements in an index array identify the input-output channels which have a higher contribution in the process. When an IA is used for control structure design, a subset S of the input-output channels is selected as the most significant. By evaluating
the closeness ofP(i,j)∈Si=1j=1tr(PjQi) to 1, the designer can comprehend the amount of the
total process dynamics that the subsystem according to S is reflecting. The objective is to find a simplified model formed by a reduced subset of the process interconnection. This model will be used for the controller design and have to capture most of the process dynamics.
The two considered index arrays in this thesis are the PM and the Σ2. Other popular index arrays are the DIOPM [10] and the equivalent HIIA [12].
The Participation Matrix (PM) The PM was defined as [11]:
φij=
tr(PjQi) X
k,l
tr(PlQk)
The unnormed PM, will be denoted as ˜φ, and it is an array collecting the values of
tr(PjQi) for each of the elementary input-output channels.
˜
φij= tr(PjQi)
In this thesis, the names unnormed PM, squared HS norm of each SISO subsystem
and tr(PjQi) are used interchangeably. Which means that also the following notation is
interchangeably used:
˜
φij= ||Gij||2HS= tr(PjQi)
Note from Equation (2.3) that the unnormed PM is proportional to the area enclosed by the Nyquist diagram of each of the input-output channels. This shows that the PM is in fact closely related to the Direct Nyquist Array (DNA) introduced by Rosenbrock in the early 1970:s (see for instance [15, 16] for an introduction to DNA analysis and
[17] for a discussion of DNA in connection with uncertain systems). In the basic DNA approach, Nyquist curves are plotted for each subsystem and the decentralized pairings corresponding to the largest Nyquist curves are selected. Obviously, the PM is a quan-titative version of this concept, however, initially derived to quantify controllability and observability of the state space.
The Σ2 Interaction Measure
The Σ2 is defined as [1]: [Σ2]ij= ||Gij(s)||2 X k,l ||Gkl(s)||2
The unnormed Σ2, will be denoted as ˜Σ2, and it is an array collecting the values of
||Gij(s)||2for each of the elementary input-output channels.
[ ˜Σ2]ij= ||Gij(s)||2
2.5
Representing model uncertainty
To represent model uncertainty in a SISO system, the following notation will be used: Π : uncertainty set. Includes all the possible plants due to uncertainty
G(s) ∈ Π : nominal plant.
Gp(s) ∈ Π : particular perturbed plant.
The nominal plant G(s) is the actual model of the plant, that is, what we assume that the plant is; Gp(s) represents a particular possible plant due to uncertainty, and Π repre-sents the set of all possible plants. The uncertainty set is described using multiplicative uncertainty as:
Π : Gp(s) = G(s)(1 + W (s)∆W(s)); |∆W(jω)| ≤ 1, ∀ω (2.5)
W (s) is the scaling factor, and it is a stable transfer function selected to represent the
uncertainty, ∆W(s) represents any stable transfer function with magnitude less or equal
than one at each frequency. G(s) 1 + W (s)∆W(s) describes at each frequency ωk a
circular region of uncertainty centered in G(jωk) with radius equal to |G(jωk) · W (jωk)|.
As illustrated in Fig. 2.1, W (s) has to be selected so that the uncertainty set Π includes
the possible uncertainty at each frequency ωk.
Another possibility is to describe the uncertainty set as additive uncertainty of the form:
Π : Gp(s) = G(s) + Wa(s)∆Wa(s); |∆Wa(jω)| ≤ 1, ∀ω (2.6)
In this case, the circular regions describing the uncertainty have as radius |Wa(jωk)| at
26 Preliminaries
Figure 2.1: Circular uncertainty regions generated by multiplicative uncertainty. W (s) has to be selected to include at each frequency all the possible perturbed plants due to uncertainty.
uncertainty as relative to the nominal case. However, both notations are equivalent if
|W (jω)| = |Wa(jω)|/|G(jω)|.
The uncertainty weights (W (s) or Wa(s)) are usually designed as rational transfer
functions. However, other options are possible, like defining the weight as a real valued
function of ω: W (ω) or Wa(ω) [18]. For the results presented in Chapter I.3 the latter is
preferred, since the results will be based on integrating expressions which are a function of the uncertainty weights, and when using computer toolboxes based on symbolic math operations it will in most cases involve less computational complexity.
For the results presented in this thesis the following conditions will be required for the weights of the used multiplicative or additive uncertainty:
lim
ω→∞W 6= ±∞ ; ω→∞lim Wa= 0
Bounds on Interaction Measures
from parametric models with
uncertainty description
In this chapter bounds for the H2and the squared HS norms are created for parametric
process models with uncertainty description. MIMO models will be considered, in which each of the SISO subsystems is independently perturbed by multiplicative uncertainty as in Equation (2.5). The effect of uncertainty in the considered norms will be determined independently for each of the SISO subsystems. Therefore, all the formulas in this chapter are applied independently to each of the SISO subsystems, but the term ij will be dropped for the sake of simplicity.
3.1
Model uncertainty in the Nyquist diagram and
connection to the bounds of tr(P
jQ
i)
It has been mentioned in Section 2.5 how the multiplicative and additive descriptions of model uncertainty affects the Nyquist diagram. It was also mentioned in Section 2.3,
that the value tr(PjQi) for each of the SISO subsystems of a MIMO system can be
computed as the area enclosed by the oriented Nyquist diagram divided by π. Therefore, by obtaining the minimum and maximum possible areas enclosed by the Nyquist diagram, the possible variation of tr(PjQi) in the uncertainty set can be obtained.
Considering a MIMO system in which each of the SISO subsystems is affected by independent multiplicative uncertainty, the uncertainty region in the Nyquist diagram can be determined for each of the SISO subsystems using the Frenet-Serret frame, which is a well know concept in differential geometry.
For a SISO model affected by multiplicative uncertainty as in Equation (2.5),
gen-erally, adding the value of|G(jω)W (jω)| in the direction and sense of the vector normal
to the nominal Nyquist curve (see Fig. 3.1), the curve enclosing the smaller area in the 27
28
Bounds on Interaction Measures from parametric models with uncertainty description
Figure 3.1: Nyquist diagram perturbed by uncertainty. The nominal Nyquist curve, and the curves enclosing the minimum and maximum areas in the uncertainty set are depicted.
uncertainty set is obtained. Similarly, by adding|G(jω)W (jω)|in the opposite sense, the
curve enclosing the larger area in the uncertainty set is obtained. Denoting G0(jω) = ∂G(jω) ∂ω , G 00(jω) = ∂2G(jω) ∂ω2 ~ GN(jω) = G 0(jω) × (G00(jω) × G0(jω)) |G0(jω)| · |G0(jω) × G00(jω)| (3.1)
where × denotes the vector product, then the the curves Glow and Ghigh enclosing the
minimum and maximum possible areas can be described as:
Glow(jω) = G(jω) + ~GN(jω) · |G(jω)| · |W (jω)| (3.2a)
Ghigh(jω) = G(jω) − ~GN(jω) · |G(jω)| · |W (jω)| (3.2b)
However, it has been observed that this statement may not hold for systems with large uncertainty. Before proceeding with the integration of the areas, a visual inspection of the validity of these bounds is therefore advised (see Fig. 3.1).
Example. Bounds on tr(PjQi)
The following model will be used for a robust analysis of interactions based on tr(PjQi).
G(s) = 5 0.5s+1 2.4 0.25s+1 2.7 0.4s+1 2.5 0.25s+1 (3.3)
where each of the input-output channels is perturbed with multiplicative uncertainty
Wij= 0.15.
The nominal value of the unnormed PM ˜φN is:
˜
φN = 6.2500 1.4400
1.8225 1.5625
The regions in the Nyquist diagram induced by uncertainty are depicted in Fig. 3.2, from
−1 0 1 2 3 4 5 6 −3 −2 −1 0 1 2 3
Nyquist Diagram From In(1) to Out(1)
Real Axis Imaginary Axis −1 0 1 2 3 4 5 6 −3 −2 −1 0 1 2 3
Nyquist Diagram From In(2) to Out(1)
Real Axis Imaginary Axis −1 0 1 2 3 4 5 6 −3 −2 −1 0 1 2 3
Nyquist Diagram From In(1) to Out(2)
Real Axis Imaginary Axis −1 0 1 2 3 4 5 6 −3 −2 −1 0 1 2 3
Nyquist Diagram From In(2) to Out(2)
Real Axis
Imaginary Axis
Figure 3.2: Uncertainty in the Nyquist diagram of the system in Equation (3.3) when each input-output channel is affected by a multiplicative uncertainty of Wlij = 0.15.
which, the bounds on tr(PjQi) for each of the scalar channels with input j and output i
are computed and collected in the matrix ˜φ.
˜
φ = [4.1435, 8.9101] [0.9549, 2.0533]
[1.2123, 2.5980] [1.0359, 2.2286]
The sum of the diagonal elements of ˜φ is, in all cases, larger than the sum of the
off-diagonal ones, and therefore this indicator suggests the off-diagonal pairing in case of using
decentralized control. However, if ˜φ is normalized to obtain the PM, the sum of the
main diagonal will vary between 0.53 and 0.84, which means that such a decentralized structure will consider between 53% and 84% of the process dynamics quantified using tr(PjQi) as indicator. An upper triangular controller will be the result of discriminating
only the subsystem G12, and it is the simpler structure which will most likely derive in
30
Bounds on Interaction Measures from parametric models with uncertainty description dynamics. However, in the worst case for such a triangular structure, the significance of the neglected input-output channel will be 24%, which is rather large considering that in a system with four equal input-output channels the contribution of each channel is 25%. The final robust decision to be taken is that a full MIMO controller should be used,
since no input-output can be discriminated considering the variations of tr(PjQi) in the
uncertainty set.
3.1.1
Computing the area enclosed by the Nyquist Diagram.
Even if Equation (3.2) gives an analytical description of the curves enclosing the minimum and maximum areas, it might not be an easy task to compute this areas. Two methods are here proposed and discussed.
Computing the area enclosed by the Nyquist Diagram using triangulation The easiest way of computing this areas is by discretizing the curves at a finite num-ber of points, and using triangulation techniques to compute the area enclosed by the
resulting polygon. Therefore, for a chosen set of increasing positive frequencies ωk =
{ω0, ω1, . . . , ωN} the area enclosed by the resulting polygon can be computed as:
Area = (real(G(jωN)) · imag(G(jω0)) − real(G(jωN)) · imag(G(jωN))) +
+ N −1 X
k=0
real(G(jωk)) · imag(G(jωk+1)) − real(G(jωk+1)) · imag(G(jωk))
(3.4)
where G can be substituted by Ghighor Glowfor computing the maximum and minimum
areas induced by uncertainty. The triangulation technique gives the desired area, since it also considers the multiplicities due to several windings (positive or negative). Equa-tion (3.4) gives the total area enclosed by the Nyquist diagram using a set of positive
frequencies by making use of the existing symmetry. It is recommended to choose ω0= 0
and ωN → ∞. It is also recommended to choose the set of frequencies as logarithmically
spaced.
This method is simple to implement and is fast to compute, and many software tools have it in a built in function, like the function polyarea in Matlab. Another advantage of the triangular integration is that it is suitable in the case of having a nonparametric model of the frequency response function, which will be discussed in Section 4.2.1.
A clear disadvantage of this method is that it requires a certain interaction from the user, for selecting the set of frequencies at which G(jω) is evaluated, and the exactitude of the result depends on the ability with which this set of frequencies is selected. Computing the area enclosed by the Nyquist Diagram using numerical inte-gration
Equations 3.2 give the analytical description of the curves enclosing the minimum and maximum areas in the Nyquist diagram due to model uncertainty. These equations are
parametrizations of curves in the complex plane with respect to the parameter ω. The area can therefore be analytically solved by evaluating the following line integral in the complex plane: Area(γ) =1 2 I γ xdy − ydx
where γ denotes the considered Nyquist curve. However, depending on the number of zeros and poles of the original transfer function and on the complexity of the functions W (jω) used to describe uncertainty, this integral can be tedious to solve, and therefore, an alternative is here presented. The enclosed area can be computed as:
Area(γ) = Z ∞ −∞ Z ∞ −∞ nυ(x + j · y)dxdy
where nυ(x + j · y) is the number of times that γ winds around the point in the complex plane with real part x and imaginary part y. This integral can numerically be solved using quadrature methods (i.e. the methods implemented in Matlab in the function dblquad ), if an algorithm which computes the winding number of a point in the complex plane is provided to the quadrature method. Such algorithm was provided in [19] as an application of general concepts in geometry applied to the Nyquist diagram. From the results in that paper, the following algorithm can be formulated for computing the
winding number for a given point z0= x0+ j · y0in the complex domain:
1.- Define the horizontal half-line which starts at the point z0 and goes in the
positive direction of the abscissa. Denote it as %.
2.- Determine the number and the sign of the crossings of γ and %. Being a crossing
considered positive if the cross product of the tangent vectors to the curves t%× tγ
is also positive.
3.- The winding number is then the number of positive crossings minus the number of negative crossings.
There are two specials cases in this algorithm which have to be clarified. For ω → ∞ the Nyquist plot approaches to a point in the real axis, and the curve is not smooth. The same situation can occur at ω = 0, and therefore, the tangent to the Nyquist curve is not well-defined. A simple approach for determining the sign of the crossings in this cases can be consulted in [19].
3.2
Model uncertainty in the Bode diagram and
con-nection to the bounds of the H
2norm
The Bode diagram of a continuous-time dynamic system is obtained by plotting inde-pendently the magnitude and the phase of the complex number G(jω) as a function of ω. Only the magnitude of the Bode diagram will be analyzed, due to its connection with
32
Bounds on Interaction Measures from parametric models with uncertainty description Considering a MIMO system in which each of the SISO subsystems is affected by independent multiplicative uncertainty, the bounds on |G(jω)| are described for each of the SISO subsystems by:
|G(jω)|min= |G(jω)|(1 − |W (jω)|)
|G(jω)|max= |G(jω)|(1 + |W (jω)|)
As described in Equation (2.4), the H2norm can be computed as the area under the
squared magnitude of the Bode diagram |G(jω)|2. Since |G(jω)|
min and |G(jω)|max in
Equation (3.2) give the minimum and maximum values of |G(jω)|, substituting |G(jω)|
in Equation (2.4) for them gives the maximum and minimum values of ||G(s)||2 due to
uncertainty.
Example. Bounds on ˜Σ2
The model in Equation (3.3) will be used for a robust analysis of interactions based on ˜Σ2.
The nominal value of the unnormed Σ2 will be denoted as ˜ΣN
2 and has been calculated
to be:
˜
ΣN2 = 5.0000 3.3941
3.0187 3.5355
Each of the input-output channels is perturbed with multiplicative uncertainty Wlij =
0.15. The uncertainty areas in the magnitude of the Bode diagram induced by uncertainty are depicted in Fig. 3.3,
By analytically evaluating the expression in for |G(jω)|min and |G(jω)|max in , the
bounds on ||Gij||2for each of the scalar channels with input j and output i are computed
and collected in the matrix ˜Σ2.
˜
Σ2=
[4.2500, 5.7500] [2.8850, 3.9032] [2.5659, 3.4715] [3.0052, 4.0659]
This indicator suggests that no robust decision on a decentralized control structure can be taken, since the dominance of the sum of the diagonal elements or the off-diagonal
ones switches within the uncertainty set. The channel G21, which is the less significant
one in the nominal case, can contribute up to 26% of the process dynamics when the
H2 norm is used as quantification. The final conclusion is that no input-output channel
should be discriminated if all the uncertainty set is considered, and therefore a full MIMO controller is to be used. This conclusion is in accordance with the one in Section 3.1 where
100 102 −20 −10 0 10 ω (rad/sec) |G 11 (j ω )| 100 102 −20 −10 0 10 ω (rad/sec) |G 21 (j ω )| 100 102 −20 −10 0 10 ω (rad/sec) |G 12 (j ω )| 100 102 −20 −10 0 10 ω (rad/sec) |G 22 (j ω )|
Figure 3.3: Uncertainty in the magnitude of the Bode diagram of the system in Equation (3.3) when each input-output channel is affected by a multiplicative uncertainty of Wij= 0.15.
34
Bounds on Interaction Measures from parametric models with uncertainty description
Confidence bounds on Interaction
Measures from estimated
nonparametric models in the
frequency domain. Applications on
weakly non-linear systems
There is a possibility of estimating a non-parametric model of the FRF of a linear system by exciting the process inputs with for example white noise or with periodic signals [3]. In the case of a SISO system, the quotient between the Fourier transform of the process output and the Fourier transform of the process input is the Maximum Likelihood
(ML) estimation GM L(jωk) of the real FRF G(jωk) at the considered frequencies ωk.
Confidence bounds on the estimation of each of the values G(jωk) can be created as a
circular complex region determined by the variance of the estimators.
An example of the resulting model from such a modeling scheme is depicted in Fig. 4.4,
were the GM L(jωj) is depicted in a Nyquist diagram, and it becomes clear that the method
described in Section 3.1 can be used to create a robust estimation of the squared HS norm by computing the maximum and minimum areas under the Nyquist diagram using the triangulation method described in Section 3.1.1.
The magnitude of the estimated FRF can also be depicted in the frequency domain, and confidence intervals for the estimation can be created from the variance of the
com-puted estimators GM L(jωk). This will result in a diagram like the depicted in Fig. 4.6.
The bounds on the H2 norm can then be computed from the maximum and minimum
obtained areas under the squared magnitude of the Bode diagram in a similar approach that the one described in Section 3.2.
It is important to clarify that the contribution of this chapter is not proposing any modeling scheme or excitation design. The contribution is to present and illustrate a
methodology for creating bounds on the estimation of the H2and the squared HS norms
36
Confidence bounds on Interaction Measures from estimated nonparametric models in the frequency domain of the SISO subsystems of a MIMO system from nonparametric models describing the its FRF. The optimal choice of the modeling scheme and the excitation signals is to be decided by a qualified engineer depending on the process to be modeled and the scenario. However, a modeling scheme using the ML estimator for linear multivariable systems is included in Section 4.1.1.
In the case of facing a nonlinear system, Section 4.1.2 includes the description of a excitation design which will quantify the contribution of the process nonlinearities [3]. If the system is found to be weakly nonlinear, then the designed experiment will provide with data for obtaining the best linear approximation ([3]) of the FRF. The variance of the estimators can then be used to create confidence regions in the Nyquist and Bode
diagrams, and therefore it will be possible to compute a robust estimation of the H2and
the squared HS norms in sections 4.2.2 and 4.2.1 respectively. If the process is found to be strongly nonlinear, then an Interaction Measure for nonlinear systems have to be used (see Chapter I.1).
4.1
Estimation of the FRF of the system
4.1.1
Estimation of the FRF for linear multivariable systems.
Periodic excitation signals will be considered. To estimate G(jωk) for a multivariable sys-tem with n inputs and m outputs, n sub-experiments are needed. The Fourier transform of the inputs for each frequency will then be collected in a matrix:
U(k) = U1[1](k) . . . U [n] 1 (k) .. . . .. ... Un[1](k) . . . Un[n](k)
where Uj[i](k) is the frequency content of the jthinput at the frequency ω
k and in the ith sub-experiment. By collecting in the same way the Fourier transform of the output in a matrix such as Y(k) = Y1[1](k) . . . Y1[n](k) .. . . .. ... Yn[1](k) . . . Yn[n](k)
then the ML estimator GM L(jωk) at the frequency ωk can be computed as:
GM L(jωk) = Y(k)U−1(k) (4.1)
For guaranteing that U(k) is regular and well conditioned, the DFT spectrum of one of the input signals will be tailored, and U(k) will be selected as U(k) = U (k)W, being W and orthogonal matrix. I.e. for the case of a 2 × 2 system:
W = 1 1
1 −1
These orthogonal input signals for TITO systems were introduced in [20] and later con-sidered for an arbitrary input dimension in [21]. They have been designed for attenuating the influence of process noise in the FRF measurements.
In case of having at disposal several successive periods of the signals, different aver-aging techniques are possible (see [3]), being one of them to average the DFT spectrum over the successive periods.
4.1.2
Excitation for nonlinear systems.
Consider a Wiener system with an output nonlinearity y = a · z + b · z2+ c · z3 excited
with a sine wave with a frequency f0. In the case of a 6= 0, b = 0 and c = 0 the system is merely a linear system, and the output spectra has a contribution only at the original frequency f0. On the other hand, nonlinear systems would create additional harmonics at
frequencies different than f0. A quadratic term (b 6= 0) will create a contribution at 2f0
(and a DC term), and a cubic term (c 6= 0) will create a contribution at f0(superposed
to the linear contribution) and another contribution at 3f0.
The previous observations can be generalized using Volterra systems for a large class of nonlinear systems [22, 3]. The implication is that, for a nonlinear system, and for periodic signals exciting only a set of selected odd frequencies (i.e. ωk = f0, 3f0, 7f0, 11f0, . . . ), the contribution of the even nonlinearities will appear at even frequencies and therefore they
will not disturb the estimation of G(jωk), whereas the odd nonlinearities will contribute
at the odd frequencies which are not excited (5f0, 9f0, . . . ), but also contribute at the excited frequencies in superposition with the linear contribution.
From the obtained data, the best linear approximation GBLA(jωk) can be estimated
using the ML estimator, obtaining as result:
GBLA(jωk) = G0(jωk) + GB(jωk)
being G0(jωk) the underlying linear system and GB(jωk) the bias errors due to nonlinear contributions, which depends only on the odd degree nonlinear distortions and the power spectrum of the input [23, 3].
If for each period of the excitation signal the FRF is estimated, and the result is averaged, then the total variance of the linear estimator will be a composition of the variance induced by the additive output noise and a zero mean stochastic nonlinear contribution from the odd degree nonlinearities [23, 3].
Example. Excitation and FRF estimation of a weakly nonlinear system The quadruple tank system was described in [24] and is depicted in Fig. 4.1. It is an
interacting systems in which two pumps u1and u2deliver their flow in four tanks, being
the output flow of the top tanks delivered to their corresponding tanks at the bottom. This interconnections form a complex system whose dynamic behavior depends on the position of two manual valves which determine how the flow rate from each pump is split into two different tanks. The nonlinear differential equations describing the process
