Optimal and Model Predictive Control of Single Phase Natural Circulation in a Rectangular Closed Loop

1. Introduction

Pipelines are basically a network of interconnected pipes linked with pumps, control valves, instruments etc. that are used to transport the fluids from one point of location to another location  [1]. The pipeline systems are integral to various industries as they enable the transportation of fluids over long distances. There are various applications of these pipelines as these are abundantly used in chemical industries, geothermal systems, solar water heating systems, oil and gas sector, water distribution, wastewater management, power generation and many more.

Various geometries and conditions were employed in such investigations, each possessing its own significance. Literature pertaining to these issues can be categorized based on the geometries and investigation methods utilized. Regarding geometry classification, it may involve either an open-loop or closed-loop configuration. Different methods for conducting these studies can be characterized into various numerical techniques, diverse analytic approaches, and various experimental methods.

This particular section encompasses the literature review on these three items, (a) Effect of using various configurations for natural circulation,(b) Effect of varying boundary and operating conditions,(c) Process dynamics and control of thermo-syphoning loops.

Natural circulation or thermo-syphoning loops are broadly categorized into two main types: single-phase thermo-syphoning loops and two-phase thermo-syphoning loops. This literature review primarily focuses on single-phase natural circulation in thermo-syphoning loops, where heat is transferred from a heat source to a heat sink. The orientation of the thermo-syphoning loop significantly influences its performance  [2]. Schematic representations of some of these loops are provided in Figure 1  [3] .

In Figure 1 (a), a typical vertical loop is depicted, comprising a heat source and a heat sink. Heat moves from the source to the sink due to density or buoyancy gradients. Figure 1 (b) illustrates toroidal loops, where constant heat flux is typically applied to the bottom half while the upper half maintains a constant wall temperature boundary condition. Once again, heat transfer occurs through natural circulation driven by density or buoyancy gradients. Figure 1(c) demonstrates parallel branches with various heat or cooling fluxes of different magnitudes, facilitating the natural circulation phenomenon. Figure 1(d) displays an open square loop with heating coils at the bottom and two cooling coils open-ended at the top. Cold water flows from cold storage tanks under gravity, passing through heating coils, and warm water flows in the same legs towards cold storage tanks. In Figure 1(e), a closed square loop is shown, featuring a heating section at the bottom and a cooling section at the top.

However, the performance of rectangular closed loop was compared with toroidal loop by Basu et al.  [4] and it was concluded that rectangular loops induce larger fluid flows than the toroidal loops undergoing identical dimensions, operating and boundary conditions, however, the stability span for toroidal geometry was found wider as compared to rectangular geometry.

Fichera et al.  [5] compared the dynamics of two rectangular loops having different tube diameters and lengths, reported the important results that the loop having smaller diameter demonstrate more periodic behavior as compared to the loop having larger dimensions whose behavior is relatively more chaotic.

The effect of varying thermal boundary condition was studied by misale et al.  [6], and it was inferred that varying the heat sink temperature could lead to crossing the stability threshold, highlighting the significance of this thermal boundary condition in influencing the system’s dynamics. In another study  [7], it was noticed that as heater power increases, there’s a notable shift in the dynamics of fluid flow within the loop. Initially, loop characteristics remain stable up to a specific heater power level. However, beyond this threshold, the system transitions from steady behavior to periodic oscillations, and eventually to chaotic patterns.

The dynamics of rectangular closed loop was compared by subjecting it to different types of input powers such as step, ramp, and exponential profiles  [8]. Exponential signals in modified form outperforms others in terms of stability both during increasing and decreasing power levels. Hashemi-Tilehnoee et al.  [9] highlighted the role of natural circulation loops in passive heat transfer across diverse applications, emphasizing advances in working fluids, configurations, and modeling techniques to improve efficiency, stability, and safety in next-generation thermal systems.

A neural network–based NARMAX model  [10] was developed to capture the unstable dynamics of a closed-loop thermosyphon from experimental input–output data, showing much closer agreement with measurements than an existing mathematical model. However, experimental and theoretical studies on thermal convection loops have shown that optimal and adaptive controllers can effectively suppress chaotic convection, outperforming neural network controllers in stability and efficiency, though the latter offer advantages in handling unknown nominal conditions and time delays [11].

Muscato et al.  [12] derives a low-order non-linear model aligned with experimental observations, validated through simulation-experiment comparisons. It was observed that the dynamics of the natural circulation becomes unstable on higher heat fluxes, however at lower heat fluxes, the flow in the rectangular loop is not so efficient in terms of heat exchange in the cooling section. This model facilitated controller design for fluid motion stabilization, ensuring effective heat dissipation. Both conventional proportional and proportional-derivative controllers, considering the fluid velocity or temperature difference as feedback, were designed for various heating powers. Simulation results demonstrated confirmed the controllers’ effectiveness in stabilizing the system at desired equilibrium points.

A. Fichera et al.  [13] developed a finite-order model of rectangular natural circulation loops by reducing a Fourier-series-based formulation of the momentum and energy equations to seven equations capturing the first three dynamic modes. Compared with Muscato et al.’s lower-order lumped-parameter model, which targeted specific nominal equilibrium points, Fichera’s higher-fidelity approach reproduced complex oscillatory behavior more accurately and supported controller design for a broader range of unstable dynamics.

