International Journal of Automation and Computing  2018, Vol. 15 Issue (6): 707-715 PDF
Model Predictive Control for Vapor Compression Cycle of Refrigeration Process
Xiao-Hong Yin, Shao-Yuan Li
Department of Automation, Shanghai Jiao Tong University, and Key Laboratory of System Control and Information Processing, Ministry of Education of China, Shanghai 200240, China
Abstract: A model predictive controller based on a novel structure selection criterion for the vapor compression cycle (VCC) of refrigeration process is proposed in this paper.Firstly, those system variables are analyzed which exert significant influences on the system performance.Then the structure selection criterion, a trade-off between computation complexity and model performance, is applied to different model structures, and the results are utilized to determine the optimized model structure for controller design.The controller based on multivariable model predictive control (MPC) strategy is designed, and the optimization problem for the reduced order models is formulated as a constrained minimization problem.The effectiveness of the proposed MPC controller is verified on the experimental rig.
Key words: Vapor compression refrigeration cycle     model structure     structure selection criterion     model reduction     model predictive control
1 Introduction

To meet the ever-increasing demand for modeling accuracy, the models of large-scale systems are high-order, computationally complex, and with complex nonlinearities. However, these characteristics greatly limit the development of effective control strategies. In order to make a trade-off between the modeling accuracy and the convenience in controller design, model structure design plays an important role and draw intensive research interests over the last two decades[1, 2].

Vapor compression cycle (VCC) is a core element in the heating, ventilating, and air-conditioning (HVAC) systems. The VCC system is a high dimensional system that has obviously nonlinear thermodynamic coupling characteristics and time-varying dynamics characteristics. As a result, direct numerical simulation for such a large scale system becomes an intractable task. Through model reduction to choose appropriate model structure is an approach to overcome this problem[3]. It aims to approximate a large scale system by low dimensional models that have similar response characteristics as the original system so that the feasible controllers can be designed.

There are a number of systematic strategies for model reduction proposed in recent year. One of the most common model reduction schemes is balanced truncation which was applied in the complex systems by Moore (1981)[4]. Another popular model reduction method is the balanced residualization introduced by Fernando and Nicholson[5]. Furthermore, this approach was extended by Rasmussen[6] to the modeling and controller design of VCC system. The reduction results of these researches were to some extent simplified the model but still presented great challenges for the controller design and implementation. Based on Laguerre polynomials, Wang and Jiang[7] proposed a model reduction method for coupled systems in time-domain. According to Laguerre coefficients, projection matrices were defined, and low order coupled systems were generated to match a desired number of these coefficients. This method retained the stability of coupled systems, however, it was not applied to an actual HVAC system. A novel model reduction strategy is based on the proper orthogonal decomposition (POD) technique, with which a reduced order model was obtained by Guha and Mishra[8] verified by a nonlinear induction heating system.

For controller design, PI (proportional integration) or PID (proportion integration differentiation) feedback control algorithm is widely used in HVAC (heating ventilation and air conditionig) fields due to its simplicity[9-13]. However, it is difficult for the multivariable systems to overcome the coupling effects among each degree of freedom resulting in poor performance. Neural network control and robust control have been introduced to deal with nonlinearities or uncertainties in HVAC processes[14-16]. In [17], a gain-scheduling approach was utilized to handle the nonlinearity. The development of model predictive control (MPC) is considered as the recent major achievement in control literature, which has been widely accepted as the next generation of a practical control technology[18-19]. The application of MPC in HVAC systems can be found in the research of Elliott and Rasmussen[20]. He designed a two-input-two-output MPC controller and a single-input-single-output PI controller to control a multiple evaporator system. Xu et al.[21] presented a linear matrix inequality-based robust MPC strategy for the temperature control of an air-conditioning system. An offset-free MPC controller, comprising of a Luenberger observer, was implemented on a VCC system by Wallace to resolve the problem of plant-model mismatch[18]. Many studies presented different multivariable control strategies to improve the system energy efficiency; however, based on the premise of simplifying the controller complexity, most of them are based on the two-input-two-output model structure.

In this paper, a structure selection criterion, which can be used to evaluate the performance of different model structures and consequently to choose the optimal model structure for MPC controller design, is proposed. In order to select the impressionable variables which have main influence on the system performance, the relationships among variables of VCC system are analyzed first. According to the computational results under different low-order models, the optimized simplified model is determined. Then MPC based controller is designed for the optimized low-order model. Experiment results indicate that the proposed method has high modeling accuracy and great control performance.

The remainder of the paper is organized as follows. Section 2 details the dynamic model of the VCC system. Section 3 proposes a structure selection criterion which is used to evaluate the performance of different reduced order models and choose the optimal simplified model. Section 4 proposes the MPC based controller design. Section 5 presents the experimental results which justify some of the conclusions discussed in the previous sections. Section 6 summarizes the main conclusions.

2 Dynamic model of VCC system

