# Neural Network Based Central Heating System Load Prediction and Constrained Control

Authors: Hongwei Wang; Fangwen Tu; Baofeng Tu; Guohui Feng; Guangming Yuan; Hao Ren; Jiarong Dong
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

A neural network (NN) based heating system load prediction and control scheme are proposed. Different from traditional physical principle based load calculation method, a multilayer NN is incorporated with selected input features and trained to predict the heating load as well as the desired supply water temperature in heating supply loop. In this manner, a complicated load calculation model can be replaced by simple but efficient data-driven scheme and the response time to outdoor temperature variation can be enhanced. Moreover, in order to handle the input and output constraints in valve opening degree control task to achieve desired supply water temperature, Barrier Lyapunov candidate function and axillary system technique are involved. An additional NN is employed to approximate the system transfer function with reliable accuracy. The stability of the system is guaranteed through rigorous mathematical analysis. The excellent performance of the novelly proposed control over traditional PID is demonstrated via extensive simulation study. A quantitative case study is also conducted to verify the flexibility and validity of proposed load prediction strategy.


---

## Body

## 1. Introduction

Nowadays, heating service is essential and significant to the wellness of human’s daily life as well as the economy development. Attributing to the improvement of technology, central heating system has been widely studied and applied to modern municipal service industry. For example, in Finland, approximately 80% heating service is provided by central heating system, which contributes a lot to the energy conservation and environment protection. In this system, sensors measure the outdoor temperature and other information to guide the adjustment of opening valve degree in primary pipe network which in turn controls the water temperature in secondary pipe network to serve the users.In practical operation, it is usually very complicated to build physical model for heating load calculation which leads to a lagged response for the outdoor temperature variation. In this regard, a data-driven method for the load prediction is able to facilitate the adjustment of heating control scheme. In literature, there are some researches that take advantage of the historical data to achieve nonlinear heating consumption of load prediction [1]. An artificial neural network based heating system mathematical model is investigated to estimate the heating energy consumption for Canadian residence in [2]. The input features’ saliency to the heating load is researched in [3]. The result shows that outdoor temperature is the key factor that affects the load of a heating system. This research paves the way for the input feature selection of training models. In [4], Support Vector Machine (SVM) is incorporated into the building energy consumption prediction using the instances of four commercial buildings in Singapore. The input features are selected as average month outdoor temperature, relative humidity, sunshine duration, and the prediction error being guaranteed below 4%. Apart from these works, some novel methods are also reported for the prediction such as [5] which applies Finite Impulse Response (FIR) method for critical point temperature prediction and [6] which adopts data-mining methods to predict heating steam. However, as far as the authors’ knowledge, there is still no research that has been done on the short-term heating load prediction for performance upgrade of heating system to achieve fast response against temperature variation. Thus, in this paper, we propose an NN based predictor to supervise the control system operation with reliable accuracy.Automatic control has been widely applied in central heating system. An optimal control strategy is proposed in [7] by building indirectly connected district heating system model. Reference [8] introduces a reginal automatic control scheme by setting temperature monitor and valve at the user side to improve the heating quality. As for the proposed heating architecture, the aforementioned load prediction module is regarded as a supervisor to control system by advising the control objective which concretely is the desired supply water temperature in the secondary pipe network. With this desired temperature trajectory, an intelligent control method is proposed with another NN to approximate the unknown nonlinear transfer model between supply water temperature and the valve opening degree in primary pipe network. But in implementation, in order to achieve accurate temperature control and at the same time consider the opening degree limitation of the valve, the output tracking error as well as the input control effort should be constrained. It brings additional challenge for the control design. Some researches have been done for the input constraint problem [9–11]. In [12, 13], an auxiliary system is novelly introduced to handle the constraint effect in a class of uncertain multiple-input multiple-output nonlinear system. In [14], a constrained-input system is tackled in combination with the optimal control to ensure a good tradeoff between control performance and energy consumption. For output constraints, artificial potential field [15], prescribed performance control [16, 17], model predictive control [18], and reference governor [19] are some of the existing strategies to handle this problem. In [20, 21], Barrier Lyapunov Function (BLF) is introduced which needs less initial conditions and does not require explicit system solution. Based on these literature reports, in this paper, the input and output constraint problem is tentatively merged with rigorous mathematical stability analysis to ensure accurate heating temperature control and limited valve opening degree. Block diagram of the overall proposed central heating system is presented in Figure 1.Figure 1
Diagram of overall system architecture.The contributions of this paper are threefold:(i)
NN based heating load prediction module is developed containing components such as data preprocessing, input variable selection, and NN architecture establishment. It plays a role of supervisor in the whole system to ensure fast response to the outdoor temperature variation.(ii)
Intelligent adaptive control is proposed in combination with additionally involved auxiliary subsystem and symmetric BLF to cope with coupled input saturation and output tracking error constraints simultaneously for accurate and efficient control.(iii)
A function approximator is incorporated into the control loop in order to handle the unknown nonlinear system parameters on NN basis. Validation of the proposed heating system is demonstrated through case study.The organization of this paper is as follows. Section2 describes the specifics of NN based heating load prediction. Section 3 presents the development of adaptive control system with multiple constraints. Numerical simulation on the proposed central heating framework is conducted to demonstrate the effectiveness in Section 4. Section 5 gives some concluding remarks. The main denotations that will be used in this paper are summarized in the Nomenclature section.


## 2. Heating Load Prediction via Neural Networks

### 2.1. System Architecture and Mathematical Model

In this paper, an indirectly connected central heating system is considered due to its high efficiency and energy conservation properties. The overall architecture of this system can be described in Figure2. Heat exchanger is the bridge that connects primary and secondary network. These two loops transfer heat through heat exchanger. This process is subject to heat balance principle. Since the secondary network is directly connected to the users and influences the heating quality, in this work, we focus on the analysis of this loop. To qualitatively analyze the heat supply system, the following mathematical model is considered [22]:(1)Q1=Q2=Q3,(2)Q1=qvVwTn-Tw,(3)Q2=a·FwTg2+Th22-Tn1+b,(4)Q3=G2cwTg2-Th2,where Q1 denotes design space-heating load of the building, Q2 is heat from the radiator, and Q3 represents the heat supplied to the users by heat supply loop. qv, Vw, a, Fw, G2, cw denote the heating index, building volume, scalar parameter, heat radiating area of the radiator, design flow of secondary network, and water specific heat capacity, respectively. b equals 0 for hot water floor radiation heating. Tn and Tw are indoor and outdoor temperature. Tg2 and Th2 denote supply and return water temperature in secondary network. According to the heat balance principle, the three heat loads Q1, Q2, Q3 are supposed to be equalled as shown in (1).Figure 2
System architecture of indirectly connected district heating system.Assumption 1.
Secondary network is usually subject to 3–5% heat loss. For simplicity purpose, the loss is ignored in this work.In this section, the daily real-time heating load and desired secondary network supply water temperatureTg2 are going to be predicted through artificial neural networks with historical information.


### 2.2. NN Based Heating Load andTg2 Prediction

To achieve the prediction, a thorough and precise data base is required to be built. Thus, various sensors are employed to collect the real-time operating system states and other related parameters such as hourly based supply water temperature, indoor/outdoor temperature, and other meteorological parameters. These parameters can act as training samples for the prediction task.


#### 2.2.1. Selection of Input and Output Variables

