# Multilinear Model of Heat Exchanger with Hammerstein Structure

Authors: Dragan Pršić; Novak Nedić; Vojislav Filipović; Ljubiša Dubonjić; Aleksandar Vičovac
Journal: Journal of Control Science and Engineering
 (2016)

---

## Abstract

The multilinear model control design approach is based on the approximation of the nonlinear model of the system by a set of linear models. The paper presents the method of creation of a bank of linear models of the two-pass shell and tube heat exchanger. The nonlinear model is assumed to have a Hammerstein structure. The set of linear models is formed by decomposition of the nonlinear steady-state characteristic by using the modified Included Angle Dividing method. Two modifications of this method are proposed. The first one refers to the addition to the algorithm for decomposition, which reduces the number of linear segments. The second one refers to determination of the threshold value. The dependence between decomposition of the nonlinear characteristic and the linear dynamics of the closed-loop system is established. The decoupling process is more formal and it can be easily implemented by using software tools. Due to its simplicity, the method is particularly suitable in complex systems, such as heat exchanger networks.


---

## Body

## 1. Introduction

Most physical systems are inherently nonlinear. Nonlinearity is especially pronounced in systems with wide ranges of desired behaviours and variable set points. There are a lot of such systems in the process industry [1, 2]: pH reactors, distillation columns, polymerization reactors, heat exchangers, and so forth. In order to solve the problem of design and high precision tracking control of such systems, it is necessary to have a better understanding of their nonlinear characteristics.Linear control theory provides many confirmed methods and tools for controller design with desired performances and robustness. Unfortunately, that theory is limited to strictly linear systems or certain classes of nonlinear systems with small deviations around the nominal operating point. In real systems with a wide operating range and multiple operating points, where nonlinearity cannot be ignored, different control strategies are necessary. In order to use the vast potential of linear control theory, at the same time taking into account nonlinearities, various modifications of classical control design are proposed.The multilinear model (or multimodal) control approach which has been shown to be suitable for strongly nonlinear systems with multiple operating points, tracking control, and wide operating ranges [3–7] has been drawing attention recently. The main idea is to present the nonlinear system as a set of linear systems where classical controller design techniques can be easily applied. The multilinear model control design framework includes four steps: determination of the global (nonlinear) mathematical model of the plant for the whole operating region; approximation of the nonlinear model by means of a set of linear models; local controller design for each particular operating region; global controller design of the whole nonlinear system.There are two important questions which should be answered in the multilinear model control design approach: how to decouple a nonlinear system into a bank of locally linear subsystems and how to design a global controller according to the desired performances of a nonlinear system.This paper gives the answer to the first question. The literature mainly uses the gap-metric based method where the minimal linear model set is determined for the given threshold value in order to span the expected operating range of a nonlinear system [4–7]. For a special class of nonlinear systems, such as nonlinear systems with Hammerstein-like structures, the authors in [8, 9] used the Included Angle Dividing (IAD) method for creating a minimum set of models. Numerous nonlinear heating, ventilating, and air-conditioning (HVAC) systems have a nonlinear static block followed by linear dynamics [10]. The paper assumes that the model of the two-pass shell and tube heat exchanger (HE) has a Hammerstein structure and a modified method of IAD is used. The major modification is seen in establishing the threshold value depending on the linear dynamics of the closed-loop system.


## 2. Model of Two-Pass Shell and Tube Heat Exchanger

### 2.1. Nonlinear Model