A typical VCC system is composed of a compressor, a condenser, an expansion valve, and an evaporator. Fig. 1 illustrates the iterative cycle of VCC system. The refrigerant with high pressure and temperature enters the condenser and is cooled into subcooled liquid phase by heat exchange with cold fluid flowing across the coil or tubes. Then the condensed refrigerant enters the expansion valve where the pressure is reduced abruptly. The two-phase refrigerant at low pressure and temperature exits the expansion valve and routes through the evaporator, in which the refrigerant evaporates and absorbs the heat from ambient air. At the outlet of evaporator, the refrigerant is a vapor at superheat state, and then enters the compressor where it is further compressed into a superheated vapor which has high pressure and temperature. After the compression process, the circulating refrigerant is routed back into the condenser to complete the refrigeration cycle. From the point of view of energy consumption, the system can be described as Fig. 2. It displays the relationship between pressure and enthalpy which is characterized as the energy parameter of VCC system. The dotted line in Fig. 2 is the saturation curve of refrigerant. The condensation process and evaporation process are considered as isobaric processes, while the expansion process is considered as isenthalpic process.

 Download: larger image Fig. 1. The schematic diagram of VCC system

 Download: larger image Fig. 2. Pressure and enthalpy diagram of VCC system

Using the lumped-parameter and moving-boundary method, the dynamic model of each component of VCC was derived by Rasmussen and Alleyan[6], and is briefed as below.

Compressor. The dynamics of the compressor are considered to be much faster than those of heat exchangers, therefore, its mass flow rate can be modeled as a static component

 \begin{align} \label{eq1} \dot {m}_k =F_k V_k \rho _k \left( {1+C_k +D_k \left( {\frac{P_{ko} }{P_{ki} }} \right)^{\frac{1}{n}}} \right) \end{align} (1)

where $\dot {m}_k$ is the mass flow rate of the refrigerant through the compressor, $F_k$ is the compressor speed, $V_k$ is the effective displacement volume of the compressor, $C_k$ and $D_k$ are volumetric efficiency coefficients for the compressor, $n$ is the polytropic coefficient, $P_{ki}$ and $P_{ko}$ are the inlet pressure and outlet pressure across the compressor, respectively.

 \begin{align*} f_c \left( {x_c, u_c } \right)= \\ \left[{{\begin{array}{*{20}l} {\dot {m}_{ci} \left( {h_{ci}-h_{cg} } \right)+A_{ci} \alpha _{ci1} \displaystyle\frac{L_{c1} }{L_{cT} }\left( {T_{cw1}-T_{cr1} } \right)-\left( {\displaystyle\frac{1}{2}\left( {\displaystyle\frac{\partial \rho _{c1} }{\partial h_{c1} }} \right)(h_{c1}-h_{cg} )+\displaystyle\frac{1}{2}\rho_{c_1}} \right)A_c L_{c1} \dot {h}_{ci} } \\ {\dot {m}_{ci} h_{cg} -\dot {m}_{co} h_{cf} +\alpha _{ci2} A_{ci} \left( {\displaystyle\frac{L_{c2} }{L_{cT} }} \right)\left( {T_{cw2} -T_{cr2} } \right)-\left( {\displaystyle\frac{1}{2}\displaystyle\frac{\partial \rho _{c1} }{\partial h_{c1} }} \right)A_c h_{cg} L_{c1} \dot {h}_{ci} } \\ \!\!\!{\begin{array}{l} \dot {m}_{co} (h_{cf} -h_{co} )+A_{ci} \alpha _{ci3} \displaystyle\frac{L_3 }{L_{cT} }\left( {T_{cw3} -T_{cr3} } \right)\quad \\ \dot {m}_{ci} -\dot {m}_{co} -\left( {\displaystyle\frac{1}{2}\displaystyle\frac{\partial \rho _{c1} }{\partial h_{c1} }} \right)A_c L_{c1} \dot {h}_{ci} \\ \end{array}} \\ {\alpha _{co} A_{co} \left( {T_{ca} -T_{cw1} } \right)-\displaystyle\frac{\alpha _{ci1} A_{ci} \left( {T_{cw1} -T_{cr1} } \right)}{L_1 }} \\ \!\!\! {\begin{array}{l} \alpha _{co} A_{co} \left( {T_{ca} -T_{cw2} } \right)-\displaystyle\frac{\alpha _{ci2} A_{ci} \left( {T_{cw2} -T_{cr2} } \right)}{L_2 } \\ \! \alpha _{co} A_{co} \left( {T_{ca} -T_{cw3} } \right)-\displaystyle\frac{\alpha _{ci3} A_{ci} \left( {T_{cw3} -T_{cr3} } \right)}{L_{c3} } \\ \end{array}} \\ \end{array} }} \right]. \end{align*}

Condenser. According to the state of refrigerant, the condenser can be divided into three regions: a subcooled liquid region, a two-phase region and a superheated vapor region. Based on the conservation of the refrigerant mass and energy, the dynamic equations of condenser can be established in a nonlinear state space form with 7 states and 5 inputs, shown as (2)-(4):

 \begin{align} &Z_{c} \left( {x_c, u_c } \right)\cdot \dot {x}_c =f_c \left( {x_c, u_c } \right) \end{align} (2)
 \begin{align} &\label{eq3} x_c =\left[{{\begin{array}{*{20}l} {\dot {L}_{c1} } \hfill&{\dot {L}_{c2} } \hfill&{\dot {P}_c } \hfill & {\dot {h}_{cro} } \hfill&{\dot {T}_{cw1} } \hfill&{\dot {T}_{cw2} } \hfill&{\dot {T}_{cw3} } \hfill \\ \end{array} }} \right]^{\rm T} \end{align} (3)
 \begin{align} &\label{eq4} u_c =\left[{{\begin{array}{*{20}l} {\dot {m}_{ci} } \hfill&{\dot {m}_{co} } \hfill&{h_{cri} } \hfill & {T_{ca} } \hfill&{\dot {m}_{ca} } \hfill \\ \end{array} }} \right]^{\rm T} \end{align} (4)