In practical engineering, the selection of input variables is significant for the performance of predictor. During this task, the variables’ relationship with heating load is investigated and the ones with strong impact to the target load are chosen as input variables. There are plenty of factors that affect the load. The most straightforward one is outdoor meteorological condition. Moreover, the inherent properties of the building such as types, heat characteristics of building envelope, and utility conditions will also influence the indoor temperatures. Nevertheless, it is usually difficult to collect historical data on all these aspects. Furthermore, too many input variables will contribute to a complicated prediction model and retard the training period. The nonindependent input variables can also degrade the prediction accuracy. Therefore, a comprehensive consideration should be conducted before decision of input variables.Through analysis, the input variables can be categorized into (i) heating system: supply water temperature, return water temperature, water supply flow; (ii) outdoor meteorological parameters: outdoor temperature, outdoor wind speed, air pressure, weather condition, relative humidity, and solar irradiation intensity; (iii) time information: workday or not, daytime, or night. Among them, water supply flow, supply temperature, and return temperature can be collected through heat exchange station record. Meteorological parameters can be collected from meteorological station. By investigating the data source and conducting trial and error experiments, a four-dimension input variables are determined, that is, time (tp), outdoor temperature (Tw), solar irradiation intensity (Rp), and the supply water temperature in previous day (Tg2p).In this paper, a temperature control based heating system is considered. Concretely, the supply water temperature in secondary networkTg2 is decided as the output. The heating load can be further calculated with predicted Tg2 using (1)–(4).


#### 2.2.2. Data Preprocessing

Due to measurement error, human caused error, data transmission interference, or equipment precision problem, bad data exists in historical data sequence. These bad data cannot reflect the variation heating load which needs to be detected before they are imported to the network. To achieve this, data of the same category forms a matrixXij, where i denotes the recording time in one day and j represents the recording date. The average value Ei and variance Di of each row are calculated first. Then the skew ratio of each data Pij is obtained. Define a threshold σ=1.1. When Pij>σ, the corresponding data is judged as bad data and removed from the data set. A linear interpolation method is applied to generate new reliable data to fill the vacancy.The data is subject to normalization after bad data detection to improve the training speed and prediction accuracy. Its objective is to fit the input and output data in the range of[-1,1] with the following operation.(5)y=x-0.5xmax+xmin0.5xmax-xmin,(6)x=0.5·xmax-xmin·y+0.5·xmax+xmin,where x is the original data before normalization and y denotes the data after normalization. xmax and xmin are the maximum and minimum data in the original data set. Equation (5) is utilized to transform original data into normalized data at input layer. The result from output layer will pass through (6) to obtain the normal-scale value.


#### 2.2.3. Neural Network for Training

In this paper, a three-layer neural network is proposed with input, hidden, and output layer. The number of nodes in input layer is set as 4 due to the dimension of input variables. The output layer node number is 1. Regarding the hidden layer node number, Kolmogorov Theorem [23] is employed as n=2k+1, where n and k denote node number of hidden and input layers. It gives 9 for the node number. Backpropagation [24] algorithm with batch training mode is applied for the training of the proposed network.


## 2.1. System Architecture and Mathematical Model

In this paper, an indirectly connected central heating system is considered due to its high efficiency and energy conservation properties. The overall architecture of this system can be described in Figure2. Heat exchanger is the bridge that connects primary and secondary network. These two loops transfer heat through heat exchanger. This process is subject to heat balance principle. Since the secondary network is directly connected to the users and influences the heating quality, in this work, we focus on the analysis of this loop. To qualitatively analyze the heat supply system, the following mathematical model is considered [22]:(1)Q1=Q2=Q3,(2)Q1=qvVwTn-Tw,(3)Q2=a·FwTg2+Th22-Tn1+b,(4)Q3=G2cwTg2-Th2,where Q1 denotes design space-heating load of the building, Q2 is heat from the radiator, and Q3 represents the heat supplied to the users by heat supply loop. qv, Vw, a, Fw, G2, cw denote the heating index, building volume, scalar parameter, heat radiating area of the radiator, design flow of secondary network, and water specific heat capacity, respectively. b equals 0 for hot water floor radiation heating. Tn and Tw are indoor and outdoor temperature. Tg2 and Th2 denote supply and return water temperature in secondary network. According to the heat balance principle, the three heat loads Q1, Q2, Q3 are supposed to be equalled as shown in (1).Figure 2
System architecture of indirectly connected district heating system.Assumption 1.
Secondary network is usually subject to 3–5% heat loss. For simplicity purpose, the loss is ignored in this work.In this section, the daily real-time heating load and desired secondary network supply water temperatureTg2 are going to be predicted through artificial neural networks with historical information.


## 2.2. NN Based Heating Load andTg2 Prediction

To achieve the prediction, a thorough and precise data base is required to be built. Thus, various sensors are employed to collect the real-time operating system states and other related parameters such as hourly based supply water temperature, indoor/outdoor temperature, and other meteorological parameters. These parameters can act as training samples for the prediction task.


### 2.2.1. Selection of Input and Output Variables

In practical engineering, the selection of input variables is significant for the performance of predictor. During this task, the variables’ relationship with heating load is investigated and the ones with strong impact to the target load are chosen as input variables. There are plenty of factors that affect the load. The most straightforward one is outdoor meteorological condition. Moreover, the inherent properties of the building such as types, heat characteristics of building envelope, and utility conditions will also influence the indoor temperatures. Nevertheless, it is usually difficult to collect historical data on all these aspects. Furthermore, too many input variables will contribute to a complicated prediction model and retard the training period. The nonindependent input variables can also degrade the prediction accuracy. Therefore, a comprehensive consideration should be conducted before decision of input variables.Through analysis, the input variables can be categorized into (i) heating system: supply water temperature, return water temperature, water supply flow; (ii) outdoor meteorological parameters: outdoor temperature, outdoor wind speed, air pressure, weather condition, relative humidity, and solar irradiation intensity; (iii) time information: workday or not, daytime, or night. Among them, water supply flow, supply temperature, and return temperature can be collected through heat exchange station record. Meteorological parameters can be collected from meteorological station. By investigating the data source and conducting trial and error experiments, a four-dimension input variables are determined, that is, time (tp), outdoor temperature (Tw), solar irradiation intensity (Rp), and the supply water temperature in previous day (Tg2p).In this paper, a temperature control based heating system is considered. Concretely, the supply water temperature in secondary networkTg2 is decided as the output. The heating load can be further calculated with predicted Tg2 using (1)–(4).


### 2.2.2. Data Preprocessing

Due to measurement error, human caused error, data transmission interference, or equipment precision problem, bad data exists in historical data sequence. These bad data cannot reflect the variation heating load which needs to be detected before they are imported to the network. To achieve this, data of the same category forms a matrixXij, where i denotes the recording time in one day and j represents the recording date. The average value Ei and variance Di of each row are calculated first. Then the skew ratio of each data Pij is obtained. Define a threshold σ=1.1. When Pij>σ, the corresponding data is judged as bad data and removed from the data set. A linear interpolation method is applied to generate new reliable data to fill the vacancy.The data is subject to normalization after bad data detection to improve the training speed and prediction accuracy. Its objective is to fit the input and output data in the range of[-1,1] with the following operation.(5)y=x-0.5xmax+xmin0.5xmax-xmin,(6)x=0.5·xmax-xmin·y+0.5·xmax+xmin,where x is the original data before normalization and y denotes the data after normalization. xmax and xmin are the maximum and minimum data in the original data set. Equation (5) is utilized to transform original data into normalized data at input layer. The result from output layer will pass through (6) to obtain the normal-scale value.


### 2.2.3. Neural Network for Training

In this paper, a three-layer neural network is proposed with input, hidden, and output layer. The number of nodes in input layer is set as 4 due to the dimension of input variables. The output layer node number is 1. Regarding the hidden layer node number, Kolmogorov Theorem [23] is employed as n=2k+1, where n and k denote node number of hidden and input layers. It gives 9 for the node number. Backpropagation [24] algorithm with batch training mode is applied for the training of the proposed network.


## 2.2.1. Selection of Input and Output Variables

