Microsoft Word - 1_R.doc HUNGARIAN JOURNAL OF INDUSTRIAL CHEMISTRY VESZPRÉM Vol. 37(1). pp. 37-49 (2009) DYNAMIC MODEL AND CONTROL OF HEAT EXCHANGER NETWORKS FOR DISTRICT HEATING L. DOBOS 1, J. JÄSCHKE2, J. ABONYI1, S. SKOGESTAD2 1University of Pannonia, Department of Process Engineering, Egyetem str. 10 Veszprém, 8200, HUNGARY E-mail: dobosl@fmt.uni-pannon.hu 2Norwegian University of Science and Technology, Department of Chemical Engineering, NO-7491 ,Trondheim, NORWAY The various governmental policies aimed at reducing the dependence on fossil fuels for space heating and the reduction in its associated emission of greenhouse gases such as CO2 demands innovative measures. District heating systems using residual industrial waste heats could provide such an efficient method for house and space heating. In such systems, heat is produced and/or thermally upgraded in a central plant and then distributed to the final consumers through a pipeline network. In this work two main objectives will be considered: the first is to create a dynamic model which can represent the main characteristics of a district heating network and the second one is to design a non-linear model predictive controller (NLMPC) to satisfy the heat demands of the consumers in the heat exchanger network. As the model predictive controller is based on minimizing an objective function, it is totally perfect to find the way to reduce the superfluous energy consumption and make the best of using the freely applicable industrial waste heats. Beside this environmental aspect, reducing the invested energy consumption can reduce the operational costs. Keywords: district heating network, modeling, non-linear model predictive control, MPC Introduction It has become natural for people today to have a network for the distribution of the electricity. However the picture is much different when it comes to heating. The majority of the buildings in western Europe are heated with individual boilers that are fed either with natural gas or with oil. Only in some cases, e.g. with waste incineration, a district heating network (DHN) is implemented to distribute the heat. District heating networks are for distributing heat generated in a centralized location for residential and commercial heating requirements. The heat can be obtained from cogeneration plants or waste incineration plants although to satisfy the periodically increased heat demand so- called heat-only/peak load boiler stations are also used. These stations can be suppliers of residential and commercial consumers for space heating and for hot tap water and if necessary it can provide heat for industrial consumers for a certain level. However, because of the numerous advantages of district heating systems, it would be beneficial to implement them in other areas too, since the main advantages of district heating systems are [12]: 1. Fewer sources of emission in densely populated areas. 2. Fewer individual boilers, thus increased the usable space in the buildings. 3. Professional and on-going operating and maintenance of the centralized heating technology. Since energy becomes a more and more competitive market nowadays, thus optimization of the energy production and distribution becomes an important task for energy companies. The district heating facility can provide higher economic and environmental efficiency compared to local boilers. This the reason why the importance of these networks are increasing, and the countries, which use local heat suppliers, such as the Nordic countries, are switching to district heating networks instead. District heating networks exist in several variations: in the district heating network reported in [1] includes several consumers are located in different areas, but there is no energy storage and just one production unit. In [13], a storage tank is added to the network. In [6], a storage tank is also considered, but there is no thermal energy supply network. So the variety of the district heating networks are numerous. The operation of a district heating network is subject to operational constraints, e.g. assure the minimum inlet temperature of consumers this way satisfying their heat demand. The aim of the control strategies therefore is to meet these restrictions and at the same time to minimize the operational costs of the heat supplier. Model predictive control methods an interesting alternative to conventional control structures since the formulation of the objective function and constraints takes both aspect into consideration. Operating a district heating network implies to assign values to integer variables (status of production units, 38 status of pumps…) and to continuous variables (amounts of energy to produce). As a result, the optimization of the production and energy supply planning appears to be a huge, mixed integer and non linear optimization issue. Consequently, most studies use a simplified model, leaving aside some of the district heating network aspects. Simplifying the model may allow the use of one of the classical optimization methods listed in [11], but the solution can be strongly suboptimal when applied to the whole district heating network. Our goal is to model a district heating network presented latter and at the same time we use the model in a model predictive controller to satisfy the heat demand of the consumers of the network. The work is organized as follows: in the modeling section the topology of the applied district heating network will be introduced and the applied model equations will be defined. In the control section first the general non-linear model predictive controllers will be described, then defining the purposes, introducing the applied MPCs and examining the control results. Modeling section The models of a district heating network found in the literature are either a physical description of the heat and mass transfer in the network [8] or they are based on a statistical description of the transfer function from the supply point to the critical point considered. A statistical modeling approach is presented [10] where an ensemble of ARMAX (Auto-Regressive Moving Average with Exogenous input) models with different fixed time delays is set up, and depending on some estimated current time the models are switched. In [5] the grey-box approach for modeling combines physical knowledge with data-based (statistical) modeling; physical knowledge provides the main structure and statistical modeling provides details on structure and the actual coefficients/ estimates. The model in this work is developed with using the method of [7] so applying the physical description of the heat and mass transfer in the network. Local models of the components of the network are established and then connected together. Topology The topology showing Fig. 1 was chosen to represent the main characteristics of a district heating network. The network contains two heat production units, three consumers, two pumps and a valve. The production unit, called Producer 1, is the base load boiler, which may be considered a waste incineration plant. The other production unit, called Producer 2, is the peak load boiler station, which has to satisfy the increased heat demand in the network, especially in case of the Consumer 3. HX1 and HX2 heat exchangers are for transfer the produced heat from the primary circles to the secondary circle that is practically for distributing the heat for the consumers. Figure 1: The topology of the examined district heating network During the modeling procedure the following simplifications and assumptions were made to avoid the excessive complexity of modeled network, while preserving important characteristics: ● Since the system containes only pressurized water thermodynamical and material properties like heat capacities and densities are assumed constant. Average values for the respective temperature intervals are used. ● Isothermal flow is assumed through the pumps and valve. This is done due to the low pressure differences in the system. ● The pressure profile of the system changes much faster than the temperature profile, so it was modeled using steady-state equations, while the heat exchangers are modeled with dynamical assumption. The following model equations were applied to describe the network: Valves The valve is modeled using the following equation: 2 2vpp inout ⋅ ⋅−= ρξ (1) where: pout – outlet pressure of the valve pin – outlet pressure of the valve ρ – density of water v – velocity of the fluid ξ means the valve coefficient that is calculated by the expression below: 100/(%)ngvalveopeni )dtotalopene(ξξ = (2) 39 As it was mentioned previously, there is no difference in value of the inlet and outlet temperature. Pumps By neglecting any temperature rise in the water during travel through the pumps, they can be described using the Bernoulli equation. The elevation difference is set to zero, and the pipe diameter was assumed equal before and after the pump, thus the following expression characterizes the pump: m Ppp inout ηρ ⋅⋅ += (3) where P means the pump duty, η means the efficiency of the pump, m means the mass flow. Mixers The mixing unit is modeled using the following expressions, under the assumption of instant and homogeneous mixing: ∑ = = N i iout mm 1 (4) out N i ii out m Tm T ∑ = ⋅ = 1 (5) In modeling the mixers’ pressure, we assume: Pin, i = Pout (6) where: T – temperature P – inlet and outlet pressure. Pipes In modeling the district heating networks taking the effects of pipelines into consideration is an important factor. The heat loss in the pipes can not be neglected, and the dead time between the ends of the pipe has to be accounted for. Assuming one-dimensional flow, this leads to the partial equation: 0)),((2 ),()(),( 02 2 =− ⋅⋅ ⋅ + ⋅⋅ + TtxT Rc tx x T R tmtx t T p ρ μ ∂ ∂ ρπ∂ ∂ (7) where: cp – specific heat capacity of the fuid R – the radius of the pipe μ – heat transfer coefficient on the wall T0 – the ambient temperature. This equation has the following solution: ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − ⋅⋅ ⋅ − ⋅−−+= ))((2 000 0 )))((()( ttt Rc inout peTtttTTtT ρ μ (8) the varying time delay t – t0(t) is defined by: Ld R mt t = ⋅⋅∫ τ ρπ τ 0 2 )( (9) L is the length of the pipe (m). As the thermal losses on pipes are assumed very low. Thus the Eq. 8 is approximated by the following expression: ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − ⋅⋅ ⋅ −⋅−= ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − ⋅⋅ ⋅ −⋅−−+= ))((21))(( ))((21)))((()( 00 0000 ttt Rc tttT ttt Rc TtttTTtT p in p inout ρ μ ρ μ (10) In Eq. 10 constant time delay is assumed regarded to the high computation demand of calculating varying time delay. This approach yields a simple non-linear dynamic system, which can be quickly solved. The mechanical losses in pipes are modeled by: D Lvp 2 2⋅ ⋅= ρξΔ (11) Δp – pressure drop of the pipe D – diameter of the pipe. Heat exchangers In order model the proper dynamic behavior of the heat exchangers an approach using a cell model with ordinary differential equations is chosen [4]. This means that the heat exchanger is divided into perfectly and instantly mixed tanks. Each cell featuring a hot side, a wall side, and a cold side element (Fig. 2). The idea is that this approach will approximate the logarithmic mean temperature difference of the heat exchanger as the number of cells increases, while showing a realistic time delay behavior. In our model five cells were used on the hot side and five cells on the cold side. It is assumed that all cells are identical, and that no back-mixing occurres. In addition we assume that the mixing is instantaneous. Figure 2: Cell model of the heat exchanger 40 The following equations are applied to a cell: Hot side: ))()(())()1(( )( iTiTAUiTiTcV dt iTcdV chhhph hph −⋅⋅−−−⋅⋅⋅ = ⋅⋅⋅ • ρ ρ (12) Cold side: ))()(())()1(( )( iTiTAUiTiTcV dt iTcdV chhcpc cpc −⋅⋅+−+⋅⋅⋅ = ⋅⋅⋅ • ρ ρ (13) To avoid excessive complexity of the model the resistance of the wall is included to the heat transfer coefficient (U). The pressure drop of a heat exchanger is usually made up of the following parts: ● Pressure drop of the inlet nozzle. ● Pressure drop caused by the friction in the shell and in the tubes. ● Pressure drop of the outlet nozzle. To model these areas separately it would be necessary to use complex equations. To reduce the number of the expressions (and because we do not have a real DHS we can fit the model to) the pressure drop of a heat exchanger is approximated with the model equation of a valve. This way it was possible to model the pressure drop as the function of the flow rate and at the same time keep the model simple. Heat production units The approach for modeling the heat production units are similar to the model of the heat exchangers, however in this case just the cold side was divided into cells, following the scheme below: Figure 3: Cell model of the production unit The following equation is for representing the model of a cell on the cold side: N QiTiTcV dt iTcdV ccpc cpc +−−⋅⋅⋅= ⋅⋅⋅ • ))()1(( )( ρ ρ (14) where: Q – invested heat N – the number of the cells. This simplification is introduced because in the aspect of the heating network it is not important how the heat was produced, just the quantity of the invested heat is significant. The previously introduced network was implemented in Simulink. An other model of this network is implemented in Matlab with the simplification of neglecting the time delay of the pipelines (in Eq. 8). This was done to compare the performance of the two models in the MPC controller. Control section Model based control concepts The development of modern model based control concepts can be traced back to the early 1960s, to the linear quadratic regulator (LQR), which is designed to operate linear dynamic system. However, even though the concept of minimizing an objective function is very simple, the complexity of the controlled systems require advanced algorithms to solve the problem. At the early stage of model predictive controllers there were no sufficient computation facilities to realize complex optimization algorithms, so linear control algorithms were preferred [3, 9]. One example is Dynamic Matrix Controller (DMC), because of the analytical solution to the objective function. With growing computation power more complex methods could be applied to solve non- linear complex optimization problems in short time, and this way the non-linear model predictive algorithms could be born. Model Predictive Controllers– theoretical basis MPC is a model based control algorithm where models are used to predict the behavior of dependent variables (i.e. outputs) of a dynamical system with respect to changes in the process independent variables (i.e. inputs). In chemical processes, independent variables are often setpoints of regulatory controllers that govern valve movement (e.g. valve positioners with or without flow, temperature or pressure controller cascades), while dependent variables might be constraints in the process (e.g. product purity, equipment safe operating limits), however it is not neccessary (e.g. in some cases temperature is measured but it is not a constraint). The MPC uses the models and current plant measurements to calculate future moves in the independent variables that will result in operation that minimizes the cost function and that satisfies all independent and dependent variable constraints. With the help of the Fig. 4 the essence of the model predictive control is easily understandable. We formulate the objective function: min Δu(k+ j ) (w(k + j) − y(k + j))2 + λ Δu2 (k + j −1) j=1 Hc ∑ j= H p1 H p 2 ∑ (15) where ∆u(k) denotes the change of the control signal, the Hp1 and Hp2 parameters are the minimum and maximum cost horizons and Hc is the control horizon, which does not necessarily have to coincide with the maximum horizon. λ is a weighting factor, it is a 41 sequence that considers future behaviors, usually constant values or exponential sequences are used. w is the set point signal following the notation of Fig. 4. Figure 4: The essence of model predictive control Generally predictive control uses the receding horizon principle. This means that after the computation of the optimal control sequence, only the first control action will be implemented, subsequently the horizon is shifted one sample and the optimization is restarted with new information about the measurements. In the presence of unmeasured disturbances and modeling errors the MPC controller can exhibit steady- state offset. One way of avoiding this is to design a disturbance estimator which gives the controller implicit integral action. The simplest method for incorporating integral action is to shift the setpoints with the disturbance estimates as depicted in Fig. 5, where the corrected setpoints w′(k) = w(k) – d(k) are modified based on differences between the output of the system and its estimated value d(k) = y(k) – y′(k). Figure 5: The IMC (Internal Model Control) scheme The scheme shown in Fig. 5 is often referred as internal model control (IMC) strategy. This disturbance model assumes that plant/model mismatch is attributable to a step disturbance in the output and that the disturbance remains constant over the prediction horizon. While these assumptions rarely hold in practice, the disturbance model does eliminate offset for asymptotically constant setpoints under most conditions. In practice all processes are subject to constraints. The actuators have a limited field of action as well as determined slew rate, as in the case of valves. Constructive reasons, safety or environmental ones or even sensor slopes themselves, can cause limits in the process variables such as levels in tanks, flows in piping of maximum temperatures and pressures. All of this leads to the introduction of constraints in the MPC problem. Usually, input constraints like umin ≤ u(k + j) ≤ umax, j = 1, ..., Hc (16) ∆umin ≤ ∆u(k + j) ≤ ∆umax, j = 1, ..., Hc (17) are hard constraints in the sense that they must be satisfied. Coversely, output constraints can be often viewed as soft constraints because their violation may be necessary to obtain a feasible optimization problem: ymin ≤ y(k + j) ≤ ymax, j = j1, ..., Hp (18) where j1 represents the lower limit for output constraint enforcement. Because of the sequential solving method, some constraints can be implemented while solving the optimization problem. These constraints are for taking into consideration for example actual physical states of actuators or valves. These constraints can be defined as input constraints, represented by Eq. 16-17. In this study the input constraint were introduced in the following form: u(k+j–1) – ∆u ≤ u(k+j) ≤ u(k+j–1) + ∆u, j = 1, ..., Hc (19) using the fact that the value of ∆u is maximized. Non-linear model based predictive controller Non-linear model-based predictive control (NLMPC) algorithms should be applied in situations where the controlled process is inherently nonlinear, or where large changes in the operating conditions can be anticipated during routine operation, such as in batch processes, or during the start-up and shut-down of continuous processes. The advantages of non-linear predictive control include the following. ● Manipulated and state variable constraints are explicitly handled. ● Nonminimum-phase processes are easily handled. (if the prediction horizon is chosen adequately) ● Knowledge of future setpoint changes is included that is useful for scheduled, coordinated operational changes. The main problem in NLMPC is that a non-linear often (non-convex) optimization problem must be solved at each sampling period in real-time. This hampers the application to fast processes where computationally expensive optimization techniques cannot be used, due to the short sampling time. Several methods can be used to solve such constrained non-linear optimization problems. The most widely studied algorithms are described on [2]. Using the sequential quadratic programming (SQP) method it is possible to minimize the value of the objective function, in each sampling period, varying the parameterized control signal values (u=[u(i)…u(i+Hc)]) on the control horizon. Hence the solution of the optimization problem is a control signal trajectory. The first element of the control signal trajectory is realized 42 in the next time sequent and the other elements are neglected. This is called recending horizon strategy. Figure 6: The scheme of the non-linear model predictive controller Model predictive control of a district heating network Energy producer are compelled to reduce their rate of polluting emissions beside fulfilling consumers power demands with the lowest global costs. Thus, technical, economical and environmental constraints have to be simultaneously dealt with. The optimization problem stated from this multi- field area can hardly be solved as it is a non-linear programming problem, consists of numerous variables. The optimal control of district heating networks, for which propagation delays can not be neglected and mechanical and thermal losses have non-linear expressions, picks up all these harsh difficulties. Managing a district heating network implies to assign values to integer variables (status of production units, status of pumps…) and to continuous variables (amounts of energy to produce). As a result, the optimization of the production and energy supply planning appears to be a huge, mixed and non linear optimization issue. However solving a non-linear mixed-integer optimization problem might have the ability to provide a control signal that may provide better performance than solving a non-linear optimization problem. In this work a non-linear SQP method with soft constraints will be introduced to avoid the complexity of mixed-integer non-linear programming. To take the different weights of the control variables into consideration the objective function is augmented with the absolute value of the control variables: ∑ ∑ ∑ = = =+ ++ +−+++−+ 2 2 1 1 2 1 22 )( )( )1())()((min p p p c H j H Hj H jjku jku jkujkyjkw α λβ Δ Δ (20) It can also be important to define weights (β) for the error, the set points and manipulated variables, because the main task – keeping the manipulated variable equal to the set point – can be easily assured. This version of the objective function will be applied in this study. Differences between the models with and without time delay In the case of the depicted (Fig. 1) district heating network the possible control variables are: ● invested heat in Production unit 1 ● invested heat in Production unit 2 ● pump duty of P1 ● pump duty of P2 ● valve opening. Since the pump P1 is chosen to compensate the pressure drop of the heat exchangers and pipelines, the P1 pump does not take part in satisfying the heat demand of consumers, so it is controlled by a local regulator. The split ratio between Consumer 2 and Consumer 3 can be adjusted using two control variables: valve opening and the pump duty of the P2 pump (as the pressure increase on the pump is the function of the pump duty (Eq 3)). These control variables determine the flow in the two directions and thus the transferred heat to the consumers. Comparing the model with and without time delay a main difference can be seen in the dynamics of the models. It is caused by the absent time delay. Fig. 7 shows the difference between the two models (full line means the model without time delay, dot line means the model with time delay). Naturally both simulations were run with the same inputs. In the initial moment there is no temperature profile in the network so the first 15 minutes show the temperature startup of the district heating network. After 15 minutes there is a step in the heat from Producer 1. After the network is in steady state, there is a new step in the heat from the Producer 2. In the 50th minute there is a step in the pump duty of Pump 2 The last change in the inputs is changing the value of the valve opening in 63rd minute. These changes can be seen on the Fig. 8. Figure 7: Transferred heat in model ‘A’ and ‘B’ model to the same input 43 Figure 8: The input parameters of the district heating network related to Fig. 7 From the Fig. 7 it can be stated that the model without time delay can provide almost the same steady state than the model with time delay. The main difference is in the dynamical behavior. Simulation scenarios The main goal is to satisfy the heat demand of the consumers. In this study the heat demand of the consumers is assumed to be known. It us to be tracked as the set points of the district heating network. To test the performance of the DHN the following set point trajectories were chosen: The set point changes are at the same time: at around 33rd minute and 60th minute. In order: ● Consumer 1 (lower line) has 25000 kW – 60000 kW – 30000 kW heat demand ● Consumer 2 (middle line) has 21000 kW – 55000 kW – 27000 kW heat demand ● Consumer 3 (upper line) has 28000 kW – 55000 kW – 32000 kW. The main goal is to minimize the transition time as possible and at the same time fulfill the heat requirements of the consumers (follow the previously presented set point trajectory) considering to minimize of the use of the Production unit 2 and the pump duty of P2 pump. This goal can be reached with minimizing the objective function presented in Eq. 20. Optimization Model A contains time delays in the pipes: so the model and the process are identical. In this case it is not necessary to use the IMC scheme, because there is no mismatch, and because there are no unmeasured disturbances. Model B has no, so the plant-model mismatch seen in Fig. 7 is expected to deteriorate control performance. Because of this model mismatch simulation experiments will be carried out to test the performance of the MPC with IMC scheme and without IMC scheme. Since the steady states of the two models should be equal the MPC was expected to reach the set point trajectories in both cases, but the dynamical behavior can not be predictable without simulations. Simulation results Control performance of MPC with A and B model The ISE (Integral Square of Error) criteria was chosen to compare the control performance of the two controllers (A and B model).Additionally plots are used to compare quality of control. The set points are the same in all cases. The tuning parameters of model predicitve controllers are the lenght of the prediction and control horizon and the value of α, β and λ parameters (Eq 20). To be able to compare the scenarios in both cases the same tuning parameters were chosen (prediction horizon: 4, control horizon: 1, sample time: 45 sec). The computation demand of the NLMPC controller is very high, since the SQP algorithm obtains gradient information via finite differences and the differential-algebraic model equations have to be solved for each perturbation. This solving method is very time-consuming and this is the reason why the tuning parameters of the NLMPC controller have not been optimized. Control results using the ‘A’ model ‘A’ model means that the applied model during the optimization is the same as the model used as the operation process. Thus there is no model mismatch is expected so there is no need to apply the IMC scheme. The constraints are formulated like Eq 19. The applied constraint can have serious effects on the computed control signal as seen the Fig. 9 and Fig. 10. Analyzing the first transition: to assure a lower flow rate in direction of Consumer 3 it was necessary to set the higher pressure drop of that section. To reach this goal infinite variations of valve opening-pump duty value pairs exist. By applying proper weight in objective function the biggest valve opening – the lowest pump duty pair can be applied. In the steady state of the system this condition is obviously determined. In the unconstrained case pump duty is decreased quickly to increase the pressure drop in the direction of flow, but at the same time the valve opening is set as high as possible – 100%. In the constrained case the controller handles the changes of these manipulated variables differently since it is not permitted to change them arbitrary. The increasing of that direction is handled by closing the valve and at the same time reduce the pump duty. Closing the valve is a necessary action since the change of the pump duty is constrained. 44 Figure 9: Comparing the control performance of the unconstrained (dashdot line) and constrained (dot line) MPC Figure 10: The manipulated variables regarding to the previous figure full line – constrained, dashed line – unconstrained MPC The quick and accurate control actions are not surprising since the most precise model was applied during the optimization, but at the same time this accuracy has enormous computational demand: simulation needs almost 8–10 hours to be finished. In the further cases only constrained MPCs will be described using the assumption that in unconstrained cases they would provide similar control action as it was recently shown. Control results using the ‘B’ model Using the ‘B’ model implies that the applied model in the optimization is not the same as the model used in as an operation process. For this scenario 3 different cases will be introduced: ● Case 1: MPC without IMC scheme ● Case 2: MPC with IMC scheme ● Case 3: MPC, combination of the previous cases. Case 1 – MPC without IMC scheme Here the control algorithm does not have any information about the model mismatch. The control performance mainly depends on the difference of the time delay. In the previous case this did not cause a problem since the model was identical, but in this case this requirement is not fulfilled. The control performance can be illustrated by Fig. 11. As the Fig. 11 shows, the existing model mismatch causes steady state offset, mainly detectable after the first transient. However this MPC has the advantage of avoiding any overshoot. The steady state offset – if it is converted to temperature difference – means 2–3 °C difference. Since the same objective function was used in the optimization section in all cases – the same penalty weights were applied for the change of inputs – the control variable trajectories kept all the characteristics that were introduced in the case of using constrained MPC using ‘A’ model. The computational time was reduced to the half of the previous simulation, the simulation needs almost 3–4 hours to be accomplished. Figure 11: The control performance of the model predictive controller in Case 1, ‘B’ model without IMC scheme 45 Figure 12: The computed control signal corresponding to Fig. 11 Figure 13: The control performance of the model predictive controller in Case 2 Figure 14: The computed control signal regarded to the control scenario depicted on Fig. 13 Figure 15: The modified set point signal (by IMC scheme) (full line) and the output of the network (dashed line) Case 2 – MPC with IMC scheme In order to eliminate steady state offset IMC scheme is applied. After running a simulation using IMC scheme, the following control result can be observed: In contrast to Case 1 and simulation using ‘A’ model in this case some overshoot and oscillation occurres. A reason for this could be that the set point signal is modified with the error of the model and the plant Fig. 15 shows the modified set point signal, which is the input of the optimization section and used in computing the minimum of the objective function. The oscillations are mainly caused by the absent time delay in the model and it is the most obvious during transitions. After the oscillations have died away at steady state there is no offset. Case 3 – combination on Case 1 and Case 2 In this case an attempt will be made to combine the advantages of Case 1 and Case 2, avoiding overshoot and eliminating the steady state offset. To reach this goal the following strategy is applied: since the IMC structure modifies the set point signals significantly during transitions, it is not advantageous to apply this scheme during the transitions. At the same time it is very useful to apply the IMC scheme to eliminate the steady state offset. So in this case a trigger is implemented in the controller to switch on the IMC scheme. The trigger is formulated with the following expression: 46 K N iMEiME ≤ −− 2))1()(( (21) where: ME – the model error vector in ith and (i-1)th sample time N – length of the model error vector K – constant. So if the change of the model error is smaller than a previously determined constant, it means that the manipulated variable is relatively close to the steady state. If this condition is fulfilled the IMC scheme is switched on and eliminate the steady state offset. By applying this method the following control result can be yielded: As the figures shows, this method can extract the advantage of IMC – no steady state offset - and at the same time exclude the disadvantage of IMC – oscillations in the set point signal. Figure 16: The control performance of the model predictive controller in Case 3 Figure 17: The computed control signal regarded to the control scenario depicted on Fig. 16 Performance evaluation of the different MPC systems In this section we compare the different MPC using a the Integral of Square Error (ISE), and graphical plots. Additionally the control performance is examined in the point of view of settling time and the existence of overshoot during the control scenarios. Furthermore we will explain the occurring differences. General comparison In this short section the controllers will be compared in some general aspects: ● Performance index ● Settling time and overshoot. Performance index The integral square error (ISE) is a measure of the control performance. It is obtained by integrating the square difference of set point signal and controlled variable over the time interval of the simulation. As the values of the controlled variables only known at sample times, the ISE can be approximated with the following expression: ∑ = −= N i ii ywISE 1 2)( (22) Where wi means the value of the set point in ith moment, yi is the output of the system, N is the number of time steps. The constrained ‘A’ model provides the best performance in terms of ISE. This is not surprising, since the controller uses the most accurate model to predict the reaction of the plant to a certain input. In cases that use ‘B’ model the controllers have worse performance because of plant – model mismatch. Case 1 and Case 3 has equal ISE value, despite the combination of IMC/noIMC scheme. It is found that introducing switching strategy benefits only for Consumer 2. Settling time and overshoot Next we consider the settling time and overshoot. In tables 2 and 3 the performance of the control scenarios is summarized. Table 1 contains the ISE value of the previously presented simulations. 47 Table 1: Comparing the performance of the applied MPCs by ISE value ISE (*109) Consumer 1 Consumer 2 Consumer 3 Mean % Case 'A' model 1.7 1.45 1.6 1.58 100% Case 1 2.54 2.31 2.72 2.52 159% Case 2 2.52 2.12 3.14 2.59 164% Case 3 2.56 2.27 2.74 2.52 159% Table 2: Comparing the applied MPCs by the existing of overshoot Overshoot Startup Transient 1 Transient 2 Case 'A' model no yes no Case 1 no no no Case 2 yes yes yes Case 3 no no no Table 3: Comparing the applied MPCs by settling time Settling time (min) Startup Transient 1 Transient 2 Case 'A' model 7 7 9 Case 1 11 12 13 Case 2 25 28 24 Case 3 15 15 13 Analysis of the performance of the ‘B’ model The main advantage of applying this model – reduced computational time. Next we highlight some important characteristics. Because of the plant – model mismatch and the lack of IMC scheme the MPC in Case 1 can not reach the set point signal as it can be seen on Fig. 12. To eliminate this phenomenon the IMC scheme was applied. Fig. 18 shows the result (for Transient 1). As mentioned before, the overshoot of IMC-MPC has been caused by the significant model error in the transients. Taking model error into consideration is not a negligible fact since the performance of the model predictive control is a function of the model parameters, and model parameters can be the function of the time (eg. fouling in the heat exchangers can change the heat transfer coefficient). Comparing the Case A model and Case 3 Fig. 19 shows the difference of these two controllers. Figure 18: Graphical comparison of control performance of Case 1 and Case 2 controller Figure 19: Graphical comparison of control performance of Case ‘A’ model and Case 2 controller 48 Figure 20: Control variables of Case ‘A’ model and Case 3 The main difference is the settling time. The Case 3 controller reaches the set point slower than the controller which uses the ‘A’ model. This happens because ‘A’ model and ‘B’ model have different dynamics (‘B’ model has no time delay). Thus the ‘B’ model can reach the set point signals faster by the effect of the same input signal than ‘A’ model. That is why the ‘B’ model considers the transient finished sooner than it is realized by the operating process (‘A’ model). In steady state both models responsed with almost the same output (for the same input). To reach the set point signal without steady state offset the IMC scheme switches on and eliminates the offset. Summary In the previous sections detailed and general comparisons of the applied controllers was given. The Case ‘A’ model was able to provide the highest accuracy that was confirmed by the lowest ISE value and settling time. The most significant disadvantage of Case ‘A’ model was the enormous computational demand that was reduced by using the ‘B’ model in the optimization. In Case 1 it was observed, that a less accurate model was also useful for control purposes, however the performance was not as good as in Case A, since there was a steady state offset and the settling time was higher than in Case A. The Case 1 controller had the drawback of the lack of the feedback of the difference of the model and plant. This feedback was realized in Case 2 by introducing IMC scheme. It had the advantage of eliminating the steady state offset of Case 1, but unfortunately it accentuated the model mismatch in the dynamics. In Case 3 the beneficial characteristics – eliminating the steady state offset caused by the model mismatch and avoiding the overshoot – of Case 1 and Case 2 were combined. This kind of solution could provide the most attractive performance. It can be difficult to predict the behavior the controller and the controlled system, compared to the Case ‘A’ model. Finally it can be stated that it is beneficial to use a model which is as accurate as possible but at the same time it is important to consider the fact that the computation demand is increasing with the increasing model accuracy. In case of plant – model mismatch the use of IMC scheme can be useful since it can handle the mismatch and the effect of unmeasured disturbances at the same time. Outlook and future work While designing the non-linear MPC framework it was very important to keep it modular, so changes are easy to implement. In this field the two most important characteristics of the control algorithms are: ● to have the opportunity to implement the control algorithm in the control system in a short time ● the ability of the control algorithm to provide a feasible control signal in a shorter time interval than the sampling time. For applying a non-linear model predictive algorithm it is very important to require that the optimization algorithm can find the (global) minimum of the objective function or at least a feasible solution in an certain time interval. Considering that numerical optimization in each sequence can be very time consuming it is necessary to implement the methods which can reduce the computational demand of the optimization process. Some possible approaches for this are: ● Model reduction, which yield less states in the model ● Applying the gradient of the objective function during the optimization. When it can be realized, these kind of control algorithms can be used widespread in the industrial practice as a real time optimization algorithm, regarding that the MPC algorithms can be applied in the advanced control level. ACKOWLEDGEMENT János Abonyi is grateful for the support of the Bolyai Research Fellowship of the Hungarian Academy of Sciences. The financial support of the TÁMOP-4.2.2- 08/1/2008-0018 project is gratefully acknowledged. REFERENCES 1. BENONYSSON A., BØHM B., RAVN H. F.: (1995) Operational optimization in a district heating system. Energy Conversion and Management Vol 36 (5) (1995), 297-314. 2. HENSON M. A.: (1998) Nonlinear model predictive control: Current status and future directions. Computers and Chemical Engineering, 23, 187-202 49 3. MARCHETTY J. L., MELLCHAMP D. A., SEBORG D. E.: (1983) Predictive Control Based on Discrete Convolution Models, Ind. Eng. Chem. Res. Dev., 22 488-495 4. MATHISEN K. W., MORARI M., SKOGESTAD S.: (1993) Dynamic models for heat exchangers and heat exchanger networks, European Symposium on Computer Aided Process Engineering 3 Graz, Austria 5. NIELSEN A. H., MADSEN H.: (2006) Modelling the heat consumption in district heating systems using a grey-box approach Energy and Buildings 38 (2006) 63-71 6. RAVN H. F., RYGAARD J. M.: (1994) Optimal scheduling of coproduction with a storage. Engineering optimization Vol 22 (1994), 267-281. 7. SANDOU G., FONT S., TEBBANI S.: Global modelling and simulation of a district heating network Supélec, Service Automatique 8. SANDOU G., FONT S., TEBBANI S., HIRET A., MONDON C.: (2005) Predictive Control of a Complex District Heating Network, Proceedings of the 44th IEEE Conference on Decision and Control, and the European Control Conference 2005, Seville, Spain, December 12-15, 2005 9. SHRIDHAR R., COOPER D. J.: (1997) A Tuning Strategy for Unconstrained SISO Model Predictive Control, Ing. Eng. Chem. Res., 36, 729-746 10. SØGAARD H. T.: (1993) Stochastic systems with embedded parameter variations – Applications to district heating. Ph.D. Dissertation. Technical University of Denmark, Institute of Mathematical Statistics and Operations Research. 11. SUBIR S., KOTHARI D. P.: (1998) Optimal thermal generating unit commitment: a review. Electrical Power & Energy Systems Vol 20 (7) (1998), 443-451. 12. WEBER C., MARECHAL F., FAVRAT D.: Network synthesis for district heating with multiple heat plants, The International Conference on ENERGY and ENVIRONMENT, CIEM 2005 13. ZHAO H., HOLST J., ARVASTSON L.: (1998) Optimal operation of coproduction with storage. Energy Vol 23 (10) (1998), 859-866.