where the state variables are: length of the two condensation regions $L_{c1}$ and $L_{c2}$; refrigerant pressure $P_c$; refrigerant enthalpy at outlet $h_{cro}$; the wall temperatures in the three regions $T_{cw1}$, $T_{cw2}$ and $T_{cw3}$, respectively. The input variables are: mass flow rate of refrigerant at inlet and outlet $\dot {m}_{ci}$ and $\dot {m}_{co}$; refrigerant enthalpy at inlet $h_{cri}$; air temperature $T_{ca}$ and air mass flow rate $\dot {m}_{ca}$. The $f_c \left( {x_c, u_c } \right)$ is expressed as the equation on the top of this page.

Expansion valve. The expansion valve is also modeled as a static component; its mass flow rate can be calculated from the orifice equation

 \begin{align} \dot {m}_v =C_v A_v \left[{\rho _v \left( {P_{vi}-P_{vo} } \right)} \right]^n \end{align} (5)

where $\dot {m}_v$ is the mass flow rate of the refrigerant through the expansion valve, $C_v$ is the orifice coefficient, $A_v$ is the opening area, $\rho _v$ is the refrigerant density, $P_{vi}$ and $P_{vo}$ are the inlet pressure and outlet pressure across the expansion valve, respectively.

Evaporator. Similar to the condenser model, the evaporator can be divided into two regions, i.e., a two-phase region with a mean void fraction, and a superheated region. According to the refrigerant mass conservation and energy balance, it can be formulated by a fifth order nonlinear model:

 \begin{align} &Z_e \left( {x_e, u_e } \right)\times \dot {x}_e =\notag\\ &\quad \left[{{\begin{array}{*{20}c} {\dot {m}_{ei} (h_{ei}-h_{eg} )+\alpha _{ei1} A_{ei} \left( {\frac{L_{e1} }{L_{eT} }} \right)\left( {T_{ew1}-T_{er1} } \right)} \hfill \\ {\dot {m}_{eo} (h_{eg}-h_{eo} )+\alpha _{ei2} A_{ei} \left( {\frac{L_{e2} }{L_{eT} }} \right)\left( {T_{ew2} -T_{er2} } \right)} \hfill \\ {\dot {m}_{ei} -\dot {m}_{eo} } \hfill \\ {\alpha _{eo} A_{eo} (T_{ea} -T_{ew1} )-\alpha _{ei1} A_{ei} (T_{ew1} -T_{er1} )} \hfill \\ {\alpha _{eo} A_{eo} (T_{ea} -T_{ew2} )-\alpha _{ei2} A_{ei} (T_{ew2} -T_{er2} )} \hfill \\ \end{array} }} \right] \end{align} (6)
 \begin{align} &\label{eq7} x_e =\left[{{\begin{array}{*{20}c} {L_{e1} } \hfill&{P_e } \hfill&{h_{eo} } \hfill&{T_{ew1} } \hfill & {T_{ew2} } \hfill \\ \end{array} }} \right]^{\rm T} \end{align} (7)
 \begin{align} &\label{eq8} u_e =\left[{{\begin{array}{*{20}c} {\dot {m}_{ei} } \hfill&{\dot {m}_{eo} } \hfill&{h_{eri} } \hfill & {T_{ea} } \hfill&{\dot {m}_{ea} } \hfill \\ \end{array} }} \right]^{\rm T} \end{align} (8)

where the state variables are: length of two phase flow $L_{e1}$; refrigerant pressure $P_e$; refrigerant enthalpy at outlet $h_{ero}$; the wall temperatures in the saturated and the superheated region, $T_{ew1}$ and $T_{ew2}$, respectively. The input variables are: mass flow rate of refrigerant at inlet and outlet, $\dot {m}_{ei}$ and $\dot {m}_{eo}$; refrigerant enthalpy at inlet $h_{eri}$; air temperature $T_{ea}$ and air mass flow rate $\dot {m}_{ea}$, respectively.

The combined cyclic model of integrated VCC system can be obtained by appropriately combining the component models according to the relations between the variables. The manipulated variables and output variables of the complete model are shown in (9) and (10). Different model structures can be obtained by different compositions of the manipulated variables and output variables. For the purposes of high accuracy and simple calculation, the optimized simplified model needs to be chosen for advancing towards further studies.

 \begin{align} &\quad u=\left[{{\begin{array}{*{20}c} {F_k } \hfill&{u_v } \hfill&{\dot {m}_{ca} } \hfill&{T_{cai} } \hfill {\dot {m}_{ea} } \hfill&{T_{eai} } \hfill \\ \end{array} }} \right]^{\rm T} \end{align} (9)
 \begin{align} &\begin{array}{l} y=[{\begin{array}{*{20}c} {L_{c1} } \hfill&{L_{c2} } \hfill&{P_c } \hfill&{h_{cro} } \hfill & {T_{wc1} } \hfill&{T_{wc2} } \hfill \\ \end{array} }\\ {\begin{array}{*{20}c} \quad {T_{wc3} } \hfill&{T_{cao} } \hfill&{T_{cro} } \hfill&{T_{crsc} } \hfill&{\dot {m}_{ca} } \hfill&{L_{e1} } \hfill&{P_e } &\hfill {h_{ero} }\hfill \\ \end{array} }\\ {\begin{array}{*{20}c} \quad \hfill {T_{crsh} } \hfill&{T_{we1} } \hfill&{T_{we2} } \hfill&{T_{eao} } \hfill&{T_{ero} } \hfill&{T_{esh} } \hfill & {\dot {m}_{ea}]^{\rm T}}. \end{array} } \end{array} \end{align} (10)