In practical engineering, the selection of input variables is significant for the performance of predictor. During this task, the variables’ relationship with heating load is investigated and the ones with strong impact to the target load are chosen as input variables. There are plenty of factors that affect the load. The most straightforward one is outdoor meteorological condition. Moreover, the inherent properties of the building such as types, heat characteristics of building envelope, and utility conditions will also influence the indoor temperatures. Nevertheless, it is usually difficult to collect historical data on all these aspects. Furthermore, too many input variables will contribute to a complicated prediction model and retard the training period. The nonindependent input variables can also degrade the prediction accuracy. Therefore, a comprehensive consideration should be conducted before decision of input variables.Through analysis, the input variables can be categorized into (i) heating system: supply water temperature, return water temperature, water supply flow; (ii) outdoor meteorological parameters: outdoor temperature, outdoor wind speed, air pressure, weather condition, relative humidity, and solar irradiation intensity; (iii) time information: workday or not, daytime, or night. Among them, water supply flow, supply temperature, and return temperature can be collected through heat exchange station record. Meteorological parameters can be collected from meteorological station. By investigating the data source and conducting trial and error experiments, a four-dimension input variables are determined, that is, time (tp), outdoor temperature (Tw), solar irradiation intensity (Rp), and the supply water temperature in previous day (Tg2p).In this paper, a temperature control based heating system is considered. Concretely, the supply water temperature in secondary networkTg2 is decided as the output. The heating load can be further calculated with predicted Tg2 using (1)–(4).


## 2.2.2. Data Preprocessing

Due to measurement error, human caused error, data transmission interference, or equipment precision problem, bad data exists in historical data sequence. These bad data cannot reflect the variation heating load which needs to be detected before they are imported to the network. To achieve this, data of the same category forms a matrixXij, where i denotes the recording time in one day and j represents the recording date. The average value Ei and variance Di of each row are calculated first. Then the skew ratio of each data Pij is obtained. Define a threshold σ=1.1. When Pij>σ, the corresponding data is judged as bad data and removed from the data set. A linear interpolation method is applied to generate new reliable data to fill the vacancy.The data is subject to normalization after bad data detection to improve the training speed and prediction accuracy. Its objective is to fit the input and output data in the range of[-1,1] with the following operation.(5)y=x-0.5xmax+xmin0.5xmax-xmin,(6)x=0.5·xmax-xmin·y+0.5·xmax+xmin,where x is the original data before normalization and y denotes the data after normalization. xmax and xmin are the maximum and minimum data in the original data set. Equation (5) is utilized to transform original data into normalized data at input layer. The result from output layer will pass through (6) to obtain the normal-scale value.


## 2.2.3. Neural Network for Training

In this paper, a three-layer neural network is proposed with input, hidden, and output layer. The number of nodes in input layer is set as 4 due to the dimension of input variables. The output layer node number is 1. Regarding the hidden layer node number, Kolmogorov Theorem [23] is employed as n=2k+1, where n and k denote node number of hidden and input layers. It gives 9 for the node number. Backpropagation [24] algorithm with batch training mode is applied for the training of the proposed network.


## 3. NN Based Adaptive Control with Input and Output Constraints

### 3.1. Preliminaries on NN Approximation

Radial Basis Function Neural Network (RBFNN) has been proven to have excellent function approximation and learning capabilities. Given a continuous nonlinear functionϕ(z):Rn→Rm on a compact set Ωz, it can be approximated by RBFNN in the following expression:(7)ϕNNz=WTSz,where z∈Ωz⊂Rn is the input vector. W∈Rp×m denotes the adjustable weight matrix; p stands for the neuron number. S(z)=[s1(z),…,sp(z)]T represents the basis function vector with si(z)=exp[-(z-μi)T(z-μi)/ψi2], i=1,2,…,p. μi=[μi1,μi2,…,μin] is the center of respective field; ψi is the width of Gaussian function. More accurately, (7) can be rewritten as(8)ϕz=W∗TSz+εz,where W∗ denotes the ideal weight; ε(z)∈Rm represents the approximation error and satisfies ε(z)≤σ, where σ is a positive constant. The ideal weight matrix W∗ is defined as(9)W∗≔argminW∈Rp×m⁡supz∈Ωzϕz-WSz.W∗ is merely an “artificial" quantity for analysis purpose. In control development, it needs to be estimated by W with designed updating law [25]. Additionally, NN approximation error ε(z) indicates the minimum deviation between optimal approximation solution W∗S(z) and the unknown function ϕ(z). It has been demonstrated that NN approximation can achieve any desired accuracy if the number of neurons is large enough; in other words, ε(z) can be arbitrarily small if p is sufficiently large [25].


### 3.2. Adaptive Control with Multiple Constraints