Nowadays, a lot of dynamic models defined to study transient responses of HE [11–14] can be encountered. This paper uses a simplified model of heat exchanger (HE) with lumped parameters. There are two reasons for that: (a) from the aspect of control, simplified models are allowed because control algorithms allow uncertainty in models; (b) HE are commonly a part of more complex systems (HE networks), so that suitable but simple models are desirable.The model of HE development by Chen et al. [12] is used in this paper. The two-pass shell and tube HE with the process water flow rate as the input and the process water exit temperature as the output is analysed (Figure 1). The output temperature also depends on the steam flow rate; however, that dependence is approximately linear. The results showed that the changes in the process flow rate have a much bigger influence on the nonlinear behaviour of HE. Therefore, it is assumed that the steam flow rate has a constant value. In addition, flow controllers are used during identification to avoid the nonlinearity of the valves. The main cause of nonlinear behaviour is the change of the heat transfer area between the process water and the steam. Namely, the amount of condensed steam which floods the exchanger varies with the change of the input flow rate. At higher flow rates, more steam is converted into liquid, so that the heat transfer area decreases, and vice versa.Figure 1
Functional view of two-pass shell and tube heat exchanger [11].It is assumed that the HE has a Hammerstein structure (Figure2). The models with this structure have a wide application in the process industry [8, 9, 11, 12]. They are block-oriented models which consist of a series connection of a nonlinear static element N(·) with a linear dynamic element Gp(s). The advantage is that the existing theory of linear systems can easily be adjusted and used for design of nonlinear systems with this structure.Figure 2
Hammerstein structure of heat exchanger model.Iterative identification with the pseudorandom binary sequence (PRBS) input in the closed loop with the switching time of 1 min is used for determination of the model [12]. The results of simulation and experiment showed a good correlation with the actual response data in a wide operating range.The collected data were the basis for obtaining the model which, after translation into the continuous domain (sampling time 12 s, zero-order hold on the input), has the following form:(1a)Gps=YsXs=-0.7825s+1154.1s2+54.84s+1,(1b)xu=Nu=-33.30u+44.47u2+48.01u3,where u denotes the changes in the set point of the process water flow controller and y denotes the change in the process water exit temperature.Model (1a)-(1b) is determined for the following operating point [11]: the steam flow rate (62% of the maximum), the water flow rate (42% of the maximum), and the inlet water temperature (30°C).Figure3 presents the estimated static nonlinearity and the steady-state operating points.Figure 3
Estimated static nonlinearity and steady-state operation points of heat exchanger [12].


### 2.2. Multilinear Model