3 Structure selection criterion

A structure selection criterion, which is used to evaluate the performance of different reduced order models and choose the optimal simplified model, is proposed in the following section. The nonlinear VCC system can be represented by the discretized time-varying state-space model which can be formulated as

 \begin{align} \begin{array}{l} x\left( {k+1} \right)=A\left( k \right)x\left( k \right)+B\left( k \right)u\left( k \right) \\ y\left( k \right)=C\left( k \right)x\left( k \right) \\ \end{array} \end{align} (11)

where $u$ and $y$ are the inputs and outputs of the VCC system, shown as (9)-(10). The dimensions of them are $b\times 1$ and $c\times 1$, respectively. $A$ is the system matrix with dimension $a\times a$, $B$ is the control input matrix with dimension $a\times b$, $C$ is the output matrix with dimension $c\times a$. $k$ is the discrete sampling time.

To maintain the model accuracy while reducing the difficulty of controller design, appropriate model structure selection is an important part of process control. A popular model reduction method is the proper orthogonal decomposition (POD) method by which the order of model is reduced. The reduced order model can be obtained using snapshots. Since low order model is simple and easy to be controlled, the POD algorithm can to some extent simplify the complex models, but it does not show whether it is the optimized model for controller design.

According to the POD theory, the states and outputs of model can be represented by the reduced states and outputs:

 \begin{align} \begin{array}{l} x\left( k \right)=\Omega _1 x_n \left( k \right) \\ y\left( k \right)=\Omega _2 ^{-1}y_n \left( k \right) \\ \end{array} \end{align} (12)

where the $\Omega _1$ is formed by a linear combination of the snapshots with a dimension of $a\times n$, which are the dominant eigenvectors of the kernel matrix. $\Omega _2$ is the selection matrix which are composed of 0$'$s and 1$'$s, with a dimension of $n\times c$.

Substituting (12) into (11), the reduced order state-space model is formulated as

 \begin{align} \begin{array}{l} x_n \left( {k+1} \right)=A_n \left( k \right)x_n \left( k \right)+B_n \left( k \right)u\left( k \right) \\ y_n \left( k \right)=C_n \left( k \right)x_n \left( k \right) \\ \end{array} \end{align} (13)

where $A_n \left( k \right)=\Omega _1^{-1} A\left( k \right)\Omega _1$, $B_n \left( k \right)=\Omega _1^{-1} B\left( k \right)$, $C_n \left( k \right)=\Omega _2 C\left( k \right)\Omega _1$. The dimension of $A_n \left( k \right)$, $B_n \left( k \right)$, and $C_n \left( k \right)$ are $n\times n$, $n\times b$, and $n\times n$, respectively.

In this section, an analysis on the system variables of VCC system is carried out first to find which variables have a significant impact on the refrigeration cycle system. Then the important variables are selected for the following structure selection criterion.

Because of the differing units and disparate scaling of the variables, the inputs and outputs of the system are normalized first. The system is excited with random inputs, and an experimentally derived linear model is created under selected operating condition. Fig. 3 shows the frequency responses of four outputs (evaporator pressure, superheat, condenser pressure, and subcool) derived by step changes of the compressor speed. Clearly, compressor speed has a strong effect on all of the four outputs, while it has a stronger effect on evaporator pressure and superheat than condenser pressure and subcool especially at lower speeds.

 Download: larger image Fig. 3. Normalized frequency responses to step changes in compressor speed

Fig. 4 details the responses to changes in the expansion valve. Compared to the subcool, the effect is stronger on the superheat, while the effect is stronger on evaporator pressure than condenser pressure.

 Download: larger image Fig. 4. Normalized frequency responses to step changes in expansion valve

According to the analysis results, the evaporator pressure and the superheat of evaporator are chosen in the proposed criterion as the crucial influencing variables for controller design of this system. The overall system model is firstly reduced from 12th-order to reduced order models using POD method, and then giving random varying inputs to the full-order model and the same varying inputs to the reduced order models, the structure selection criterion is defined as

 \begin{align} J(i)\mbox{=}(\bar {P}_e (i)-P_e (i))^2+(\bar {T}_{esh} (i)-T_{esh} (i))^2, ~~i=1, \cdots, Z \end{align} (14)

where $\bar {P}_e$ and $\bar {T}_{esh}$ are the values of the evaporator pressure and the superheat of evaporator computed by the full-order model, $P_e$ and $T_{esh}$ are the values of the evaporator pressure and the superheat of evaporator computed by the reduced order model. $Z$ is the whole simulation time. This criterion evaluates the deviations between the full-order model and the reduced order model at each sample point $i$.

Denote the cost value $J$ in (14) of the reduced order model, the dimension of which is $n$, as $J_{n}$. The value of $J_{n}$ at time $i$ is obtained from the reduced order model correspondingly, and denoted as $J_n (i)$. Similar to $J_n (i)$, the value of $J_{n}$ at time $i+1$ is obtained from the reduced order model with dimension $n$, which is denoted as $J_n (i+1)$. Then we can define difference over the whole working time $Z$ as

 \begin{align} \label{eq12} \Delta J_{n} =\sum\limits_{i=1}^{Z-1} \left[-\frac{( {J_n (i+1)}J_n (i))^2}{J_n^2 (i)} \right]. \end{align} (15)