With the predicted supply water temperatureTg2 in secondary network, we intend to control the valve opening degree in primary network to ensure the supply water temperature tracking the prediction value. Herein, the prediction value Tg2 in the previous section acts as the desired signal denoted as Tg2d. The transfer function in quality adjustment is presented as follows [26].(10)G1s=K1T1s+1T2s+1e-τ1s, where K1 is a scalar coefficient, T1 and T2 are inertial time constant, and τ1 is the lagging time. This model describes the relationship between system input Vg1, that is, opening degree of the valve in primary pipe network, and the system output Tg2, that is, actual water temperature in secondary pipe network. To proceed, we transfer the above expression into time-domain formulation.(11)T1T2K1T¨g2+T1+T2K1T˙g2+Tg2K1=Vg1t-τ1.Equation (11) can be concisely expressed with(12)MT¨g2+CT˙g2+gTg2=Vg1t-τ1,where M=T1T2/K1, C=T1+T2/K1, g(Tg2)=Tg2/K1. Without loss of generality, we make the following remarks and assumptions.Remark 2.
Although the parametersK1,T1, and T2 in [26] are assumed as constants, it is obvious and straightforward that the transfer process is time-varying; that is, the aforementioned parameters are state-dependent and time-varying. Subsequently, M,C, and g(Tg2) should be more reasonably depicted as M(T˙g2,Tg2,t),C(T˙g2,Tg2,t), and g(Tg2,t).Remark 3.
A bounded disturbance term should be complemented in model (12) to consider the system noise and other disturbance.Assumption 4.
τ
1 is added to consider the input delay of the system. In this work, the delay is removed since we only would like to investigate the development of control to tackle multiple constraint problem. The method to cope with input delay is very trivial and many state-of-the-art works can be found [27–29].Apart from the remarks and assumption above, we involve another variableTdg2=T˙g2. Therefore, the system model of (12) can be rewritten as(13)T˙g2=Tdg2MTdg2,Tg2,tT˙dg2+CTdg2,Tg2,tTdg2+gTg2,t=Vg1t+dTdg2,Tg2,t.Assumption 5.
In most literatures, only one primary/secondary pipe network is considered; in other words, the above model is an SISO system. In this research, we extend the system to an MIMO system which can generate different water temperatures in secondary pipe network to serve different operating conditions.With Assumption5, the input and output variables become Tg2=[Tg21,Tg22]T and Vg1=[Vg11,Vg12]T. Herein, without loss of generality, a two-dimension variable for input and output is considered. Moreover, we slightly abuse the denotation in model (13) by directly defining the original SISO system formulation as MIMO system without any additional denotation redefinition.In practical implementation, the valve opening degree must have an available range such that the rate of inflow is limited under a certain threshold. With this statement, we need to consider a nonlinear input saturation for the opening degreeVg1.(14)Vg1iVg1i0=V¯g1,ifVg1i0>V¯g1Vg1i0,ifV_g1≤Vg1i0≤V¯g1V_g1,ifVg1i0<V_g1,i=1,2,where Vg1i0 is the designed control law. V_g1 and V¯g1 are the lower bound and upper bound for the control input. Since the minimum value of valve opening degree is the situation when it completely turns off, the lower bound V_g1 is set as zero.Assumption 6.
The temperatureTg2 in secondary pipe network can be measured and its time derivative Tdg2 can be estimated with temperature variation between consecutive sampling intervals. Moreover, the desired tracking trajectory of the temperature Tg2d which is determined by the previous load prediction section and its nth (n=1,2) order derivatives are piecewise continuous, known, and bounded.In this tracking control, the temperature in secondary pipe network is investigated and the tracking errorTg2-Tg2d is used for backstepping design. Besides, since we would like to more precisely control the temperature inside the target space, the output error limitation of ±Tlim is involved, which describes the tracking error constraints to the control system as [-Tlim,-Tlim]T≤Tg2-Tg2d≤[Tlim,Tlim]T. For conciseness, the constraint is rewritten as -NTlim≤Tg2-Tg2d≤NTlim. In this paper, we apply the symmetric barrier Lyapunov function (SBLF) [30] incorporated with backstepping method to tackle this output constraint problem.Step 1.
Define error variablesz1=Tg2-Tg2d and z2=Tdg2-α1, where α1 is a stabilizing function that is to be designed. Consider the derivative of z1.(15)z˙1=T˙g2-T˙g2d=z2+α1-T˙g2d.Assumption 7.
The initial condition ofz1 satisfies z10<NTlim.Choose a positive definite SBLF candidate as(16)Vc1=∑i=1212log⁡NTlim2iNTlim2i-z1i=12log⁡NTlimTIxNTlimNTlimTIxNTlim-z1TIxz1+12log⁡NTlimTIyNTlimNTlimTIyNTlim-z1TIyz1,where(17)Ix=1000,Iy=0001.The selected SBLF candidate (16) is positive definite and C1 continues for -NTlim≤z1≤NTlim. The time derivative of Vc1 is(18)V˙c1=z1TIxz˙1NTlimTIxNTlim-z1TIxz1+z1TIyz˙1NTlimTIyNTlim-z1TIyz1. In (15), the stabilizing function α1 is designed as(19)α1=-NTlimTNTlim-z1Tz1Kc1z1+T˙g2d, where Kc1∈R2×2 is a symmetric positive definite design parameter matrix. With above design, z˙1 becomes(20)z˙1=z2-NTlimTNTlim-z1Tz1Kc1z1.The derivative of Vc1 gives(21)V˙c1=-2z1TKc1z1+z1TIxz2NTlimTIxNTlim-z1TIxz1+z1TIyz2NTlimTIyNTlim-z1TIyz1.Step 2.
To analyze the input saturation, the following auxiliary system is introduced to the system [12].(22)ζ˙=-Kζ1ζ-z2TΔVg1+0.5ΔVg1TΔVg1ζ2ζ+ΔVg1,ζ≥μζ0,ζ<μζ,where ζ∈R2 is the state of auxiliary system. ΔVg1=Vg1-Vg10. Kζ1 is a symmetric positive definite matrix. μζ is a small positive design parameter. For the temperature variation velocity of the system, no limitation is needed. Thus, a quadratic form Lyapunov function is selected as(23)Vc2∗=Vc1+12ζTζ+12z2TMz2. The time derivative of Vc2∗ is given by(24)V˙c2∗=-2z1TKc1z1+z1TIxz2NTlimTIxNTlim-z1TIxz1+z1TIyz2NTlimTIyNTlim-z1TIyz1+ζTζ˙+z2TMz˙2.The time derivative of z2 is(25)z˙2=T˙dg2-α˙1.By substituting z˙2 into (13), we get the following expression.(26)MT˙dg2=Mz˙2+α˙1=-CTdg2,Tg2-gTg2+Vg1+dTdg2,Tg2.Further, it yields(27)Mz˙2=-CTdg2,Tg2-gTg2+Vg1+dTdg2,Tg2-Mα˙1.Substituting (27) into (24),(28)V˙c2∗=-2z1TKc1z1+z1TIxz2NTlimTIxNTlim-z1TIxz1+z1TIyz2NTlimTIyNTlim-z1TIyz1-ζTKζ-12Iζ-z2TΔτ+z2T-CTdg2,Tg2-gTg2+Vg1+dTdg2,Tg2-Mα˙1.The desired control law is designed as follows:(29)Vg10=CTdg2,Tg2+gTg2-dTdg2,Tg2+Mα˙1-Kc2z2-ζ-Ixz1NTlimTIxNTlim-z1TIxz1-Iyz1NTlimTIyNTlim-z1TIyz1,where Kc2=Kc2T∈R2×2. To approximate the unknown terms ρc(Zc)=C(Tdg2,Tg2)+g(Tg2)-d(Tdg2,Tg2)+Mα˙1, a neural network is used.(30)W^cSZc=ρ^cZcwith(31)Wc∗SZc+ϵc=ρcZc, where W^c, Wc∗, ϵc are the weights, optimal weights, and approximation error, respectively. Zc=[Tg2T,Tdg2T,α1T,α˙1T]T is the input vector to the neural network. With this approximation, the control law Vg10 is proposed as(32)Vg10=W^cSZc-Kc2z2-ζ-Ixz1NTlimTIxNTlim-z1TIxz1-Iyz1NTlimTIyNTlim-z1TIyz1-NTlimTIxNTlimz2T+z1TKc1z1NTlimTIxNTlim-z1TIxz1-NTlimTIyNTlimz2T+z1TKc1z1NTlimTIyNTlim-z1TIyz1,where z2T+ denotes the Moore-Penrose pseudoinverse of z2T. Substituting (32) into (28), V˙c2∗ becomes(33)V˙c2∗≤-2z1Kc1z1-ζTKζ-12Iζ-z2TKc2z2-ζ-z2Tϵc+z2TW~cSZc-NTlimTIxNTlimz2Tz2T+z1TKc1z1NTlimTIxNTlim-z1TIxz1-NTlimTIyNTlimz2Tz2T+z1TKc1z1NTlimTIyNTlim-z1TIyz1≤-ζTKζ-12Iζ-z2TKc2z2-ζ-z2Tϵc+z2TW~cSZc-λminKc1Ξx+Uxz1Tz1Ξx-λminKc1Ξy+Uyz1Tz1Ξy≤-ζTKζ-12Iζ-z2TKc2z2-ζ-z2Tϵc+z2TW~cSZc-λminKc1NTlimTIxNTlimz1TIxz1Ξx-λminKc1NTlimTIyNTlimz1TIyz1Ξy,where(34)Ξx=NTlimTIxNTlim-z1TIxz1,Ξy=NTlimTIyNTlim-z1TIyz1,Ux=NTlimTIxNTlimz2Tz2T+,Uy=NTlimTIyNTlimz2Tz2T+ and λmin(•), λmax(•) denote the minimum and maximum eigenvalues of •. Lyapunov function candidate Vc2∗ is augmented with NN deviations to ensure the convergence of the neural weights.(35)Vc2=Vc1+12ζTζ+12z2TMz2+12∑i=12W~ciΓci-1W~ciT. The updating law of the NN is designed as(36)W^˙ci=-ΓciSiZcz2i+θ2iW^ci,where Γci>0 and θ2i>0(i=1,2) are tuning learning rate. z2i denotes the ith element of z2. Incorporate the following inequalities:(37)z2TKc2ζ≤σcz2Tz2+σc-1ζTKc2TKc2ζ,-z2Tϵc≤12z2Tz2+12ϵc2,-W~ciTW^ci≤-12W~ci2-12Wci∗2.The derivative of V˙c2 yields(38)V˙c2≤-λminKc1NTlimTIxNTlimz1TIxz1Ξxx-λminKc1NTlimTIyNTlimz1TIyz1Ξy-ζTKζ-12I-σc-1Kc2TKc2ζ-z2TKc2-σc+0.5Iz2-∑i=12θ2i2W~ci2+∑i=12θ2i2W~ci∗2+12ϵc2.The following inequalities hold if the output constraint is met, that is, -NTlim≤z1≤NTlim [31].(39)-Ixz1TIxz1IxNTlimTTIxNTlim-Ixz1TIxz1≤-logIxNTlimTTIxNTlimIxNTlimTTIxNTlim-Ixz1TIxz1-Iyz1TIyz1IyNTlimTTIyNTlim-Iyz1TIyz1≤-log⁡IyNTlimTTIyNTlimIyNTlimTTIyNTlim-Iyz1TIyz1.Invoking (39), (38) becomes(40)V˙c2≤-λminKc1NTlimTIxNTlimlog⁡NTlimTIxNTlimΞx-λminKc1NTlimTIyNTlimlog⁡NTlimTIyNTlimΞy-ζTKζ-12I-σc-1Kc2TKc2ζ-z2TKc2-σc+0.5Iz2-12∑i=12θ2iλmaxΓci-1W~ciΓci-1W~ciT+∑i=12θ2i2W~ci∗2+12ϵc2≤-αVc2+β,where α,β>0 and α=min[2λmin(Kc1)NTlimTIxNTlim,2λmin(Kc1)NTlimTIyNTlim, 2λmin(Kζ-1/2I-σc-1Kc2TKc2),2λmin(Kc2-(σc+0.5)I)/λmax(M),θ2i/λmax(Γci-1)(i=1,2)], β=∑i=12‍θ2i/2W~ci∗2+1/2ϵc2.Lemma 8 (see [32, 33]).
Given bounded initial conditions, if there exists aC1 continuous and positive definite Lyapunov function V(x) satisfying v1(x)≤V(x)≤v2(x), such that V˙≤-αV(x)+β, where v1,v2:Rn→R are class K functions and α,β>0, then the solution x(t) is uniformly bounded.Theorem 9.
Considering the central heating system (13) satisfying Assumptions 4–7, with the auxiliary analysis system (22), neural network weight update laws (36), and control input (32) (14), signals z1, z2, ζ, W~c are semi-globally uniformly bounded.Proof.
Incorporating Lemma8 and rigorous Lyapunov stability analysis (15)–(40), the Semi-Globally Uniform Boundedness (SGUB) of z1, z2, ζ, W~c is proven. The tracking error is limited as z1≤NTlim. Moreover, it should be noted that the above proof merely considers the result when the states of auxiliary system (22) satisfy ζ≥μζ which indicates the case where input saturation occurs. If there is no input saturation, that is, ζ<μζ, the control gives Vg1=Vg10. Since Vg10 is bounded, Vg1 must be bounded in this situation. In this fashion, Theorem 9 still holds. The proof is completed.Remark 10.
The research in this work is motivated by intelligently controlling a city central heating system to more rapidly satisfy user’s demand. But the proposed method can be tailored to other industry systems with similar formulation.