In the next stage of multilinear model control design approach, the goal is to approximate the nonlinear plant (1a)-(1b) by a set M of piecewise linear approximations Mi:(2a)M=Mi,i=1,2,3,…,n,where Mi is the ith linearized model within a specified operating range:(2b)Mi:KHEiGps=KHEi-0.7825s+1154.1s2+54.84s+1,i=1,…,n,(2c)KHEmin≤KHEi≤KHEmax,where KHEmin and KHEmax are the minimum and maximum values of gain of the heat exchanger, respectively.The change of gain of the heat exchanger as a function of the control signalu is presented in Figure 4.Figure 4
The change of gain of the heat exchanger.The main question which should be answered is the following: how many and which linear models are required to span the expected operating region of a nonlinear system? An approach based on the IAD method [8, 9] is used in this paper for decomposition of the nonlinear model. The reason is its simplicity. The method is primarily used in SISO nonlinear systems with Hammerstein and Wiener structures. It is particularly appropriate for systems with input multiplicity and noninvertible nonlinearity. It is assumed that the nonlinear steady-state characteristic is known is advance, either from identification or from the first-principle modelling. Creation of the bank of linear models is based on decoupling the nonlinear characteristic into linear segments. The algorithm of the method is here briefly stated for completeness of the presentation (Algorithm 1), where γa is the threshold value, Na is the number of steady-state points, Pi is the ith operating point, uo,xoi are coordinates of ith operating point, θi is the slope angle of ith operating point, θij is the difference between two slope angles at two operating points, and Si is the ith linearized segment.Algorithm 1:Algorithmic view of the Included Angle Dividing method.
(1) Prescribeγa, Na (2) DistributeNa steady-state points Pi, i=1,2,…,Na along the steady-state IO curve evenly (3) Computeki, θi, i=1,2,…,Na (4) fori=1:Na (5)   forj=i:Na (6)    ifθij<γa (7)Sij=Pj (8)    else (9)j=Na (10)    end (11)   end (12) end (13) compression subregion of a small range (14) In the middle ofSi choose operating point uo,xoiAfter the previous procedure of 14 steps, a set of linearized segments and the corresponding operating points,Si, uo,xo,yoi, i=1,2,…,Nm, respectively, are obtained. Also,(3)⋃i=1NmSi=ORholds, where OR is the expected operating region of the nonlinear system.In order to obtain a set of linear models, it is necessary to identify the local linear model or linearize the nonlinear first-principle model for each segment, around its operating point.The method is simple and easily realized. It is only necessary to know the nonlinearity of the steady-state characteristic.The paper proposes two modifications of the IAD method. The first modification refers to the algorithm shown in Algorithm1. As it can be seen, the inner loop of the algorithm (lines 5–11) tests the closeness (line 6) of the current operating point, along the main diagonal, and the points to the right of the angle-matrix θijNa×Na. Our proposal is to introduce another inner loop which would test the closeness of the current operating point and the points on the left side on the same row of the angle-matrix (Algorithm 2: lines 12–18).Algorithm 2:Algorithmic view of the modified IAD method.
(1) Prescribeγa, Na (2) DistributeNa steady-state points Pi, i=1,2,…,Na along the steady-state IO curve evenly (3) Computeki, θi, i=1,2,…,Na (4) fori=1:Na (5)   forj=i:Na (6)    ifθij<γa (7)Sij=Pj (8)    else (9)j=Na (10)    end (11)   end (12)   forj=i:1:step  −1 (13)    ifθij<γa (14)Sij=Pj (15)    else (16)j=1 (17)    end (18)   end (19) end (20) compression subregion of a small range (21) In the middle ofSi choose operating point uo,xoiThe header of the new loop has the form fori=i:1:step  −1. The segments Si which include a larger number of operating points are thus obtained.The second modification refers to two weaknesses of the IAD method:(i)
Selection of the threshold value requires certain experience and a priori knowledge. Therefore, it is difficult to make a procedure (e.g., software tool) for a systemized approach.(ii)
Decomposition of the plant model is performed only on the basis of the nonlinear steady-state characteristic; that is, the dynamics of the linear part are not taken into account.The starting point in solving these problems is the desired dynamics of closed-loop systems in order to make the basis for defining the criterion for segmentation of the nonlinear steady-state characteristic. In other words, the threshold valueγa of the nonlinear plant is defined based on the desired behaviour of the closed-loop system.In [8, 9], the threshold value is defined depending on the plant (HE) gain. However, it is known that the closeness of the two open-loop processes does not guarantee the closeness of their closed-loop systems [15, 16]. The step responses (u(t)=-0.02h(t)) for four different values of linearized gain of plant (1a)-(1b) are shown in Figure 5. It could be noticed that all four responses have basically the same character.Figure 5
Step response of the HE for different values of linearized gains.On the other hand, step responses (u(t)=h(t)) of the closed-loop system with the unit gain error feedback for the same values of the linearized gain of plant (1a)-(1b) are shown in Figure 6. It can be seen that in the closed loop the character of the transitional process changes depending on the gain of plant (1a)-(1b).Figure 6
Step response of the closed-loop system with unit gain.


## 2.1. Nonlinear Model