This criterion is to ensure a smoothly varying profile of the output variables. The structure selection criterion aims to find the optimal structure which minimizes the sum of quadratic partial variances of the evaporator pressure and the superheat of evaporator. With different dimension $n$, different $\Delta J_{n}$ can be computed. The model with dimension $n$ which minimizes $\Delta J_{n}$ should be chosen as the optimized model structure.

Thus, the structure selection criterion algorithm performs the following steps:

Step 1. Reduce the overall system model first from 12th-order model to a particular low order model using POD method, and then further reduce from that reduced order model to 2nd-order model.

Step 2. Give same various random inputs to the full-order model and the reduced order models, respectively.

Step 3. Compute the criterion $J$ between the reduced order models and the full-order model as described in (14) at each sampling time.

Step 4. Calculate the differences $\Delta J_{n}$ of the reduced order models expressed as (15) over the whole simulation time.

Step 5. Compare the values of $\Delta J_{n}$ and choose the model which minimums $\Delta J_{n}$ as the optimized model.

4 MPC controller design

MPC is a control algorithm which computes a sequence of control inputs based on an explicit prediction of outputs within some future horizon[22]. One of the most important advantages of MPC is that it accounts for the constraints of input and output variables that can be inherent to the real industrial systems, e.g., a valve cannot open past 100% open or close past 0% open. Another advantage of MPC is that additional constraints can be defined by the user to keep the system operating in a safe range, e.g., keeping evaporator superheat above a desired minimum.

The performance of MPC depends on a number of design parameters, such as length of the control time interval, the number of future moves for the manipulated variable, and the number of time intervals in the output prediction. In order to define how well the predicted process tracks the set points, an objective function $J_{MPC} (k)$ for the predictive control needs to be optimized as

 \begin{align} \label{eq13} \begin{array}{l} J_{MPC} (k)= \\ \sum\limits_{i=1}^P {\left\{ {\left( {w(k+i)-y(k+i)} \right)^{\rm T}Q\left( {w(k+i)-y(k+i)} \right)} \right\}}+ \\ \, \sum\limits_{i=1}^M {\left\{ {\left( {\Delta u(k+i-1)} \right)^{\rm T}{\pmb R}\left( {\Delta u(k+i-1)} \right)} \right\}} \\ \end{array} \end{align} (16)

where $k$ is the current sampling interval, $k+i$ is the future sampling interval, $P$ is the prediction horizon, $Q$ is the weight matrix of outputs, $w(k+i)$ is the desired output at instant $k+i$, $y(k+i)$ is the actual output at instant $k+i$, $M$ is the control horizon, ${\pmb R}$ is the weight matrix of inputs and $\Delta u(k+i-1)$ is the predicted adjustment of input $u$ at future instant $k+i-1$.

Equation (16) computes the weighted sum of squared deviations of the outputs from the set points and the weighted sum of squared incremental manipulated variables.

The optimization algorithm searches for values of $\Delta u$ over the control horizon that minimize the objective function. Although $M$ control moves are calculated, only the first move is implemented. At the next sampling time, the weights are updated and the optimization is done again. Input and output values have a limited range to operate, therefore $\Delta u$, $u$ and $y$ are subjected to constraints introduced as (17)-(19).

 $\Delta u_{\min } \le \Delta u(k)\le \Delta u_{\max }$ (17)
 $u_{\min } \le u(k)\le u_{\max }$ (18)
 $y_{\min } \le y(k)\le y_{\max }.$ (19)

The objective function $J_{MPC} (k)$ in (16) is for the integrated model of VCC system, which is not suitable for the reduced order models. Therefore, the objective function for the reduced order model with dimension $n$ is formulated as follows:

 \begin{align} &J_{MPC, n} (k)= \notag\\ &\quad \sum\limits_{i=1}^P \left\{ \left( {w(k+i)-y_n (k+i)}\right)^{\rm T}\Omega _2^{-\rm T}\right.\notag\\ &\quad \left.Q\Omega _2^{-1} \left( {w(k+i)-y_n (k+i)} \right) \right\} +\notag\\ &\quad \sum\limits_{i=1}^M {\left\{ {\left( {\Delta u_n (k+i-1)} \right)^{\rm T}{\pmb R}\left( {\Delta u_n (k+i-1)} \right)} \right\}} \end{align} (20)
 \begin{align} &\Delta u_{\min, n} \le \Delta u_n (k)\le \Delta u_{\max, n} \notag\\ &u_{\min, n} \le u_n (k)\le u_{\max, n} \notag\\ &y_{\min, n} \le y_n (k)\le y_{\max, n} \end{align} (21)

where $J_{MPC, n} (k)$ is the objective function for the reduced order model with dimension $n$. $w(k+i)$ is the desired output for the reduced order model, while $y_n (k+i)$ is the actual output for the reduced order model.

The MPC control algorithm for VCC system is as follows:

Step 1. Compute the objective function $J_{MPC} (k)$ or $J_{MPC, n} (k)$ as (16) or (20) by iteratively modifying $\Delta u\left( k \right)$.

Step 2. Calculate the best control signal $u(k)$ or $u_n (k)$.

Step 3. Apply $u(k)$ or $u_n (k)$ to the VCC models and calculate the actual value of outputs.

Step 4. At the next sampling time, go back to calculate $\Delta u$.

5 Experiment results