## 3.1. Preliminaries on NN Approximation

Radial Basis Function Neural Network (RBFNN) has been proven to have excellent function approximation and learning capabilities. Given a continuous nonlinear functionϕ(z):Rn→Rm on a compact set Ωz, it can be approximated by RBFNN in the following expression:(7)ϕNNz=WTSz,where z∈Ωz⊂Rn is the input vector. W∈Rp×m denotes the adjustable weight matrix; p stands for the neuron number. S(z)=[s1(z),…,sp(z)]T represents the basis function vector with si(z)=exp[-(z-μi)T(z-μi)/ψi2], i=1,2,…,p. μi=[μi1,μi2,…,μin] is the center of respective field; ψi is the width of Gaussian function. More accurately, (7) can be rewritten as(8)ϕz=W∗TSz+εz,where W∗ denotes the ideal weight; ε(z)∈Rm represents the approximation error and satisfies ε(z)≤σ, where σ is a positive constant. The ideal weight matrix W∗ is defined as(9)W∗≔argminW∈Rp×m⁡supz∈Ωzϕz-WSz.W∗ is merely an “artificial" quantity for analysis purpose. In control development, it needs to be estimated by W with designed updating law [25]. Additionally, NN approximation error ε(z) indicates the minimum deviation between optimal approximation solution W∗S(z) and the unknown function ϕ(z). It has been demonstrated that NN approximation can achieve any desired accuracy if the number of neurons is large enough; in other words, ε(z) can be arbitrarily small if p is sufficiently large [25].


## 3.2. Adaptive Control with Multiple Constraints