Cammarata et al.  [14] implemented an optimal control in the rectangular configuration and reported that the adopted strategy for the optimal controller is quite capable to minimizes a functional considering both system stability and control cost, weighted appropriately. Simulations and experimental implementation confirmed the effectiveness of the optimal controller in suppressing the hazardous oscillations and flow reversals, guiding the system to the desired steady state.

In geothermal pipeline or solar water heating systems based on Natural Circulation Loops (NCLs), where the primary objective is to extract maximum heat with a cooling water, stable and efficient operation is fundamental to ensuring long-term sustainability of energy extraction. Due to their reliance on buoyancy-driven flow without mechanical pumping, NCLs are inherently sensitive to variations in thermal and hydraulic conditions, making them vulnerable to instabilities, oscillations, and thermal stresses. While advanced control methods have been investigated for various renewable energy systems, the application of Model Predictive Control (MPC) to NCL-driven geothermal pipelines has not, to the best of our knowledge, been reported in the literature. This study addresses this research gap by introducing a predictive, optimization-based control strategy that utilizes an NCL dynamic model to forecast system responses and implement optimal control actions. The proposed approach enhances thermal stability, mitigates flow oscillations to prevent the excessive material stresses, thereby supporting both operational efficiency and the sustainable exploitation of geothermal or solar resources.

A key strength of the proposed MPC framework lies in its ability to explicitly incorporate operational and physical constraints into the control optimization problem. These constraints include maximum allowable temperatures and pressures to preserve material integrity, flow rate limits to minimize erosion and scaling, and bounds on heat extraction rates to prevent reservoir over-exploitation. Unlike conventional control methods, MPC ensures that these safety and sustainability criteria are respected at all times, even under variable operating conditions. By integrating predictive control with constraint handling, this work offers a novel pathway to improve process reliability to extend infrastructure lifespan, and enhance the environmental compatibility of geothermal energy systems. In this particular study, the test assembly consists of a rectangular loop with a heating section in the bottom horizontal leg and a cooling section in the top horizontal leg, connected by adiabatic vertical riser and downcomer tubes. The primary objective is to develop an MPC framework capable of achieving effective and stable heat removal by regulating the cooling rate, while explicitly incorporating physical constraints such as maximum and minimum coolant flow capacity. By leveraging MPC’s predictive capabilities, the proposed approach aims to enhance system stability, improve operational safety, and optimize equipment sizing, ultimately contributing to energy efficiency, reduced maintenance requirements, and the broader adoption of sustainable passive thermal management technologies in industrial applications.

2. Configuration of Experimental Setup

The Figure 2 illustrates the experimental arrangement for a single-phase natural circulation system  [12]. This setup primarily comprises heating and cooling sections positioned in the lower and upper horizontal parts of the loop, respectively, interconnected via adiabatic vertical glass tubes. The heating section utilizes a variable heater with adjustable heating power, while the cooling section incorporates a double-pipe heat exchanger. Here, the main water flow passes through the inner copper tube, while tap water is injected into the annular pipe to extract heat from the inner fluid. The flow rate of the tap water is regulated by a controller.

The dimensions of the entire experimental setup are outlined in Table 1.

3. Mathematical Modeling for Natural Circulation Loop (NCL)

The assumptions that have been made for the natural circulation model are;(i) The bulk fluid flow is only in axial direction, (ii) The system behaves iso-thermally in the entire loop except heating and cooling sections. Therefore, the heat supplied by the heating section will be extracted in the cooling section, (iii) The cross-sectional area of loop remains constant throughout the axial direction, (iv) The fluid in the loop remains single phase during the entire operation. The mathematical model for natural circulation based on the fundamental conservation laws of momentum and energy can be described  [12] as follows.

$$T_t(\zeta ,t)+v(t)T_{\zeta }(\zeta ,t)=h(\zeta )+\alpha T_{\zeta \zeta }(\zeta ,t)$$

(1)

$$v_t(t)+{\displaystyle \frac{1}{D}}{\left({\displaystyle \frac{v_c}{D}}\right)}^d{v}^{2-d}={\displaystyle \frac{g\beta }{2(L+H)}}\int \mathsf{\Phi }(T-T_0)f(\zeta )d\zeta$$

(2)

$$T(\zeta ,0)=T_0(\zeta )$$

(3)

$$v(0)=v_0$$

(4)

Where $\zeta \in [0,2(H+L)]$ and $t\in [0,\infty )$ denote the spatial and temporal domains respectively. Here, T, v, $\alpha$, D, g, $\beta$, $T_0$, and $v_0$ represent the fluid temperature in the loop, fluid velocity, thermal diffusivity, inner tube diameter, gravitational acceleration, volumetric expansion coefficient, reference temperature, and initial velocity, respectively.

The functions $f(\zeta )={\displaystyle \frac{dz}{d\zeta }}$ and $h(\zeta )$ depend on the geometry of the rectangular closed-loop system and are given as:

$$f(\zeta )=\left\{\begin{array}{cc}1,\hfill & 0<\zeta <L\hfill \\ 0,\hfill & L<\zeta <L+H\hfill \\ -1,\hfill & L+H<\zeta <2L+H\hfill \\ 0,\hfill & 2L+H<\zeta <2(L+H)\hfill \end{array}\right.h(\zeta )=\left\{\begin{array}{cc}0,\hfill & 0<\zeta <L\hfill \\ -Q_L\hfill \\ 0,\hfill & L+H<\zeta <2L+H\hfill \\ +Q_H,\hfill & 2L+H<\zeta <2(L+H)\hfill \end{array}\right.$$

(5)

where $Q_L$ and $Q_H$ and in Equation (5) can be defined as follows;

$$Q_L={\displaystyle \frac{4}{{\rho }_o{C}_{p}D}}.{\displaystyle \frac{\dot{m}{c}_p\mathsf{\Delta }{T}_c}{\pi DH}};Q_H={\displaystyle \frac{4}{{\rho }_o{C}_{p}D}}.q$$

(6)

Where q is the heat flux, $\mathsf{\Delta }{T}_c$ is the temperature gradient between the entrance and exit of the cooling section, m is the mass flow rate, and $C_p$ is the heat capacity of the cooling water.

A key physical assumption that ensures the presence of equilibrium points is that the total heat energy transferred to the fluid as it passes through the heat source is exactly balanced by the heat energy removed as it flows through the heat sink. Based on this principle, we have:

$$\oint h(\zeta )=0$$

(7)

which in turn implies:

$$-{\displaystyle \frac{\dot{m}{c}_p\mathsf{\Delta }{T}_c}{\pi DH}}=q$$

(8)

The main model equations (1) and (2) can be further simplified by introducing Fourier Series expansion [15] of state variables T as follows:

$$T(\zeta ,t)-T_0=\sum _{k\in K}{a}_k(t)e^{i{\displaystyle \frac{2\pi }{2(L+H)}}k\zeta }$$

(9)

$$h(\zeta )=\sum _{k\in K}{b}_k{e}^{i{\displaystyle \frac{2\pi }{2(L+H)}}k\zeta }$$

(10)

$$f(\zeta )=\sum _{k\in K}{c}_k{e}^{i{\displaystyle \frac{2\pi }{2(L+H)}}k\zeta }$$

(11)

Here, $a_k(t)={\alpha }_k(t)+j{\beta }_k(t)$ is a complex number. From equation (5) and geometrical considerations, where $K,J\subset \mathbb{Z}\setminus \left\{0\right\}$ (i.e., the set of all nonzero integers), substituting expressions (9)–(11) into (1) and (2) yields the following system of infinite ordinary differential equations:

$${\displaystyle \frac{\mathrm{d}v(t)}{\mathrm{d}t}}+{\displaystyle \frac{1}{r}}b{\left({\displaystyle \frac{v_c}{2r}}\right)}^d{v}^{2-d}=g\beta \sum _{k\in K\cap J}{a}_k(t)·{\overline{c}}_k$$

(12)

$${\displaystyle \frac{\mathrm{d}{a}_k(t)}{\mathrm{d}t}}+\left[{\displaystyle \frac{\pi }{(L+H)}}ikv(t)+a{\displaystyle \frac{{\pi }^2}{{(L+H)}^2}}{k}^2\right]·a_k(t)=b_k,k\subset K\cap J$$

(13)

where ${\overline{a}}_k(t)=a_{-k}(t)$, ${\overline{b}}_k=b_{-k}$, and ${\overline{c}}_k=c_{-k}$. Note that modes not in K decay, while modes in K but not in J are dependent variables.

In the case of the rectangular loops considered here, the coefficients $c_k$ from the expansion of $f(\zeta )$ are obtained from the following expression:

$$c_k={\displaystyle \frac{1}{2(L+H)}}\oint f(\zeta )e^{-i{\displaystyle \frac{2\pi }{2(L+H)}}k\zeta }d\zeta$$

(14)

For simplification, we introduce the following substitutions:

$$\eta ={\displaystyle \frac{\pi \zeta }{L+H}},\gamma ={\displaystyle \frac{\pi L}{L+H}}$$

(15)

With the substitutions applied and considering the specific form of $f(\zeta )$, the coefficients $c_k$ for the rectangular geometry simplify to:

$$c_k={\displaystyle \frac{1}{2\pi ki}}\left[(1-cos(k\gamma )+isin(k\gamma ))(1-cos(k\pi ))\right]$$

(16)

The coefficients $c_k$ vanish when k is even. As a result, the set J consists only of odd integers. Furthermore, some coefficients $c_k$ may also vanish for odd k if the circuit’s aspect ratio is rational.

In a similar fashion, the coefficients $b_k$ arising from the expansion of $h(\zeta )$ can be computed using:

$$b_k={\displaystyle \frac{1}{2(L+H)}}\oint h(\zeta )·e^{-i{\displaystyle \frac{2\pi }{2(L+H)}}k\zeta }\mathrm{d}\zeta$$

(17)

By again applying substitution (15) and introducing the following definitions:

$$\mathsf{\Gamma }={\displaystyle \frac{4}{{\rho }_{0}cD}}\left({\displaystyle \frac{\dot{m}{c}_p\mathsf{\Delta }{T}_c}{\pi DH}}+q\right),{\mathsf{\Gamma }}_1={\displaystyle \frac{4}{{\rho }_{0}cD}}\left({\displaystyle \frac{\dot{m}{c}_p\mathsf{\Delta }{T}_c}{\pi DH}}-q\right)$$

(18)

the expressions for $b_k$ in the considered case are obtained as:

$$b_k={\displaystyle \frac{\mathsf{\Gamma }}{2\pi ki}}(-1-cosk\gamma +isink\gamma )\mathrm{for}\mathrm{odd}k=2n+1,(n=\pm 1,\pm 2,\dots )$$

(19)

$$b_k={\displaystyle \frac{{\mathsf{\Gamma }}_1}{2\pi ki}}(1-cosk\gamma +isink\gamma )\mathrm{for}\mathrm{even}k=2n,(n=\pm 1,\pm 2,\dots )$$

(20)

Notably, since from the earlier result ${\mathsf{\Gamma }}_1=0$, this leads to the conclusion that all even-mode coefficients $b_k$ vanish. Therefore, the set K again consists solely of odd integers.

Now, after plugging eqs. (10)–(12) into the main model equations (1) and (2), and subsequently calculating the respective coefficients, truncating the Fourier Series expansion at its third mode (k=3) and applying the method of residuals for the separation of the real and imaginary parts of same order [12], we might be able to get the following seven ODEs to represent the overall system:

$$v_t=-{\displaystyle \frac{2}{D}}b{\left({\displaystyle \frac{v_c}{D}}\right)}^d{v}^{2-d}+{\displaystyle \frac{2g\beta }{\pi }}\left[{\alpha }_{1}sin\gamma -{\beta }_1(1-cos\gamma )+{\displaystyle \frac{1}{3}}{\alpha }_{3}sin3\gamma -{\displaystyle \frac{1}{3}}{\beta }_3(1-cos3\gamma )\right]$$

(21)

$${\alpha }_{1t}=-{\displaystyle \frac{a{\pi }^2}{{(L+H)}^2}}{\alpha }_1+{\displaystyle \frac{\pi }{(L+H)}}v{\beta }_1+{\displaystyle \frac{\mathsf{\Gamma }}{2\pi }}sin\gamma$$

(22)

$${\beta }_{1t}=-{\displaystyle \frac{a{\pi }^2}{{(L+H)}^2}}{\beta }_1+{\displaystyle \frac{\pi }{(L+H)}}v{\alpha }_1+{\displaystyle \frac{\mathsf{\Gamma }}{2\pi }}(1+cos\gamma )$$

(23)

$${\alpha }_{2t}=-a{\displaystyle \frac{4{\pi }^2}{{(L+H)}^2}}{\alpha }_2+{\displaystyle \frac{2\pi }{(L+H)}}v{\beta }_2+{\displaystyle \frac{{\mathsf{\Gamma }}_1}{4\pi }}sin2\gamma$$

(24)

$${\beta }_{2t}=-a{\displaystyle \frac{4{\pi }^2}{{(L+H)}^2}}{\beta }_2+{\displaystyle \frac{2\pi }{(L+H)}}v{\alpha }_2+{\displaystyle \frac{{\mathsf{\Gamma }}_1}{4\pi }}(1+cos2\gamma )$$

(25)

$${\alpha }_{3t}=-a{\displaystyle \frac{9{\pi }^2}{{(L+H)}^2}}{\alpha }_3+{\displaystyle \frac{3\pi }{(L+H)}}v{\beta }_3+{\displaystyle \frac{\mathsf{\Gamma }}{6\pi }}sin3\gamma$$

(26)

$${\beta }_{3t}=-a{\displaystyle \frac{9{\pi }^2}{{(L+H)}^2}}{\beta }_3+{\displaystyle \frac{3\pi }{(L+H)}}v{\alpha }_3+{\displaystyle \frac{\mathsf{\Gamma }}{6\pi }}(1+cos3\gamma )$$

(27)

where,

$$\gamma ={\displaystyle \frac{\pi L}{L+H}};\mathsf{\Gamma }={\displaystyle \frac{4}{\rho cD}}\left({\displaystyle \frac{m{C}_p\mathsf{\Delta }T}{\pi DH}}+q\right);{\mathsf{\Gamma }}_1={\displaystyle \frac{4}{\rho cD}}\left({\displaystyle \frac{m{C}_p\mathsf{\Delta }T}{\pi DH}}-q\right)$$

According to our second assumption mentioned earlier in this section, the existence of equilibrium point is ensured. Usually, the first step after formulating the model equations, is to find the equilibrium points followed by linearization of the system around these equilibrium points. The properties that are considered for this study are presented in Table 2.

The Jacobian matrix $\mathbf{A}$ of the system is a $7\times 7$ square matrix, with each element $A_{kj}$ being the partial derivative of the $k^{\mathrm{th}}$ equation with respect to the $j^{\mathrm{th}}$ variable.

$$\mathbf{A}=\left[\begin{array}{ccccccc}-(2-d){\displaystyle \frac{2b}{D}}{\left({\displaystyle \frac{v_c}{D}}\right)}^d{v}_{ss}^{1-d}& {\displaystyle \frac{2g\beta }{\pi }}sin\gamma & -{\displaystyle \frac{2g\beta }{\pi }}(1-cos\gamma )& 0& 0& {\displaystyle \frac{2g\beta }{3\pi }}sin3\gamma & -{\displaystyle \frac{2g\beta }{\pi }}(1-cos3\gamma )\\ {\displaystyle \frac{\gamma }{L}}{\beta }_{1ss}& -a{\left({\displaystyle \frac{\gamma }{L}}\right)}^2& {\displaystyle \frac{\gamma }{L}}{v}_{ss}& 0& 0& 0& 0\\ -{\displaystyle \frac{\gamma }{L}}{\alpha }_{1ss}& -{\displaystyle \frac{\gamma }{L}}{v}_{ss}& -a{\left({\displaystyle \frac{\gamma }{L}}\right)}^2& 0& 0& 0& 0\\ {\displaystyle \frac{2\gamma }{L}}{\beta }_{2ss}& 0& 0& -a{\left({\displaystyle \frac{2\gamma }{L}}\right)}^2& {\displaystyle \frac{2\gamma }{L}}{v}_{ss}& 0& 0\\ -{\displaystyle \frac{2\gamma }{L}}{\alpha }_{2ss}& 0& 0& -{\displaystyle \frac{2\gamma }{L}}{v}_{ss}& -a{\left({\displaystyle \frac{2\gamma }{L}}\right)}^2& 0& 0\\ {\displaystyle \frac{3\gamma }{L}}{\beta }_{3ss}& 0& 0& 0& 0& -a{\left({\displaystyle \frac{3\gamma }{L}}\right)}^2& {\displaystyle \frac{3\gamma }{L}}{v}_{ss}\\ -{\displaystyle \frac{3\gamma }{L}}{\alpha }_{3ss}& 0& 0& 0& 0& -{\displaystyle \frac{3\gamma }{L}}{v}_{ss}& -a{\left({\displaystyle \frac{3\gamma }{L}}\right)}^2\end{array}\right]$$

(28)

Similarly, the matrix $\mathbf{B}$ was calculated as follows:

$$\mathbf{B}={\left[0{\displaystyle \frac{sin\gamma }{2\pi }}{\displaystyle \frac{1+cos\gamma }{2\pi }}00{\displaystyle \frac{sin3\gamma }{6\pi }}{\displaystyle \frac{1+cos3\gamma }{6\pi }}\right]}^T$$

(29)

So, the linearized system in continuous time settings can be written as:

$${\dot{\mathbf{X}}}^{\prime }=\mathbf{A}{\mathbf{X}}^{\prime }+\mathbf{B}{\mathbf{U}}^{\prime }$$

(30)

The above system can be discretized using suitable time step, and represented as follows;

$${\mathbf{X}}_{k+1}^{\prime }={\mathbf{A}}_d{\mathbf{X}}_k^{\prime }+{\mathbf{B}}_d{\mathbf{U}}_k^{\prime }$$

(31)

The controllability matrix is similarly defined as:

$$\mathcal{C}=\left[\begin{array}{ccccc}{B}_d& A_d{B}_d& A_d^2{B}_d& \dots & A_d^{n-1}{B}_d\end{array}\right]$$

(32)

Controllability Condition

The system is said to be controllable if and only if[16]:

$$rank(\mathcal{C})=n$$

where n is the dimension of the state vector x; here in our settings, the system was not found fully controllable, as the rank of the controllability matrix was 5. However, for such systems which are not fully controllable, it is advised to check the stabilizability. From PBH test (procedure details can be seen in Appendix A.1), it was found that all the uncontrolled eigenvalues are stable. The details of eigenvalues through controllability and stability analysis are shown in the following Figure 3. Three of the eigenvalues (green) are stable and controllable, two (pink) are unstable but controllable and rest of the two eigenvalues (blue) are uncontrollable but stable.

4. Discrete-Time Linear Quadratic Regulator (dLQR)

To develop an optimal control strategy for the discrete system, a discrete-time Linear Quadratic Regulator (dLQR) was designed. The discrete-time linearized model, characterized by the state-space matrices $A_d$ and $B_d$ were used for the controller design.

The performance index for the dLQR is defined as:

$$J=\sum _{k=0}^{\infty }\left[{\tilde{x}}_k^{T}Q{\tilde{x}}_k+{\tilde{u}}_k^{T}R{\tilde{u}}_k\right]$$

(33)

Here, Q and R, the weighting matrices are used to penalize state deviations and control effort respectively. In this work, they were selected as:

$$Q=\mathrm{diag}(1000,10,10,1,1,1000,1000),R=10>0.$$

The dLQR gain is obtained by solving the Discrete-Time Algebraic Riccati Equation (DARE) [17]:

$$P=A_d^{T}P{A}_d-A_d^{T}P{B}_d{(R+B_d^{T}P{B}_d)}^{-1}{B}_d^{T}P{A}_d+Q$$

(34)

The matrix P can be computed using the "dare" command in MATLAB. The optimal state-feedback gain matrix K is then given by:

$$K={(R+B_d^{T}P{B}_d)}^{-1}{B}_d^{T}P{A}_d$$

(35)

The corresponding optimal control input is:

$${\tilde{u}}_k^{*}=-K{\tilde{x}}_k$$

(36)

When this control law is applied to the discrete-time state-space model, the closed-loop dynamics become:

$${\tilde{x}}_{k+1}=(A_d-B_{d}K){\tilde{x}}_k$$

(37)

Once the model states have been calculated, the final step will be the reconstruction of temperature function from these model states i.e. ${\alpha }_k(t)$ and ${\beta }_k(t)$ using the following equation.

$$T(\zeta ,t)-T_0=2{\alpha }_1(t)cos\zeta -2{\beta }_1(t)sin\zeta +2{\alpha }_2(t)cos2\zeta -2{\beta }_2(t)sin2\zeta +2{\alpha }_3(t)cos3\zeta -2{\beta }_3(t)sin3\zeta$$

(38)

5. Linear Model Predictive Control (LMPC)

To extend the optimal control framework for this multivariable system, a Linear Model Predictive Control (LMPC) strategy was implemented. MPC is particularly useful for multiple-input multiple-output (MIMO) systems and for directly handling system constraints [18]. These features make MPC a natural choice for this case study.

5.1. Discrete-Time Linear Model

The continuous-time nonlinear dynamics were linearized and discretized with a sampling time of $\mathrm{dt}=0.1$ s, producing the discrete-time state-space model:

$${\mathbf{x}}_{k+1}=A_d{\mathbf{x}}_k+B_d{\mathbf{u}}_k$$

(39)

5.2. MPC Cost Function with Terminal Penalty

The finite-horizon quadratic cost function is defined as:

$$J=\sum _{i=0}^{N-1}\left[{({\mathbf{x}}_i-{\mathbf{x}}_s)}^{T}Q({\mathbf{x}}_i-{\mathbf{x}}_s)+{({\mathbf{u}}_i-{\mathbf{u}}_s)}^{T}R({\mathbf{u}}_i-{\mathbf{u}}_s)\right]+{({\mathbf{x}}_N-{\mathbf{x}}_s)}^{T}P({\mathbf{x}}_N-{\mathbf{x}}_s)$$

(40)

Here: $Q>0$ penalizes state deviations, $R>0$ penalizes control effort and P is the terminal cost matrix obtained from the Discrete Algebraic Riccati Equation (DARE) for the unconstrained infinite-horizon LQR problem. This choice of P guarantees closed-loop stability when constraints are inactive.

5.3. Prediction Model for Quadratic Programming (QP) Formulation

For LMPC framework, built-in program "quadprog" of MATLAB was used [19]. However, lets define the stacked control sequence and the predicted states over the horizon (N) as::

$$\mathbf{U}=\left[\begin{array}{c}{\mathbf{u}}_0\\ {\mathbf{u}}_1\\ ⋮\\ {\mathbf{u}}_{N-1}\end{array}\right],\mathbf{X}=\left[\begin{array}{c}{\mathbf{x}}_1\\ {\mathbf{x}}_2\\ ⋮\\ {\mathbf{x}}_N\end{array}\right]$$

(41)

The states can be expressed in terms of the initial condition ${\mathbf{x}}_k$ and $\mathbf{U}$ as:

$$\mathbf{X}=\mathcal{A}{\mathbf{x}}_k+\mathcal{B}\mathbf{U}$$

(42)

where:

$$\mathcal{A}=\left[\begin{array}{c}{A}_d\\ A_d^2\\ ⋮\\ A_d^N\end{array}\right],\mathcal{B}=\left[\begin{array}{cccc}{B}_d& 0& \dots & 0\\ A_d{B}_d& B_d& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ A_d^{N-1}{B}_d& A_d^{N-2}{B}_d& \dots & B_d\end{array}\right]$$

5.4. Quadratic Program (QP) Matrices H and f

The cost function can be rewritten as [19,20] :

$$J={\displaystyle \frac{1}{2}}{\mathbf{U}}^{T}H\mathbf{U}+{\mathbf{f}}^T\mathbf{U}+\mathrm{const}$$

(43)

where the Hessian matrix H and the linear cost vector $\mathbf{f}$ are given by

$$H=2\left({\mathcal{B}}^T{Q}_{\mathrm{bar}}\mathcal{B}+R_{\mathrm{bar}}\right),\mathbf{f}=2{\mathcal{B}}^T{Q}_{\mathrm{bar}}\left(\mathcal{A}{\mathbf{x}}_k-{\mathbf{x}}_s^{\mathrm{stack}}\right)-2R_{\mathrm{bar}}{\mathbf{u}}_s^{\mathrm{stack}}$$

Here, $\mathcal{A}$ and $\mathcal{B}$ define the stacked state prediction model, and $Q_{\mathrm{bar}}$, $R_{\mathrm{bar}}$ are block-diagonal matrices formed from $Q,R$, and terminal cost P.

where:

$$Q_{\mathrm{bar}}=\mathrm{blkdiag}(Q,Q,\dots ,Q,P),R_{\mathrm{bar}}=\mathrm{blkdiag}(R,R,\dots ,R)$$

and:

$$H=2\left({\mathcal{B}}^T{Q}_{\mathrm{bar}}\mathcal{B}+R_{\mathrm{bar}}\right)$$

$$\mathbf{f}=2{\mathcal{B}}^T{Q}_{\mathrm{bar}}\left(\mathcal{A}{\mathbf{x}}_k-{\mathbf{x}}_s^{\mathrm{stack}}\right)-2R_{\mathrm{bar}}{\mathbf{u}}_s^{\mathrm{stack}}$$

Here, ${\mathbf{x}}_s^{\mathrm{stack}}$ and ${\mathbf{u}}_s^{\mathrm{stack}}$ are the stacked state and input setpoints over the horizon. The details of derivation of Quadratic Cost Function for MPC can be found in Appendix A.2.

5.5. Constraints

In this study, only input constraints were considered:

$$u_{min}\le u_k\le u_{max}$$

(44)

These are converted to a compact inequality constraint:

$$G\mathbf{U}\le \mathbf{w}$$

where G and $\mathbf{w}$ are constructed to represent the bounds for each input over the prediction horizon.

6. Results

This section will be used to display results obtained from the discrete Linear Quadratic Control (dLQR) and Linear Model Predictive Control (LMPC) with and without constraints implemented on the natural circulation pipeline system.

The estimated states and control input, represented as deviation variables, obtained using LQR, unconstrained MPC, and constrained linear MPC are illustrated in Figure 4. The dLQR and unconstrained MPC yield identical responses. However, due to the imposed input constraints (U_min = -0.784 and U_max = 1.518), the constrained MPC exhibits an initial overshoot in the temperature state. Despite this, the system reaches steady state within 65 seconds, and the overall response remains smooth and free from oscillations. $\mathsf{\Gamma }$ from Equation (18) used as control input obtained all the controllers was also plotted in Figure 4.

6.1. Remarks on dLQR vs. MPC

When the constraints are inactive, the optimal solution from Linear MPC coincides with the discrete-time LQR (dLQR) control law [21]:

$${\mathbf{u}}_k=-K{\mathbf{x}}_k$$

where K is obtained from DARE. The primary advantage of LMPC is realized when constraints become active, as the optimizer adjusts future inputs to satisfy them while minimizing the cost, enabling safe and feasible operation under system limits.

The performance of MPC is highly sensitive to initial conditions, with greater deviations from the steady state requiring correspondingly higher control effort. This effect was examined by initiating control after different durations of open-loop operation. Due to the inherent instability of the system, longer delays resulted in the system state moving further from the steady state, thereby increasing the control effort demanded from the controller. The resulting comparison is presented in Figure 5. As a result, physical constraints be violated by conventional LQR or unconstrained MPC, and constrained MPC is very useful tool in considering such constraints.

$\mathsf{\Gamma }$ used as control input, which consists of both coolant flow rate and the power heater. Here in our case, it was assumed that power of the heating in the heating section will remain fixed and coolant flow rate will be manipulated as control input. Such applications with fixed available heat source can be found in geothermal systems. Here, in our case, the estimated flow rate of the coolant as "control input" with fixed heater power (Q) =1800 W is shown in the Figure 5.

In this example, the primary physical constraint is the maximum flow capacity of the coolant pump. As shown in Figure 5, when the initial conditions deviate further from the steady state, the controller must exert greater control effort to stabilize the system. This increased demand can lead to the physical constraint being reached or exceeded, effectively surpassing the pump’s maximum flow capacity. Consequently, a pump with a higher capacity may be required to meet the control objectives under such conditions. One of the significant advantages of employing constrained MPC is that it enables the proper selection of pump capacity, ensuring maximum heat removal from the natural circulation loop while respecting system limitations. Moreover, by explicitly accounting for such constraints during controller design, MPC can prevent operational violations, improve system safety, and optimize equipment sizing for both performance and cost-effectiveness.

7. Discussion

Pipeline systems are fundamental to a wide range of industrial operations, and ensuring their reliability is essential for maintaining both safety and efficiency. Among these systems, natural circulation loops present distinct challenges due to their dependence on buoyancy-driven flow rather than mechanical pumping. This inherent characteristic makes them particularly susceptible to disturbances, where flow instabilities can significantly impact both performance and operational safety.

This study investigates the use of Model Predictive Control (MPC) as an effective approach for removal of heat in natural circulation pipelines used in geothermal systems. Findings show that MPC not only stabilizes system dynamics but also optimizes control efforts, which are closely tied to heating or cooling costs. When compared with conventional methods such as the Linear Quadratic Regulator (LQR), MPC gives an advantage of setting physical constraints on control action.

One of MPC’s key advantages is its capacity to manage input constraints explicitly. This allows for more practical and safer control actions that respect the physical and operational boundaries of system components like valves, actuators, sensors , maximum voltage capacity of the motor/pump, which may exhibit nonlinear behaviors or operational limits. Furthermore, MPC’s predictive capabilities enable it to respond proactively to dynamic conditions such as process fluctuations, thermal feedback, and variable flow rates.

From an implementation perspective, the control framework proposed in this research emphasizes both computational efficiency and ease of deployment, making it well-suited for industrial applications. By utilizing scalable control logic and adaptable system models, the approach effectively translates theoretical control concepts into practical, field-ready solutions.

In conclusion, incorporating MPC into natural circulation pipeline systems represents a promising strategy for improving process reliability and operational safety. Beyond its role as a control method, MPC has the potential to serve as a core component of advanced pipeline management systems.

Appendix A.1. Details of PBH Test for stabilizability [13]

Consider the discrete-time linear system:

$$x_{k+1}=A_d{x}_k+B_d{u}_k,$$

where $A_d\in {\mathbb{R}}^{n\times n}$, $B_d\in {\mathbb{R}}^{n\times m}$.

• Compute the eigenvalues of the system matrix $A_d$: denote them as ${\lambda }_i$.
• For each eigenvalue ${\lambda }_i$ such that $Re({\lambda }_i)\ge 1$, verify controllability using the Popov–Belevitch–Hautus (PBH) test:

  $$rank\left[\begin{array}{cc}{\lambda }_{i}I-A_d& B_d\end{array}\right]=n$$

  where n is the number of states in the system.
• If the PBH test is satisfied for all unstable (or marginally unstable) eigenvalues, then the system is said to be stabilizable.

Appendix A.2. Derivation of Quadratic Cost Function in MPC

The stage cost at time step k is given by

$${\ell }_k={({\mathbf{x}}_k-{\mathbf{x}}_s)}^{T}Q({\mathbf{x}}_k-{\mathbf{x}}_s)+{({\mathbf{u}}_k-{\mathbf{u}}_s)}^{T}R({\mathbf{u}}_k-{\mathbf{u}}_s)$$

(A1)

where $Q⪰0$ and $R\succ 0$ are weighting matrices for states and inputs, respectively.

Define the stacked input vector over the prediction horizon N:

$$\mathbf{U}=\left[\begin{array}{c}{\mathbf{u}}_0\\ {\mathbf{u}}_1\\ ⋮\\ {\mathbf{u}}_{N-1}\end{array}\right]$$

(A2)

and the stacked predicted states:

$$\mathbf{X}=\left[\begin{array}{c}{\mathbf{x}}_1\\ {\mathbf{x}}_2\\ ⋮\\ {\mathbf{x}}_N\end{array}\right]$$

(A3)

The system dynamics in discrete-time linear form are

$${\mathbf{x}}_{k+1}=A_d{\mathbf{x}}_k+B_d{\mathbf{u}}_k$$

(A4)

Stacked over the horizon, the states can be predicted from the initial state ${\mathbf{x}}_0$ and the input sequence $\mathbf{U}$ as

$$\mathbf{X}=\mathcal{A}{\mathbf{x}}_0+\mathcal{B}\mathbf{U}$$

(A5)

where $\mathcal{A}$ and $\mathcal{B}$ are appropriately defined matrices capturing system evolution.

Construct block-diagonal weighting matrices:

$$Q_{\mathrm{bar}}=\mathrm{blkdiag}(Q,Q,\dots ,Q,P),R_{\mathrm{bar}}=\mathrm{blkdiag}(R,R,\dots ,R)$$

(A6)

and stacked reference vectors

$${\mathbf{x}}_s^{\mathrm{stack}}=\left[\begin{array}{c}{\mathbf{x}}_s\\ {\mathbf{x}}_s\\ ⋮\\ {\mathbf{x}}_s\end{array}\right],{\mathbf{u}}_s^{\mathrm{stack}}=\left[\begin{array}{c}{\mathbf{u}}_s\\ {\mathbf{u}}_s\\ ⋮\\ {\mathbf{u}}_s\end{array}\right]$$

(A7)

The total cost over the horizon is

$$J={(\mathbf{X}-{\mathbf{x}}_s^{\mathrm{stack}})}^T{Q}_{\mathrm{bar}}(\mathbf{X}-{\mathbf{x}}_s^{\mathrm{stack}})+{(\mathbf{U}-{\mathbf{u}}_s^{\mathrm{stack}})}^T{R}_{\mathrm{bar}}(\mathbf{U}-{\mathbf{u}}_s^{\mathrm{stack}})$$

(A8)

Substituting $\mathbf{X}=\mathcal{A}{\mathbf{x}}_0+\mathcal{B}\mathbf{U}$ gives

$$J={\left(\mathcal{A}{\mathbf{x}}_0+\mathcal{B}\mathbf{U}-{\mathbf{x}}_s^{\mathrm{stack}}\right)}^T{Q}_{\mathrm{bar}}\left(\mathcal{A}{\mathbf{x}}_0+\mathcal{B}\mathbf{U}-{\mathbf{x}}_s^{\mathrm{stack}}\right)+{(\mathbf{U}-{\mathbf{u}}_s^{\mathrm{stack}})}^T{R}_{\mathrm{bar}}(\mathbf{U}-{\mathbf{u}}_s^{\mathrm{stack}})$$

(A9)

Expanding and grouping terms yields a quadratic form in $\mathbf{U}$:

$$J={\mathbf{U}}^T\left({\mathcal{B}}^T{Q}_{\mathrm{bar}}\mathcal{B}+R_{\mathrm{bar}}\right)\mathbf{U}+2{\mathbf{U}}^T\left({\mathcal{B}}^T{Q}_{\mathrm{bar}}(\mathcal{A}{\mathbf{x}}_0-{\mathbf{x}}_s^{\mathrm{stack}})-R_{\mathrm{bar}}{\mathbf{u}}_s^{\mathrm{stack}}\right)+\mathrm{const}$$

(A10)

Identifying the Hessian matrix and linear term:

$$\begin{array}{|c|}\hline \begin{array}{cc}\hfill H& =2\left({\mathcal{B}}^T{Q}_{\mathrm{bar}}\mathcal{B}+R_{\mathrm{bar}}\right)\hfill \\ \hfill \mathbf{f}& =2\left({\mathcal{B}}^T{Q}_{\mathrm{bar}}(\mathcal{A}{\mathbf{x}}_0-{\mathbf{x}}_s^{\mathrm{stack}})-R_{\mathrm{bar}}{\mathbf{u}}_s^{\mathrm{stack}}\right)\hfill \end{array}\\ \hline\end{array}$$

(A11)

Thus, the MPC cost is equivalently expressed as

$$J={\displaystyle \frac{1}{2}}{\mathbf{U}}^{T}H\mathbf{U}+{\mathbf{f}}^T\mathbf{U}+\mathrm{const}$$

(A12)

which forms the objective function of the quadratic program solved at each control step.