Nowadays, a lot of dynamic models defined to study transient responses of HE [11–14] can be encountered. This paper uses a simplified model of heat exchanger (HE) with lumped parameters. There are two reasons for that: (a) from the aspect of control, simplified models are allowed because control algorithms allow uncertainty in models; (b) HE are commonly a part of more complex systems (HE networks), so that suitable but simple models are desirable.The model of HE development by Chen et al. [12] is used in this paper. The two-pass shell and tube HE with the process water flow rate as the input and the process water exit temperature as the output is analysed (Figure 1). The output temperature also depends on the steam flow rate; however, that dependence is approximately linear. The results showed that the changes in the process flow rate have a much bigger influence on the nonlinear behaviour of HE. Therefore, it is assumed that the steam flow rate has a constant value. In addition, flow controllers are used during identification to avoid the nonlinearity of the valves. The main cause of nonlinear behaviour is the change of the heat transfer area between the process water and the steam. Namely, the amount of condensed steam which floods the exchanger varies with the change of the input flow rate. At higher flow rates, more steam is converted into liquid, so that the heat transfer area decreases, and vice versa.Figure 1
Functional view of two-pass shell and tube heat exchanger [11].It is assumed that the HE has a Hammerstein structure (Figure2). The models with this structure have a wide application in the process industry [8, 9, 11, 12]. They are block-oriented models which consist of a series connection of a nonlinear static element N(·) with a linear dynamic element Gp(s). The advantage is that the existing theory of linear systems can easily be adjusted and used for design of nonlinear systems with this structure.Figure 2
Hammerstein structure of heat exchanger model.Iterative identification with the pseudorandom binary sequence (PRBS) input in the closed loop with the switching time of 1 min is used for determination of the model [12]. The results of simulation and experiment showed a good correlation with the actual response data in a wide operating range.The collected data were the basis for obtaining the model which, after translation into the continuous domain (sampling time 12 s, zero-order hold on the input), has the following form:(1a)Gps=YsXs=-0.7825s+1154.1s2+54.84s+1,(1b)xu=Nu=-33.30u+44.47u2+48.01u3,where u denotes the changes in the set point of the process water flow controller and y denotes the change in the process water exit temperature.Model (1a)-(1b) is determined for the following operating point [11]: the steam flow rate (62% of the maximum), the water flow rate (42% of the maximum), and the inlet water temperature (30°C).Figure3 presents the estimated static nonlinearity and the steady-state operating points.Figure 3
Estimated static nonlinearity and steady-state operation points of heat exchanger [12].


## 2.2. Multilinear Model

In the next stage of multilinear model control design approach, the goal is to approximate the nonlinear plant (1a)-(1b) by a set M of piecewise linear approximations Mi:(2a)M=Mi,i=1,2,3,…,n,where Mi is the ith linearized model within a specified operating range:(2b)Mi:KHEiGps=KHEi-0.7825s+1154.1s2+54.84s+1,i=1,…,n,(2c)KHEmin≤KHEi≤KHEmax,where KHEmin and KHEmax are the minimum and maximum values of gain of the heat exchanger, respectively.The change of gain of the heat exchanger as a function of the control signalu is presented in Figure 4.Figure 4
The change of gain of the heat exchanger.The main question which should be answered is the following: how many and which linear models are required to span the expected operating region of a nonlinear system? An approach based on the IAD method [8, 9] is used in this paper for decomposition of the nonlinear model. The reason is its simplicity. The method is primarily used in SISO nonlinear systems with Hammerstein and Wiener structures. It is particularly appropriate for systems with input multiplicity and noninvertible nonlinearity. It is assumed that the nonlinear steady-state characteristic is known is advance, either from identification or from the first-principle modelling. Creation of the bank of linear models is based on decoupling the nonlinear characteristic into linear segments. The algorithm of the method is here briefly stated for completeness of the presentation (Algorithm 1), where γa is the threshold value, Na is the number of steady-state points, Pi is the ith operating point, uo,xoi are coordinates of ith operating point, θi is the slope angle of ith operating point, θij is the difference between two slope angles at two operating points, and Si is the ith linearized segment.Algorithm 1:Algorithmic view of the Included Angle Dividing method.
(1) Prescribeγa, Na (2) DistributeNa steady-state points Pi, i=1,2,…,Na along the steady-state IO curve evenly (3) Computeki, θi, i=1,2,…,Na (4) fori=1:Na (5)   forj=i:Na (6)    ifθij<γa (7)Sij=Pj (8)    else (9)j=Na (10)    end (11)   end (12) end (13) compression subregion of a small range (14) In the middle ofSi choose operating point uo,xoiAfter the previous procedure of 14 steps, a set of linearized segments and the corresponding operating points,Si, uo,xo,yoi, i=1,2,…,Nm, respectively, are obtained. Also,(3)⋃i=1NmSi=ORholds, where OR is the expected operating region of the nonlinear system.In order to obtain a set of linear models, it is necessary to identify the local linear model or linearize the nonlinear first-principle model for each segment, around its operating point.The method is simple and easily realized. It is only necessary to know the nonlinearity of the steady-state characteristic.The paper proposes two modifications of the IAD method. The first modification refers to the algorithm shown in Algorithm1. As it can be seen, the inner loop of the algorithm (lines 5–11) tests the closeness (line 6) of the current operating point, along the main diagonal, and the points to the right of the angle-matrix θijNa×Na. Our proposal is to introduce another inner loop which would test the closeness of the current operating point and the points on the left side on the same row of the angle-matrix (Algorithm 2: lines 12–18).Algorithm 2:Algorithmic view of the modified IAD method.
(1) Prescribeγa, Na (2) DistributeNa steady-state points Pi, i=1,2,…,Na along the steady-state IO curve evenly (3) Computeki, θi, i=1,2,…,Na (4) fori=1:Na (5)   forj=i:Na (6)    ifθij<γa (7)Sij=Pj (8)    else (9)j=Na (10)    end (11)   end (12)   forj=i:1:step  −1 (13)    ifθij<γa (14)Sij=Pj (15)    else (16)j=1 (17)    end (18)   end (19) end (20) compression subregion of a small range (21) In the middle ofSi choose operating point uo,xoiThe header of the new loop has the form fori=i:1:step  −1. The segments Si which include a larger number of operating points are thus obtained.The second modification refers to two weaknesses of the IAD method:(i)
Selection of the threshold value requires certain experience and a priori knowledge. Therefore, it is difficult to make a procedure (e.g., software tool) for a systemized approach.(ii)
Decomposition of the plant model is performed only on the basis of the nonlinear steady-state characteristic; that is, the dynamics of the linear part are not taken into account.The starting point in solving these problems is the desired dynamics of closed-loop systems in order to make the basis for defining the criterion for segmentation of the nonlinear steady-state characteristic. In other words, the threshold valueγa of the nonlinear plant is defined based on the desired behaviour of the closed-loop system.In [8, 9], the threshold value is defined depending on the plant (HE) gain. However, it is known that the closeness of the two open-loop processes does not guarantee the closeness of their closed-loop systems [15, 16]. The step responses (u(t)=-0.02h(t)) for four different values of linearized gain of plant (1a)-(1b) are shown in Figure 5. It could be noticed that all four responses have basically the same character.Figure 5
Step response of the HE for different values of linearized gains.On the other hand, step responses (u(t)=h(t)) of the closed-loop system with the unit gain error feedback for the same values of the linearized gain of plant (1a)-(1b) are shown in Figure 6. It can be seen that in the closed loop the character of the transitional process changes depending on the gain of plant (1a)-(1b).Figure 6
Step response of the closed-loop system with unit gain.