With the predicted supply water temperatureTg2 in secondary network, we intend to control the valve opening degree in primary network to ensure the supply water temperature tracking the prediction value. Herein, the prediction value Tg2 in the previous section acts as the desired signal denoted as Tg2d. The transfer function in quality adjustment is presented as follows [26].(10)G1s=K1T1s+1T2s+1e-τ1s, where K1 is a scalar coefficient, T1 and T2 are inertial time constant, and τ1 is the lagging time. This model describes the relationship between system input Vg1, that is, opening degree of the valve in primary pipe network, and the system output Tg2, that is, actual water temperature in secondary pipe network. To proceed, we transfer the above expression into time-domain formulation.(11)T1T2K1T¨g2+T1+T2K1T˙g2+Tg2K1=Vg1t-τ1.Equation (11) can be concisely expressed with(12)MT¨g2+CT˙g2+gTg2=Vg1t-τ1,where M=T1T2/K1, C=T1+T2/K1, g(Tg2)=Tg2/K1. Without loss of generality, we make the following remarks and assumptions.Remark 2.
Although the parametersK1,T1, and T2 in [26] are assumed as constants, it is obvious and straightforward that the transfer process is time-varying; that is, the aforementioned parameters are state-dependent and time-varying. Subsequently, M,C, and g(Tg2) should be more reasonably depicted as M(T˙g2,Tg2,t),C(T˙g2,Tg2,t), and g(Tg2,t).Remark 3.
A bounded disturbance term should be complemented in model (12) to consider the system noise and other disturbance.Assumption 4.
τ
1 is added to consider the input delay of the system. In this work, the delay is removed since we only would like to investigate the development of control to tackle multiple constraint problem. The method to cope with input delay is very trivial and many state-of-the-art works can be found [27–29].Apart from the remarks and assumption above, we involve another variableTdg2=T˙g2. Therefore, the system model of (12) can be rewritten as(13)T˙g2=Tdg2MTdg2,Tg2,tT˙dg2+CTdg2,Tg2,tTdg2+gTg2,t=Vg1t+dTdg2,Tg2,t.Assumption 5.
In most literatures, only one primary/secondary pipe network is considered; in other words, the above model is an SISO system. In this research, we extend the system to an MIMO system which can generate different water temperatures in secondary pipe network to serve different operating conditions.With Assumption5, the input and output variables become Tg2=[Tg21,Tg22]T and Vg1=[Vg11,Vg12]T. Herein, without loss of generality, a two-dimension variable for input and output is considered. Moreover, we slightly abuse the denotation in model (13) by directly defining the original SISO system formulation as MIMO system without any additional denotation redefinition.In practical implementation, the valve opening degree must have an available range such that the rate of inflow is limited under a certain threshold. With this statement, we need to consider a nonlinear input saturation for the opening degreeVg1.(14)Vg1iVg1i0=V¯g1,ifVg1i0>V¯g1Vg1i0,ifV_g1≤Vg1i0≤V¯g1V_g1,ifVg1i0<V_g1,i=1,2,where Vg1i0 is the designed control law. V_g1 and V¯g1 are the lower bound and upper bound for the control input. Since the minimum value of valve opening degree is the situation when it completely turns off, the lower bound V_g1 is set as zero.Assumption 6.
The temperatureTg2 in secondary pipe network can be measured and its time derivative Tdg2 can be estimated with temperature variation between consecutive sampling intervals. Moreover, the desired tracking trajectory of the temperature Tg2d which is determined by the previous load prediction section and its nth (n=1,2) order derivatives are piecewise continuous, known, and bounded.In this tracking control, the temperature in secondary pipe network is investigated and the tracking errorTg2-Tg2d is used for backstepping design. Besides, since we would like to more precisely control the temperature inside the target space, the output error limitation of ±Tlim is involved, which describes the tracking error constraints to the control system as [-Tlim,-Tlim]T≤Tg2-Tg2d≤[Tlim,Tlim]T. For conciseness, the constraint is rewritten as -NTlim≤Tg2-Tg2d≤NTlim. In this paper, we apply the symmetric barrier Lyapunov function (SBLF) [30] incorporated with backstepping method to tackle this output constraint problem.Step 1.
Define error variablesz1=Tg2-Tg2d and z2=Tdg2-α1, where α1 is a stabilizing function that is to be designed. Consider the derivative of z1.(15)z˙1=T˙g2-T˙g2d=z2+α1-T˙g2d.Assumption 7.
The initial condition ofz1 satisfies z10<NTlim.Choose a positive definite SBLF candidate as(16)Vc1=∑i=1212log⁡NTlim2iNTlim2i-z1i=12log⁡NTlimTIxNTlimNTlimTIxNTlim-z1TIxz1+12log⁡NTlimTIyNTlimNTlimTIyNTlim-z1TIyz1,where(17)Ix=1000,Iy=0001.The selected SBLF candidate (16) is positive definite and C1 continues for -NTlim≤z1≤NTlim. The time derivative of Vc1 is(18)V˙c1=z1TIxz˙1NTlimTIxNTlim-z1TIxz1+z1TIyz˙1NTlimTIyNTlim-z1TIyz1. In (15), the stabilizing function α1 is designed as(19)α1=-NTlimTNTlim-z1Tz1Kc1z1+T˙g2d, where Kc1∈R2×2 is a symmetric positive definite design parameter matrix. With above design, z˙1 becomes(20)z˙1=z2-NTlimTNTlim-z1Tz1Kc1z1.The derivative of Vc1 gives(21)V˙c1=-2z1TKc1z1+z1TIxz2NTlimTIxNTlim-z1TIxz1+z1TIyz2NTlimTIyNTlim-z1TIyz1.Step 2.
To analyze the input saturation, the following auxiliary system is introduced to the system [12].(22)ζ˙=-Kζ1ζ-z2TΔVg1+0.5ΔVg1TΔVg1ζ2ζ+ΔVg1,ζ≥μζ0,ζ<μζ,where ζ∈R2 is the state of auxiliary system. ΔVg1=Vg1-Vg10. Kζ1 is a symmetric positive definite matrix. μζ is a small positive design parameter. For the temperature variation velocity of the system, no limitation is needed. Thus, a quadratic form Lyapunov function is selected as(23)Vc2∗=Vc1+12ζTζ+12z2TMz2. The time derivative of Vc2∗ is given by(24)V˙c2∗=-2z1TKc1z1+z1TIxz2NTlimTIxNTlim-z1TIxz1+z1TIyz2NTlimTIyNTlim-z1TIyz1+ζTζ˙+z2TMz˙2.The time derivative of z2 is(25)z˙2=T˙dg2-α˙1.By substituting z˙2 into (13), we get the following expression.(26)MT˙dg2=Mz˙2+α˙1=-CTdg2,Tg2-gTg2+Vg1+dTdg2,Tg2.Further, it yields(27)Mz˙2=-CTdg2,Tg2-gTg2+Vg1+dTdg2,Tg2-Mα˙1.Substituting (27) into (24),(28)V˙c2∗=-2z1TKc1z1+z1TIxz2NTlimTIxNTlim-z1TIxz1+z1TIyz2NTlimTIyNTlim-z1TIyz1-ζTKζ-12Iζ-z2TΔτ+z2T-CTdg2,Tg2-gTg2+Vg1+dTdg2,Tg2-Mα˙1.The desired control law is designed as follows:(29)Vg10=CTdg2,Tg2+gTg2-dTdg2,Tg2+Mα˙1-Kc2z2-ζ-Ixz1NTlimTIxNTlim-z1TIxz1-Iyz1NTlimTIyNTlim-z1TIyz1,where Kc2=Kc2T∈R2×2. To approximate the unknown terms ρc(Zc)=C(Tdg2,Tg2)+g(Tg2)-d(Tdg2,Tg2)+Mα˙1, a neural network is used.(30)W^cSZc=ρ^cZcwith(31)Wc∗SZc+ϵc=ρcZc, where W^c, Wc∗, ϵc are the weights, optimal weights, and approximation error, respectively. Zc=[Tg2T,Tdg2T,α1T,α˙1T]T is the input vector to the neural network. With this approximation, the control law Vg10 is proposed as(32)Vg10=W^cSZc-Kc2z2-ζ-Ixz1NTlimTIxNTlim-z1TIxz1-Iyz1NTlimTIyNTlim-z1TIyz1-NTlimTIxNTlimz2T+z1TKc1z1NTlimTIxNTlim-z1TIxz1-NTlimTIyNTlimz2T+z1TKc1z1NTlimTIyNTlim-z1TIyz1,where z2T+ denotes the Moore-Penrose pseudoinverse of z2T. Substituting (32) into (28), V˙c2∗ becomes(33)V˙c2∗≤-2z1Kc1z1-ζTKζ-12Iζ-z2TKc2z2-ζ-z2Tϵc+z2TW~cSZc-NTlimTIxNTlimz2Tz2T+z1TKc1z1NTlimTIxNTlim-z1TIxz1-NTlimTIyNTlimz2Tz2T+z1TKc1z1NTlimTIyNTlim-z1TIyz1≤-ζTKζ-12Iζ-z2TKc2z2-ζ-z2Tϵc+z2TW~cSZc-λminKc1Ξx+Uxz1Tz1Ξx-λminKc1Ξy+Uyz1Tz1Ξy≤-ζTKζ-12Iζ-z2TKc2z2-ζ-z2Tϵc+z2TW~cSZc-λminKc1NTlimTIxNTlimz1TIxz1Ξx-λminKc1NTlimTIyNTlimz1TIyz1Ξy,where(34)Ξx=NTlimTIxNTlim-z1TIxz1,Ξy=NTlimTIyNTlim-z1TIyz1,Ux=NTlimTIxNTlimz2Tz2T+,Uy=NTlimTIyNTlimz2Tz2T+ and λmin(•), λmax(•) denote the minimum and maximum eigenvalues of •. Lyapunov function candidate Vc2∗ is augmented with NN deviations to ensure the convergence of the neural weights.(35)Vc2=Vc1+12ζTζ+12z2TMz2+12∑i=12W~ciΓci-1W~ciT. The updating law of the NN is designed as(36)W^˙ci=-ΓciSiZcz2i+θ2iW^ci,where Γci>0 and θ2i>0(i=1,2) are tuning learning rate. z2i denotes the ith element of z2. Incorporate the following inequalities:(37)z2TKc2ζ≤σcz2Tz2+σc-1ζTKc2TKc2ζ,-z2Tϵc≤12z2Tz2+12ϵc2,-W~ciTW^ci≤-12W~ci2-12Wci∗2.The derivative of V˙c2 yields(38)V˙c2≤-λminKc1NTlimTIxNTlimz1TIxz1Ξxx-λminKc1NTlimTIyNTlimz1TIyz1Ξy-ζTKζ-12I-σc-1Kc2TKc2ζ-z2TKc2-σc+0.5Iz2-∑i=12θ2i2W~ci2+∑i=12θ2i2W~ci∗2+12ϵc2.The following inequalities hold if the output constraint is met, that is, -NTlim≤z1≤NTlim [31].(39)-Ixz1TIxz1IxNTlimTTIxNTlim-Ixz1TIxz1≤-logIxNTlimTTIxNTlimIxNTlimTTIxNTlim-Ixz1TIxz1-Iyz1TIyz1IyNTlimTTIyNTlim-Iyz1TIyz1≤-log⁡IyNTlimTTIyNTlimIyNTlimTTIyNTlim-Iyz1TIyz1.Invoking (39), (38) becomes(40)V˙c2≤-λminKc1NTlimTIxNTlimlog⁡NTlimTIxNTlimΞx-λminKc1NTlimTIyNTlimlog⁡NTlimTIyNTlimΞy-ζTKζ-12I-σc-1Kc2TKc2ζ-z2TKc2-σc+0.5Iz2-12∑i=12θ2iλmaxΓci-1W~ciΓci-1W~ciT+∑i=12θ2i2W~ci∗2+12ϵc2≤-αVc2+β,where α,β>0 and α=min[2λmin(Kc1)NTlimTIxNTlim,2λmin(Kc1)NTlimTIyNTlim, 2λmin(Kζ-1/2I-σc-1Kc2TKc2),2λmin(Kc2-(σc+0.5)I)/λmax(M),θ2i/λmax(Γci-1)(i=1,2)], β=∑i=12‍θ2i/2W~ci∗2+1/2ϵc2.Lemma 8 (see [32, 33]).
Given bounded initial conditions, if there exists aC1 continuous and positive definite Lyapunov function V(x) satisfying v1(x)≤V(x)≤v2(x), such that V˙≤-αV(x)+β, where v1,v2:Rn→R are class K functions and α,β>0, then the solution x(t) is uniformly bounded.Theorem 9.
Considering the central heating system (13) satisfying Assumptions 4–7, with the auxiliary analysis system (22), neural network weight update laws (36), and control input (32) (14), signals z1, z2, ζ, W~c are semi-globally uniformly bounded.Proof.
Incorporating Lemma8 and rigorous Lyapunov stability analysis (15)–(40), the Semi-Globally Uniform Boundedness (SGUB) of z1, z2, ζ, W~c is proven. The tracking error is limited as z1≤NTlim. Moreover, it should be noted that the above proof merely considers the result when the states of auxiliary system (22) satisfy ζ≥μζ which indicates the case where input saturation occurs. If there is no input saturation, that is, ζ<μζ, the control gives Vg1=Vg10. Since Vg10 is bounded, Vg1 must be bounded in this situation. In this fashion, Theorem 9 still holds. The proof is completed.Remark 10.
The research in this work is motivated by intelligently controlling a city central heating system to more rapidly satisfy user’s demand. But the proposed method can be tailored to other industry systems with similar formulation.