The experimental platform used in this research is developed at the process instrumentation laboratory of Nanyang Technological University of Singapore. The photograph and schematic diagram of the experimental platform are shown in Figs. 5 and 6.

 Download: larger image Fig. 5. Photograph of the experimental platform in this research

 Download: larger image Fig. 6. Schematic of the experimental platform in this research

The system includes a variable speed compressor, an electronic expansion valve, an air-cooled condenser and evaporator, a liquid receiver after the condenser, an accumulator after the evaporator, and the fans of condenser and evaporator with variable frequency. In addition, the pressure measurement devices are installed on the system, and the measurement range and measuring error of these devices are 0$\sim$1600 kPa and $\pm$0.5%. The temperature sensors are also installed with the measurement range of $-$40${^\circ}\sim$ 200${^\circ}$ and measuring error of $\pm$0.3${^\circ}$. The locations of their installation are shown in Fig. 6. The R134a is selected as the refrigerant in this research. The measuring error of the refrigerant mass flow in this system is $\pm$1.6%.

5.1 Simulation results of the model structure selection

To verify that the structure selection criterion can choose optimized model structure effectively, a comparison study has been carried out. The model has been firstly linearized around a steady state operating point indicated in Table 1. Then with the POD method, the system model is reduced to 4th-order model with the selection of input variables as $u_4 =\left[ {{\begin{array}{*{20}c} {F_k } \hfill&{u_v } \hfill&{\dot {m}_{ca} } \hfill&{\dot {m}_{ea} } \hfill \\ \end{array} }} \right]^{\rm T}$and output variables as $y_4 =\left[ {{\begin{array}{*{20}c} {P_c } \hfill&{T_{crsc} } \hfill&{P_e } \hfill&{T_{esh} } \hfill \\ \end{array} }} \right]^{\rm T}$. While the 4th-order model is still complicated to control, according to the Step 1 of structure selection criterion, it is further reduced to 3rd-order and 2nd-order model with the selection of input variables as $u_3 =\left[ {{\begin{array}{*{20}c} {F_k } \hfill&{u_v } \hfill&{\dot {m}_{ca} } \hfill \\ \end{array} }} \right]^{\rm T}$, $u_2 =\left[{{\begin{array}{*{20}c} {F_k } \hfill&{u_v } \hfill \\ \end{array} }} \right]^{\rm T}$ and output variables as $y_3 =\left[ {{\begin{array}{*{20}c} {P_c } \hfill&{P_e } \hfill&{T_{esh} } \hfill \\ \end{array} }} \right]^{\rm T}$, $y_2 =\left[{{\begin{array}{*{20}c} {P_e } \hfill&{T_{esh} } \hfill \\ \end{array} }} \right]^{\rm T}$, respectively.

Table 1
The steady state operating conditions for the model linearization

The $\Delta J_{n}$ among the 4th-order, 3rd-order, and 2nd-order models are compared in this section. The same square wave variations in compressor speed within a certain range from 1100 rpm to 1150 rpm are firstly given to the full-order model and three reduced order models. Then the system outputs profiles are shown in Figs. 7 and 8, which compare the outputs performance among the four model structures under the same variations of compressor speed.

 Download: larger image Fig. 7. Pressure of evaporator for random variations in compressor speed

 Download: larger image Fig. 8. Superheat of the evaporator for random variations in compressor speed

Through the simulation, the structure selection criteria of the three simplified structures are calculated according to (14)-(15) respectively, and the corresponding structure selection criteria are shown in (22)-(24).

 \begin{align} &J_2 (i)\mbox{=}(\bar {P}_e (i)-P_{e2} (i))^2+(\bar {T}_{esh} (i)-T_{esh2} (i))^2 \notag\\ &\Delta J_2 =\sum\limits_{i=1}^{Z-1} {\left[{\frac{\left( {J_2 (i+1)-J_2 (i)} \right)^2}{J_2^2 (i)}} \right]}~~i=1, \cdots, Z \end{align} (22)
 \begin{align} &J_3 (i)\mbox{=}(\bar {P}_e (i)-P_{e3} (i))^2+(\bar {T}_{esh} (i)-T_{esh3} (i))^2\notag \\&\Delta J_3 =\sum\limits_{i=1}^{Z-1} {\left[{\frac{\left( {J_3 (i+1)-J_3 (i)} \right)^2}{J_3^2 (i)}} \right]} ~~ i=1, \cdots, Z \end{align} (23)
 \begin{align} &J_4 (i)\mbox{=}(\bar {P}_e (i)-P_{e4} (i))^2+(\bar {T}_{esh} (i)-T_{esh4} (i))^2 \notag\\ &\Delta J_4 =\sum\limits_{i=1}^{Z-1} {\left[{\frac{\left( {J_4 (i+1)-J_4 (i)} \right)^2}{J_4^2 (i)}} \right]}~~i=1, \cdots, Z \end{align} (24)

where $P_{e2}$, $P_{e3}$, $P_{e4}$, $T_{esh2}$, $T_{esh3}$, $T_{esh4}$, $J_2$, $J_3$, $J_4$, $\Delta J_2$, $\Delta J_3$, $\Delta J_4$ are the evaporator pressure, superheat, criterion and criterion difference computed with the 2nd-order model, 3th-order model, and 4th-order model, respectively. $\bar {P}_e$ and $\bar {T}_{esh}$ are evaporator pressure and superheat of evaporator computed by the full-order model. $Z$ is simulation time and $i$ is the sample instant.