## 3. Determination of the Threshold Value

Let us observe the multilinear model based control system shown in Figure7. The control is realized at two levels: the first level uses the local controller with the fixed linear block Gc(s) and the variable segment gain (1/Kc).Figure 7
Multilinear model based control system.At the second level, the local controllers are switched by means of the global controller (GC) using hard-switching or soft-switching methods [17, 18]. The idea is to control the linear dynamics of the plant (Gp(s)) by the linear part of the controller (Gc(s)), while the output gain (1/Kc) is adjusted according to the linearized parts of static nonlinearity (N·).It should be noted thatKc belongs to a set of discrete values(4)Kc∈KHEi,i=1,2,3,…,n,where n is the number of linear models.The linear modelsMi differ only by the gain value KHEi:(5)MiKHEi=KHEiGps.SinceGc(s) is projected based on Gp(s), in order to have the desired behaviour of the closed-loop system, the following should hold:(6)1-δ-≤KHEKc≤1+δ+,where KHE is the current value of the heat exchanger gain (Figure 4).The parametersδ- and δ+ represent the allowed tolerance in the change of the open-loop gain of the linearized system L(s)=Gc(s)Gp(s) for the case when KHE/Kc<1 and KHE/Kc>1, respectively. They can be determined based on the desired behaviour of the closed loop, where (7)0≤δ-≤1;0≤δ+holds.If the value(8)ΔKc=KHE-Kcwhich represents the difference between the current value of gain of the nonlinear plant and the current local value of the controller gain is introduced, then, based on (6), it can be written that(9)-δ-≤ΔKcKc≤δ+.Inequality (9) defines the relative change of gain of the plant relative to the gain of the controller from the aspect of behaviour of the closed loop.If the current gain is now defined by using the slope angle of the steady-state characteristic (φ) (Figure 3), (10a)Kc=tan⁡φ,that is, when Δφ→0,(10b)ΔKcΔφ=1cos2⁡φ,then, based on (9), it can be written that(11)-δ-≤2Δφsin⁡2φ≤δ+.Relation (11) defines the allowed tolerance in the slope of the steady-state characteristic of the Hammerstein system. The characteristics of the closed system (δ-,δ+) and the static nonlinearity (φ) are taken into account. In [8], that tolerance has a constant value along the steady-state characteristic.If(12)δ-=δ+=δ,then relation (11) can be written in the following form:(13)2Δφsin⁡2φ≤δ.In the case of heat exchanger with the steady-state characteristic shown in Figure 3, the change of the slope angle of the steady-state characteristic is shown in Figure 8.Figure 8
Slope angle of the HE steady-state characteristic.It can be seen from the figure that(14)sin⁡2φ<0∀u∈-0.4,0.3holds and then condition (13) for the HE can be written in the following form:(15)Δφ≤γδ,φ=12δsin⁡2φ,0≤δ≤1,where γδ,φ represents the threshold value in creating the set of multilinear models of HE.Figure9 shows the change of the value γ by several values of the parameter δ. It can be seen from the figure that the value γδ,φ has different values for acceptable changes of the open-loop gain depending on the control value u.Figure 9
Threshold value for different values ofu.In selecting values for the parameterδ, the following should be taken into account:(i)
Parameters of a linearized model depend on the operating conditions.(ii)
Hammerstein model is just an approximation of a real system. The parameterδ should guarantee robust stability in the presence of uncertainty in dynamical systems (parameter variations and neglected dynamics).(iii)
Transfer from one linear model to another (hard-switching or soft-switching) can result in oscillations in the system.(iv)
Linear dynamics of the system depend not only on the dynamics of the plant but also on the parameters of the local controller (Figure7).(v)
Value of the parameterδ can be estimated by simulation for a wide range of values of the parameters.TheD-decomposition and the parameter plane can be used for establishing the connection between the linear dynamics and the threshold value, as in [19]. The Nyquist plot with the stability margin (Sm) (Figure 10) as a measure of system robustness can be used for approximate estimation.Figure 10
Nyquist plot with the stability margin.Figure11 presents the change of the number of linearized segments depending on the allowed tolerance in the change of the open-loop gain (δ). The dependence is presented for different values of the number of steady-state points (Na).Figure 11
Number of segments depending onδ (parameter Na).The modified IAD method was used. It can be seen that in this method the number of linear models practically depends to a small extent on the number of steady-state points. That dependence is more pronounced in lower values of the parameterδ. For δ>0.5, that dependence can be neglected, so that the number of segments is determined by the allowed tolerance in the change of the open-loop gain.Figure12 presents the linearized segments of static nonlinearity of the heat exchanger for (a) 2 segments (δ=0.5) and (b) 5 segments (δ=0.25).Figure 12
Linearization of the heat exchanger static nonlinearity.
(a)
Two segments (b)
Five segmentsThe corresponding set points are shown in Tables1 and 2.Table 1
Operating points for 2 segments.
u o −0.09 0.21 y o −2.72 −4.61 K c −23.53 −7.84Table 2
Operating points for 5 segments.
u o −0.03 0.11 0.17 0.21 0.22 y o 0.93 −3.15 −4.25 −4.60 −4.68 K c −35.58 −21.23 −12.82 −8.16 −6.21


## 4. Conclusions

The paper presents the method for creating a set of linear models of heat exchanger described by the Hammerstein model type. Decomposition of the model is based on linearization of the nonlinear steady-state characteristic which is assumed to be a priori known. The method is intuitive and can be easily implemented by using software tools. It is shown that the number of linear segments practically depends only on the a priori given threshold value but that it does not depend on the number of steady-state points. The threshold value can be determined based on the allowed changes of the open-loop gain, which establishes the correlation between decomposition of the nonlinear model and the linear dynamics. The proposed method can easily be extended to other SISO Hammerstein-like systems with a memoryless nonlinearity.



---

*Source: 7856510-2016-05-25.xml*