## 4. Simulation Study

In this simulation, a numerical study is conducted to demonstrate the validation of proposed load prediction and control scheme.


### 4.1. NN Based Heating Load Prediction

In this study, the historical recorded data from a block exchanger station in local city is collected as the training and testing samples. The period is half month from December 5, 2014 to December 20, 2014. Among them, data from 5 to 19 December are employed as training data. Information of 20 December acts as testing samples. As mentioned in Section2.2.1, four-dimension features are imported into the approximation network, that is, time, outdoor temperature, solar radiation intensity, and Tg2 in previous day. Tg2 in current day is regarded as ground truth to produce the prediction error for weights updating.As for the NN architecture, a 4-9-1 layout is selected with learning rate fixed at 0.2, momentum factor 0.8, maximum training iteration 1500, and maximum square error10-2. Gradient descent method with momentum [34] is employed as training function. In this manner, NN training can be guaranteed with a fast convergence speed and avoid local minima problem.After training the network, data from 20 December is adopted to test the predictor. The testing input is plotted in Figure3. Actually, this testing data is very challenging; an obvious temperature decreasing can be observed in Figure 3(a). In the midnight of 19 December, that is, when the time is zero, outdoor temperature is as high as −8°C. But during the morning, a sudden temperature decrease occurs at 6–8 am and the lowest temperature can be less than −18°C. At the midnight of 20 December (23-24), the temperature is as low as −20°C. This severe temperature variation brings additional challenges to the prediction task since the input of previous day’s Tg2 (Figure 3(c)) cannot provide an efficient hint or reference to the prediction. Robustness and generalization capability of the proposed algorithm are required to tackle this situation. The prediction results are shown in Figure 4.Figure 3
Test data plots on December 20th.
(a)
Outdoor temperature on Dec 20th
(b)
Solar radiation intensity on Dec 20th
(c)
T
g
2 record on Dec 19thFigure 4
NN prediction performance.
(a)
Actual and predictedTg2 on Dec 20th
(b)
Prediction error rate
(c)
Heating load prediction resultsFigure4(a) reveals that the NN predicted Tg2 can generally follow the actual requirement with acceptable accuracy under such harsh temperature variation condition. The prediction error rate in Figure 4(b) demonstrates that the error rate can be guaranteed within 2%. It should be mentioned that although the output from proposed predictor cannot perfectly match the actual data, the accuracy is still sufficient to ensure a good heating service since there always exists allowable range or bias in practical engineering problem so that values within the range can result in acceptable performance. The load prediction is shown in Figure 4(c) with mathematical deduction using Tg2 as introduced in Section 2.2.1.


### 4.2. Adaptive NN Control

In this subsection, the performance of proposed control scheme is investigated. As aforementioned, we consider a more complicated nonlinear system model. Some of the system parameters refer to [26]. Without loss of generality, we make M(Tdg2,Tg2,t) as a constant matrix and C(Tdg2,Tg2,t), g(Tg2), t as time-varying and nonlinear. The specifications of these matrices are M=diag([149.27,132.7]), C=diag([26.24∗Tdg11,23.26∗Tdg22]), g=0.8197∗Tg2, and d=[0.05∗Tdg21∗cos⁡0.5t,0.05∗Tdg22∗cos⁡(0.5t+π/2)]T. The overall simulation time is 40 s; the desired secondary network temperature Tg2d is assumed to be [47.7+(50-47.7)∗sin(π/2∗40t+π/2),45.6+(48-45.6)∗sin(π/2∗40t+π/2)]T which depicts a temperature variation from 50/48 to 47.7/45.6 for the two subnetworks, respectively. Systematically, the predicted result Tg2 in previous subsection is supposed to be the reference signal Tg2d in the control task. But, for more efficient and concise demonstration purpose, a simplified tracking trajectory is designed instead as mentioned in afore contents. The input and output constraints are defined as V¯g1=65, V_g1=0, and NTlim=0.008. It should be noted that all the simulation setting parameters above can be selected arbitrarily based on the practical implementation requirements. The values in this paper are merely a test case to verify the control performance.As for the proposed NN control settings, the coefficientsKc1, Kc2, Kϵ1, θ21, θ22, μϵ are tuned as diag([750,720]), diag([30.8,33]), diag([2,1]), 1, 1, 0.005, respectively. The learning rate is Γc1=Γc2=5. The number of nodes in the proposed control is eight which is obtained through trial and error in the simulation. To achieve comparative study, a PID control is implemented with properly tuned parameters Kp=diag([8100,6750]), Kd=diag([50,50]), Ki=diag([120,90]). The simulation results are shown in Figures 5–9.Figure 5
Tracking performance forTg2 with proposed control.
(a)
Tracking temperatureTg21 in subnetwork 1
(b)
Tracking temperatureTg22 in subnetwork 2Figure 6
Tracking performance forTg2 with PID control.
(a)
Tracking temperatureTg21 in subnetwork 1
(b)
Tracking temperatureTg22 in subnetwork 2Figure 7
Tracking error comparative study.
(a)
Tracking error with proposed control
(b)
Tracking error with PID controlFigure 8
Control input comparative study.
(a)
InputVg1 with proposed control
(b)
InputVg1 with PID controlFigure 9
Temperature variation speed comparative study.
(a)
Temperature variation speed with proposed control
(b)
Temperature variation speed with PID controlFigures5–7 reveal that although both the proposed and PID control can achieve acceptable performance as the first two figures show, PID control will result in larger tracking errors that exceed the constraint 0.008. On the contrary, the errors with proposed control can be well limited within the constraints throughout the whole simulation process. As for the control input in Figure 8, attributing to the auxiliary system (22), the control effort is successfully guaranteed below the input constraint, that is, 65 without compromising the system stability and tracking performance as presented in Figure 8(a). By comparison, control input with PID control can be as high as 72 which goes beyond the constraint as shown in Figure 8(b). Moreover, if we would like to further reduce the tracking error of PID control, the parameters Kp, Ki, Kd will increase accordingly, which subsequently will cause the increase of control input. Therefore, PID control cannot achieve the control objective in the presence of multiple constraints and system nonlinearity. Figure 9 compares the temperature variation speed of the two control schemes. Both of them can ensure the boundedness of all the closed-loop signals in this case study. Moreover, system stability with proposed control is also determined by the selection of constraints. If too tight constraints are involved, the system will probably goes instable. Thus, a good tradeoff between the constraint values should be guaranteed to ensure a safe and smooth operation of the control system. Finally, it can be observed that the output constraint in this simulation is very tight (NTlim=0.008). In reality, due to the high nonlinearity of actual system, the control difficulty increases and this constraint should be looser. In that case, the proposed control (nonlinear control in essence) will show more obvious advantage over traditional PID scheme because of the specific designed constraint-handle mechanism.