The results of the proposed structure selection criterions are shown in Table 2 which indicate that the 3rd-order model structure has the minimum structure selection criterion. Thus the 3rd-order model structure is finally chosen as optimal model structure, which is convenient to the controller design of the VCC system. The coefficient matrices of the 3rd-order model are expressed as (25).

Table 2
The structure selection criterion

 \begin{align} &A_3 \mbox{=}\left[{{\begin{array}{*{20}c} {\mbox{-$0.152 5}} \hfill&{\mbox{$-$0.010 6}} \hfill&{\mbox{0.046 6}} \hfill \\ {\mbox{$-$0.000 0}} \hfill&{\mbox{$-$0.171 8}} \hfill&{\mbox{0.244 1}} \hfill \\ {\mbox{$-$0.000 0}} \hfill&{\mbox{$-$0.128 1}} \hfill&{\mbox{$-$0.010 3}} \hfill \\ \end{array} }} \right]\notag\\ &B_3 \mbox{=}\left[{{\begin{array}{*{20}c} {\mbox{$-$0.007 2}} \hfill&{\mbox{0.131 5}} \hfill&{\mbox{2.233 9}} \hfill \\ {\mbox{$-$0.011 5}} \hfill&{\mbox{3.810 7}} \hfill&{\mbox{$-$4.729 4}} \hfill \\ {\mbox{0.018 9}} \hfill&{\mbox{$-$0.507 7}} \hfill&{\mbox{$-$1.412 1}} \hfill \\ \end{array} }} \right]\notag\\ &C_3 \mbox{=}\left[{{\begin{array}{*{20}c} {\mbox{3.416 6}} \hfill&{\mbox{$-$1.365 6}} \hfill&{\mbox{$-$1.998 3}} \hfill \\ {\mbox{0.937 6}} \hfill&{\mbox{0.173 5}} \hfill&{\mbox{$-$1.145 9}} \hfill \\ {\mbox{$-0.069 0}} \hfill&{\mbox{0.273 5}} \hfill&{\mbox{0.522 1}} \hfill \\ \end{array} }} \right]\notag\\ &D_3 \mbox{=}\left[{{\begin{array}{*{20}c} \mbox{0} \hfill&\mbox{0} \hfill&\mbox{0} \hfill \\ \mbox{0} \hfill&\mbox{0} \hfill&\mbox{0} \hfill \\ \mbox{0} \hfill&\mbox{0} \hfill&\mbox{0} \hfill \\ \end{array} }} \right]. \end{align} (25)
5.2 Experimental results of MPC controller

In this section, the MPC controller for the proposed 3rd-order model is designed based on objective function $J_{MPC, 3} (k)$ computing by (20). The condenser pressure, evaporator pressure, and superheat are the regulated outputs, and the compressor speed, expansion valve opening, and the air mass flow rate through the condenser are controllable inputs of the 3rd-order model. In the MPC controller, weights of 1, 1 and 100 are placed on the condenser pressure, evaporator pressure and superheat, respectively. Rate weights of 0.1, 0.01 and 0.001 are placed on the compressor speed, expansion valve opening and the air mass flow rate through the condenser, respectively. A control interval of 1 s, a control horizon of 15 intervals, and a prediction horizon of 50 intervals are determined as the tuning parameters of the MPC controller.

To demonstrate the effectiveness of the designed control system, an MPC controller for the 2nd-order model is designed as follows: the evaporator pressure and superheat are the system outputs, meanwhile the compressor speed, expansion valve opening are inputs. The weight of 1 and 100 are placed on the evaporator pressure and superheat, and rate weights of 0.1 and 0.01 on the compressor speed, expansion valve in MPC controller, respectively. Other tuning parameters are the same with controller for 3rd-order model. The performance comparison between the two controllers are carried out, and the output and input profiles are shown in Figs. 9 and 10.

 Download: larger image Fig. 9. Outputs comparison of VCC system under the two MPC controllers

 Download: larger image Fig. 10. Inputs comparison of VCC system under the two MPC controllers

Fig. 9 shows the output performance in response to the disturbances due to the changes in the air mass flow rate of evaporator at 1000s. The 3rd-order controller tracks the setting value faster than the 2nd-order controller when the disturbance occurs, which demonstrate that the proposed control system has satisfactory tracking performance and robustness against disturbance.

The corresponding input variables are given as Fig. 10. When the disturbance is changed, the compressor speed and the expansion valve opening respond to change quickly in order to reject the disturbance and keep the superheat and pressure as before.

6 Conclusions

This paper presented a model predictive control system based on a structure selection criterion which is used to simplify model and select the optimal reduced order model for MPC controller. The overall system model is reduced from 12th-order to reduced order models using POD method, and the model which minimizes the cost value is chosen as the optimized simplified model. Based on the proposed model, the MPC based controller is designed, and the objective function for the simplified model is optimized. To validate the effectiveness of the proposed control system, comparison experiments have been carried out, and the experimental results indicate that the proposed controller has excellent advantages of tracking performance and disturbance rejection.

Acknowledgements

This work was supported by National Natural Science Foundation of China (Nos. 61233004, 61221003, 61374109 and 61473184), National Basic Research Program of China (973 Program)(No. 2013CB035500), and partly sponsored by the Higher Education Research Fund for the Doctoral Program of China (No. 20120073130006), and National Research Foundation of Singapore (No. NRF2011 NRF-CRP001-090).