## 4.1. NN Based Heating Load Prediction

In this study, the historical recorded data from a block exchanger station in local city is collected as the training and testing samples. The period is half month from December 5, 2014 to December 20, 2014. Among them, data from 5 to 19 December are employed as training data. Information of 20 December acts as testing samples. As mentioned in Section2.2.1, four-dimension features are imported into the approximation network, that is, time, outdoor temperature, solar radiation intensity, and Tg2 in previous day. Tg2 in current day is regarded as ground truth to produce the prediction error for weights updating.As for the NN architecture, a 4-9-1 layout is selected with learning rate fixed at 0.2, momentum factor 0.8, maximum training iteration 1500, and maximum square error10-2. Gradient descent method with momentum [34] is employed as training function. In this manner, NN training can be guaranteed with a fast convergence speed and avoid local minima problem.After training the network, data from 20 December is adopted to test the predictor. The testing input is plotted in Figure3. Actually, this testing data is very challenging; an obvious temperature decreasing can be observed in Figure 3(a). In the midnight of 19 December, that is, when the time is zero, outdoor temperature is as high as −8°C. But during the morning, a sudden temperature decrease occurs at 6–8 am and the lowest temperature can be less than −18°C. At the midnight of 20 December (23-24), the temperature is as low as −20°C. This severe temperature variation brings additional challenges to the prediction task since the input of previous day’s Tg2 (Figure 3(c)) cannot provide an efficient hint or reference to the prediction. Robustness and generalization capability of the proposed algorithm are required to tackle this situation. The prediction results are shown in Figure 4.Figure 3
Test data plots on December 20th.
(a)
Outdoor temperature on Dec 20th
(b)
Solar radiation intensity on Dec 20th
(c)
T
g
2 record on Dec 19thFigure 4
NN prediction performance.
(a)
Actual and predictedTg2 on Dec 20th
(b)
Prediction error rate
(c)
Heating load prediction resultsFigure4(a) reveals that the NN predicted Tg2 can generally follow the actual requirement with acceptable accuracy under such harsh temperature variation condition. The prediction error rate in Figure 4(b) demonstrates that the error rate can be guaranteed within 2%. It should be mentioned that although the output from proposed predictor cannot perfectly match the actual data, the accuracy is still sufficient to ensure a good heating service since there always exists allowable range or bias in practical engineering problem so that values within the range can result in acceptable performance. The load prediction is shown in Figure 4(c) with mathematical deduction using Tg2 as introduced in Section 2.2.1.


## 4.2. Adaptive NN Control

In this subsection, the performance of proposed control scheme is investigated. As aforementioned, we consider a more complicated nonlinear system model. Some of the system parameters refer to [26]. Without loss of generality, we make M(Tdg2,Tg2,t) as a constant matrix and C(Tdg2,Tg2,t), g(Tg2), t as time-varying and nonlinear. The specifications of these matrices are M=diag([149.27,132.7]), C=diag([26.24∗Tdg11,23.26∗Tdg22]), g=0.8197∗Tg2, and d=[0.05∗Tdg21∗cos⁡0.5t,0.05∗Tdg22∗cos⁡(0.5t+π/2)]T. The overall simulation time is 40 s; the desired secondary network temperature Tg2d is assumed to be [47.7+(50-47.7)∗sin(π/2∗40t+π/2),45.6+(48-45.6)∗sin(π/2∗40t+π/2)]T which depicts a temperature variation from 50/48 to 47.7/45.6 for the two subnetworks, respectively. Systematically, the predicted result Tg2 in previous subsection is supposed to be the reference signal Tg2d in the control task. But, for more efficient and concise demonstration purpose, a simplified tracking trajectory is designed instead as mentioned in afore contents. The input and output constraints are defined as V¯g1=65, V_g1=0, and NTlim=0.008. It should be noted that all the simulation setting parameters above can be selected arbitrarily based on the practical implementation requirements. The values in this paper are merely a test case to verify the control performance.As for the proposed NN control settings, the coefficientsKc1, Kc2, Kϵ1, θ21, θ22, μϵ are tuned as diag([750,720]), diag([30.8,33]), diag([2,1]), 1, 1, 0.005, respectively. The learning rate is Γc1=Γc2=5. The number of nodes in the proposed control is eight which is obtained through trial and error in the simulation. To achieve comparative study, a PID control is implemented with properly tuned parameters Kp=diag([8100,6750]), Kd=diag([50,50]), Ki=diag([120,90]). The simulation results are shown in Figures 5–9.Figure 5
Tracking performance forTg2 with proposed control.
(a)
Tracking temperatureTg21 in subnetwork 1
(b)
Tracking temperatureTg22 in subnetwork 2Figure 6
Tracking performance forTg2 with PID control.
(a)
Tracking temperatureTg21 in subnetwork 1
(b)
Tracking temperatureTg22 in subnetwork 2Figure 7
Tracking error comparative study.
(a)
Tracking error with proposed control
(b)
Tracking error with PID controlFigure 8
Control input comparative study.
(a)
InputVg1 with proposed control
(b)
InputVg1 with PID controlFigure 9
Temperature variation speed comparative study.
(a)
Temperature variation speed with proposed control
(b)
Temperature variation speed with PID controlFigures5–7 reveal that although both the proposed and PID control can achieve acceptable performance as the first two figures show, PID control will result in larger tracking errors that exceed the constraint 0.008. On the contrary, the errors with proposed control can be well limited within the constraints throughout the whole simulation process. As for the control input in Figure 8, attributing to the auxiliary system (22), the control effort is successfully guaranteed below the input constraint, that is, 65 without compromising the system stability and tracking performance as presented in Figure 8(a). By comparison, control input with PID control can be as high as 72 which goes beyond the constraint as shown in Figure 8(b). Moreover, if we would like to further reduce the tracking error of PID control, the parameters Kp, Ki, Kd will increase accordingly, which subsequently will cause the increase of control input. Therefore, PID control cannot achieve the control objective in the presence of multiple constraints and system nonlinearity. Figure 9 compares the temperature variation speed of the two control schemes. Both of them can ensure the boundedness of all the closed-loop signals in this case study. Moreover, system stability with proposed control is also determined by the selection of constraints. If too tight constraints are involved, the system will probably goes instable. Thus, a good tradeoff between the constraint values should be guaranteed to ensure a safe and smooth operation of the control system. Finally, it can be observed that the output constraint in this simulation is very tight (NTlim=0.008). In reality, due to the high nonlinearity of actual system, the control difficulty increases and this constraint should be looser. In that case, the proposed control (nonlinear control in essence) will show more obvious advantage over traditional PID scheme because of the specific designed constraint-handle mechanism.


## 5. Conclusion

In this paper, an NN based central heating system has been proposed which consists of load prediction and adaptive control modules. By employing proposed data-driven prediction, the utilization of complicated physical model can be avoided and response time to temperature variation shrinks. The proposed solution is accurate and concise. With the predicted load demand, an intelligent adaptive control is proposed to handle multiple constraints on temperature control errors and valve opening degree in practical implementation. NN is applied to approximate the complicated system dynamic model which allows users to achieve accurate control without exact knowledge of the system parameters. Boundedness of all the signals in the closed loop has been guaranteed via rigorous Lyapunov verification. Finally, numerical simulation study has been conducted and feasibility of the proposed framework is demonstrated.



---

*Source: 2908608-2018-02-07.xml*