References
 [1] S. Skogestad, I. Postlethwaite. Multivariable Feedback Control:Analysis and Design, New York, USA: Wiley, 1996. [2] A. C. Antoulas. An overview of approximation methods for large-scale dynamical systems. Annual Reviews in Control, vol.29, no.2, pp.181-190, 2005. DOI:10.1016/j.arcontrol.2005.08.002 [3] Y. E. An, C. Q. Gu. Model reduction for large-scale dynamical systems via equality constrained least squares. Journal of Computational and Applied Mathematics, vol.234, no.8, pp.2420-2431, 2010. DOI:10.1016/j.cam.2010.03.006 [4] B. C. Moore. Principal component analysis in linear systems:Controllability, observability, and model reduction. IEEE Transactions on Automatic Control, vol.26, no.1, pp.17-32, 1981. DOI:10.1109/TAC.1981.1102568 [5] K. V. Fernando, H. Nicholson. Singular perturbational model reduction of balanced systems. IEEE Transactions on Automatic Control, vol.27, no.2, pp.466-468, 1982. DOI:10.1109/TAC.1982.1102932 [6] B. P. Rasmussen, A. G. Alleyne. Dynamic Modeling and Advanced Control of Air Conditioning and Refrigeration Systems, Ph.dissertation D., USA: University of Illinois, 2006. [7] X. L. Wang, Y. L. Jiang. Model order reduction methods for coupled systems in the time domain using Laguerre polynomials. Computers & Mathematics with Applications, vol.62, no.8, pp.3241-3250, 2011. [8] P. Guha, S. Mishra. Reduced order modeling & controller design for mass transfer in a grain storage system. International Journal of Automation and Computing, vol.11, no.4, pp.399-403, 2014. DOI:10.1007/s11633-014-0805-6 [9] K. M. B. Legweel, D. V. Lazić, M. R. Ristanović, J. V. Lozanović Šajić. The performance of pip-cascade controler in HVAC system. Thermal Science, vol.18, no.Suppl.1, pp.S213-S220, 2014. DOI:10.2298/TSCI130812183L [10] M. Xu, S. Y. Li, W. J. Cai, L. Lu. Effects of a GPC-PID control strategy with hierarchical structure for a cooling coil unit. Energy Conversion and Management, vol.47, no.1, pp.132-145, 2006. DOI:10.1016/j.enconman.2005.03.012 [11] Y. L. Shen, W. J. Cai, S. Y. Li. Normalized decoupling control for high-dimensional MIMO processes for application in room temperature control HVAC systems. Control Engineering Practice, vol.18, no.6, pp.652-664, 2010. DOI:10.1016/j.conengprac.2010.03.006 [12] J. N. Lee, T. M. Lin, C. C. Chen. Modeling validation and control analysis for controlled temperature and humidity of air conditioning system. The Scientific World Journal, vol.2014, pp.903032, 2014. [13] M. S. Elliott, B. P. Rasmussen. On reducing evaporator superheat nonlinearity with control architecture. International Journal of Refrigeration, vol.33, no.3, pp.607-614, 2010. DOI:10.1016/j.ijrefrig.2009.12.013 [14] M. Mohanraj, S. Jayaraj, C. Muraleedharan. Applications of artificial neural networks for refrigeration, air-conditioning and heat pump systems-A review. Renewable and Sustainable Energy Reviews, vol.16, no.2, pp.1340-1358, 2012. DOI:10.1016/j.rser.2011.10.015 [15] A. Preglej, J. Rehrl, D. Schwingshackl, I. Steiner, M. Horn, I. Škrjanc. Energy-efficient fuzzy model-based multivariable predictive control of a HVAC system. Energy and Buildings, vol.82, pp.520-533, 2014. DOI:10.1016/j.enbuild.2014.07.042 [16] S. Hussain, H. A. Gabbar, D. Bondarenko, F. Musharavati, S. Pokharel. Comfort-based fuzzy control optimization for energy conservation in HVAC systems. Control Engineering Practice, vol.32, pp.172-182, 2014. DOI:10.1016/j.conengprac.2014.08.007 [17] B. P. Rasmussen, A. G. Alleyne. Gain scheduled control of an air conditioning system using the Youla parameterization. IEEE Transactions on Control Systems Technology, vol.18, no.5, pp.1216-1225, 2010. DOI:10.1109/TCST.2009.2035104 [18] M. Wallace, B. Das, P. Mhaskar, J. House, T. Salsbury. Offset-free model predictive control of a vapor compression cycle. Journal of Process Control, vol.22, no.7, pp.1374-1386, 2012. DOI:10.1016/j.jprocont.2012.06.011 [19] M. Wallace, R. McBride, S. Aumi, P. Mhaskar, J. House, T. Salsbury. Energy efficient model predictive building temperature control. Chemical Engineering Science, vol.69, no.1, pp.45-58, 2012. DOI:10.1016/j.ces.2011.07.023 [20] M. S. Elliott, B. P. Rasmussen. Model-based predictive control of a multi-evaporator vapor compression cooling cycle. In Proceedings of the American Control Conference, IEEE, Seattle, USA, pp.1463-1468, 2008. [21] X. H. Xu, S. W. Wang, G. S. Huang. Robust MPC for temperature control of air-conditioning systems concerning on constraints and multitype uncertainties. Building Services Engineering Research & Technology, vol.31, no.1, pp.39-55, 2010. [22] M. L. Wang, N. Li, S. Y. Li. Model-based predictive control for spatially-distributed systems using dimensional reduction models. International Journal of Automation and Computing, vol.8, no.1, pp.1-7, 2011. DOI:10.1007/s11633-010-0547-z