Constraint Optimal Model-Based Disturbance Predictive and Rejection Control Method of a Parabolic Trough Solar Field

Abstract

The control of the field outlet temperature of a parabolic trough solar field (PTSF) is crucial for the safe and efficient operation of the solar power system but with the difficulties arising from the multiple disturbances and constraints imposed on the variables. To this end, this paper proposes a constraint optimal model-based disturbance predictive and rejection control method with a disturbance prediction part. In this method, the steady-state target sequence is dynamically corrected in the presence of constraints, the lumped disturbance, and its future dynamics predicted by the least-squares support vector machine. In addition, a maximum controlled allowable set is constructed in real time to transform an infinite number of constraint inequalities into finite ones with the integration of the corrected steady-state target sequence. On this basis, an equivalent quadratic programming constrained optimization problem is constructed and solved by the dual-mode control law. The simulation results demonstrate the setpoint tracking and disturbance rejection performance of our design under the premise of constraint satisfaction.

1. Introduction

With the promotion of carbon neutralization worldwide, the sustainable development of renewable energy has become a promising direction in China [1]. Solar energy is one of the environmentally friendly renewable energy sources [2]. Parabolic trough concentrating solar power plants stand out among the renewable energy crowd as solar energy is one of the most abundant renewable energy sources [3]. As an energy absorption and conversion center, the effective control of the field outlet temperature of a parabolic trough solar field (PTSF) is crucial to the safe and efficient operation of the whole solar power system [4].

The control objective of the high inertia PTSF is to maintain the field outlet temperature close to the desired setpoint by adjusting the flow rate of heat transfer fluid (HTF) in the tube. However, this control performance suffers from several factors. On the one hand, there are multiple time-varying disturbances in the PTSF, such as direct normal irradiation (DNI), heat transfer fluid fluctuation, and parameter perturbation. On the other hand, the constraints imposed on the control variables limit the variation in the variables. Therefore, how to eliminate the effects of the field disturbances under the constraint scenarios is a difficult but urgent task.

In recent years, model predictive control (MPC) has gained popularity for solar power plant control due to its natural advantages in dealing with large inertia and constraints [5,6]. The MPC strategies applied to the solar power plant are robust [7,8], nonlinear [9,10,11], scheduling [12,13], and machine learning-based form [14]. The traditional MPC strategies do possess good disturbance rejection performance. Nevertheless, the stability of the closed-loop system cannot be guaranteed because optimality is not the same as stability. To this end, Rossiter proposed an optimal model predictive control (OMPC) method in the monograph [15]. The OMPC strategy adopts infinite time domain prediction and dual-mode control law. Thus, the closed-loop stability of the nominal system can be guaranteed. Some variants of OMPC have been applied to PTSF for disturbance rejection [16,17,18,19]. Specifically, the references [16,17] mainly focus on eliminating the effect of direct normal irradiation, while the multiple disturbances are comprehensively considered in [18,19] from the perspective of lumped disturbance estimation and feedforward compensation. The results show that the disturbance compensation property is better in the lumped disturbance aspect, especially when the future dynamics of the lumped disturbance on the prediction horizon are feedforward compensated.

However, there is no mechanism for dealing with constraints imposed on the variables in [16,17,18,19]. In practice, there are upper and lower limits to the heat transfer fluid flow rate for safety purposes. It is well known that constraint satisfaction and disturbance rejection are conflicting control goals. Disturbance rejection usually requires a large fluctuation of the flow rate, which leads to the violation of the upper and lower constraints. Moreover, In the OMPC framework, the system dynamics evolve during the infinite prediction horizon, and the resulting control variables should satisfy the constraint requirement at each prediction horizon, which will generate an infinite number of constraint inequalities. Accordingly, Rossiter developed a maximum controlled allowable set to deal with the infinite number of constraint inequalities produced by an infinite prediction horizon [15]. However, the traditional MCAS strategy fails when there are multiple disturbances under constraints. To this end, how to deal with the infinite number of constraint inequalities by organically fusing the future behavior of the lumped disturbance, and then how to solve the infinite constraint optimization problem are essential to effectively reject the time-varying disturbances under constraints.

Motivated by the above concerns, this paper proposes a constraint optimal model-based disturbance predictive and rejection controller (C-ODPRC) to suppress disturbances while satisfying the given constraints. The main contributions of the composite controller are summarized as follows:

• A real-time correction strategy for steady-state target sequence is designed to modify the equilibrium by compromising constraint satisfaction and disturbance rejection with the integration of the future dynamics of the lumped disturbance.

• A dynamic construction method of a maximum controlled allowable set is proposed with the disturbance prediction part based on the obtained steady-state target sequence to convert an infinite number of constraint inequalities into finite ones.

• The infinite horizon constrained optimization problem is converted into an equivalent quadratic programming problem, which is then solved by dual-mode control law and the constructed maximum controlled allowable set.

The proposed design exhibits superiorities in target tracking and disturbance rejection performance under the constrained scenario as the future dynamics of the lumped disturbance are organically integrated to provide more feedforward signals.

The rest of the paper is organized as follows. The analytical and control model of the PTSF is given in Section 2. Section 3 presents the methods of the proposed composite controller. Simulation studies are carried out in Section 4 to verify the superiorities of our proposed controller in disturbance rejection and set-point tracking. The last section concludes this paper.

2. Modeling of Parabolic Trough Solar Field

2.1. Analytical Model and Model Validation

The schematic diagram of the studied parabolic trough solar field (PTSF) is shown in Figure 1.

As shown in Figure 1, the PTSF consists of numerous parallel trough solar collectors (PTCs), each of which consists of an absorber tube, a glass cover, and a reflector. In practice, the direct normal irradiation (DNI) passes through the glass cover and is then absorbed by the heat transfer fluid (HTF) flowing inside the absorber tube. The heated heat transfer fluid is then sent to the power generation module or the heat storage subsystem.

In general, the analytical model can be represented by a single PTC (including an absorber tube, a glass cover, and a reflector), as the PTSF is made up of many PTCs of the same structure. Moreover, the length of each PTC can be divided into n segments on average. Thus, a distributed analytical parameter model can be established for each segment, in which the heat transfer process of the jth segment can be expressed by [20] and shown as follows:

(1)$\left\{\begin{array}{l}{\rho }_H{c}_H{A}_{a,i}{\displaystyle \frac{d{T}_H(j)}{dt}}=h_H{P}_{a,i}[T_a(j)-T_H(j-1)]-{\displaystyle \frac{q_H{c}_H}{\Delta l}}[T_H(j)-T_H(j-1)]\\ {\rho }_a{c}_a{A}_a{\displaystyle \frac{d{T}_a(j)}{dt}}=I_n\cdot \eta \cdot w-h_H{P}_{a,i}[T_a(j)-T_H(j)]-Q_{rad}\\ {\rho }_g{c}_g{A}_g{\displaystyle \frac{d{T}_g(j)}{dt}}=Q_{rad}-\sigma {\epsilon }_g{P}_{g,o}\{{[T_g(j)+273.15]}^4-{[T_{sky}+273.15]}^4\}-h_g{P}_{g,o}[T_g(j)-T_{atm}]\end{array}\right.$

where the subscripts H, g, a, o, i, sky, and atm represent heat transfer fluid, glass cover, absorber tube, the tube outside, the tube inside, sky, and atmosphere, respectively; T_H(j), T_a(j), T_g(j) are the temperature of the heat transfer fluid, absorber tube, and glass cover at the outlet of the jth section, respectively; T_H(n) = T_H,out (T_H(0) = T_H,in) is the field outlet (inlet) temperature; q_H, the flow rate of heat transfer fluid, I_n, the heat of direct normal irradiation, ρ, the density, c, the heat capacity, A, the cross-section area, P, the diameter, Δl, axial length of each section, h, convective heat transfer coefficient, w, aperture width, ε_a, absorber wall emissivity, ε_g, glass emissivity, Q_rad, radiant heat between the absorber tube and glass cover, η, optical efficiency. Q_rad and η can be computed by

(2)$\left\{\begin{array}{l}{Q}_{rad}=P_{a,o}\{{[T_a(j)+273.15]}^4-{[T_g(j)+273.15]}^4\}\sigma /[{\displaystyle \frac{1}{{\epsilon }_a}}+{\displaystyle \frac{1-{\epsilon }_g}{{\epsilon }_g}}({\displaystyle \frac{r_{a,o}}{r_{g,i}}})]\\ {\eta }_{opt}={\prod }_{k=1}^8{\eta }_k\end{array}\right.$

Parameters in Equations (1) and (2) are listed in Table 1.

In practice, the field outlet temperature, T_H,out, is crucial to the balance between the economy and safety of the PTSF. Operators usually maintain the T_H,out close to a desired setpoint by adjusting the heat transfer fluid flow rate, q_H. Hence, the T_H,out and q_H are defined as the control output and control input, respectively. However, the control performance is often affected by multiple disturbances. The first type is external disturbances, like I_n, the fluctuation of q_H and T_H,in. The second kind is internal disturbances, including parameter perturbation owing to the varying h_H and h_g and model mismatch. For technical safety, q_H should be limited to a safe range.

The analytical PTSF model (1) is supposed to be validated for its accuracy on dynamic characteristics for simulation purposes. A commonly used method is to compare the analytical model output and the measured output of a real thermal power plant under the same input variables. Guided by this, the model (1) has been validated in [19,20].

Borrowed from [19,20], it can be concluded that the field outlet temperature curve produced by the analytical model (1) can track the measured output with the same direct normal irradiation, field inlet temperature, and heat transfer fluid flow rate, which means that the dynamic characteristics of the analytical model (1) are close to the real PTSF. Thus, the analytical model (1) is accurate enough to be used as a simulation model and to verify the proposed controller.

2.2. Control Model of PTSF

In order to simplify the controller design process, a state space model with only one lumped disturbance is expected to be constructed as the control model based on the analytical model (1). According to model (1), only the nth dynamic characteristic of the PTSF describes the relationship between the control output and input, T_H,out and q_H, which can be expressed as

(3)${\displaystyle \frac{d{T}_{H,out}}{dt}}=a_2[T_a(n)-T_{H,out}]-a_1{q}_H[T_{H,out}-T_H(n-1)]$

where a₂ = h_HP_a,i/(ρ_Hc_HA_a,i), a₁ = 1/(ρ_HA_a,iΔl).

Thus, only (3) is used to construct the control model. The dynamic characteristic (3) can be further linearized based on the typical operating conditions exhibited in Table 2.

(4)$\left\{\begin{array}{l}\dot{x}=A_{m}x+B_{m}u+B_{dm}[(B_{m,o}-B_m)u+g_o]\\ y=C_{m}x\end{array}\right.$

where x = T_H,out − T_H,out,r, u = q_H − q_H,r, A_m = −a₂ − a₁q_H,r, B_dm = 1, C_m = 1, g_o is unknown dynamic, B_m,o = −a₁[T_H,out,r − T_H,r(n − 1)]. Since it is difficult to measure the T_H,r(n − 1) in B_m,o, so its approximation, B_m, is given here instead. The error (B_m,o − B_m)u as well as g_o can be reconstructed into a lumped single one, namely d = (B_m,o − B_m)u + g_o.

Subtracting the linearized model (4) from (3) gives d as

(5)$d=a_2[T_a(n)-T_{H,out}]-a_1{q}_H[T_{H,out}-T_H(n-1)]-A_{m}x-B_{m}u$

According to the zero-order holding principle [21], model (4) can be discretized as

(6)$\left\{\begin{array}{l}{x}_{k+1}=A{x}_k+B{u}_k+B_d{d}_k\\ y_k=C{x}_k\end{array}\right.$

where x_k, u_k, y_k, and d_k are the discretized values of x, u, y, and d, $A=e^{T_s}$, $B={\displaystyle {\int }_0^{T_s}{e}^{A_m{T}_s}{B}_{m}dt}$, $B_d={\displaystyle {\int }_0^{T_s}{e}^{A_m{T}_s}{B}_{dm}dt}$, $C=C_m$, T_s is sampling time.

3. Methods

3.1. Lumped Disturbance Estimation and Prediction

The future dynamics of the lumped disturbance are aimed to be predicted in this subsection. For this purpose, its estimate has to be obtained first, to which we turn next.

3.1.1. Lumped Disturbance Estimation by ESO

An extended state observer (ESO) is designed to estimate the lumped disturbance, d. Rewriting d as the new state x₂, one can obtain an expanded form of model (4) as

(7)$\left\{\begin{array}{l}{\dot{x}}_1=A_m{x}_1+B_{m}u+x_2\\ {\dot{x}}_2=g\\ y=x_1\end{array}\right.$

where g is the first derivative of d. On the basis of model (7), ESO can be designed as

(8)$\left\{\begin{array}{l}{\dot{z}}_1=A_m{z}_1+B_{m}u+z_2+{\beta }_{01}(y-z_1)\\ {\dot{z}}_2={\beta }_{02}(y-z_1)\end{array}\right.$

where z₁ and z₂ are the estimates of y and d; β₀₁ and β₀₂ are the observer gains and can be determined by the pole placement method [22,23] to hold the stability of the ESO. With the ESO (8) at hand, one can obtain the estimate of d at each sampling instant, namely ${\widehat{d}}_k$ = z₂, which is the benchmark for predicting the future dynamics of d. The prediction process will be examined in the next subsection.

3.1.2. Lumped Disturbance Prediction Based on LS-SVM

The least-squares support vector machine (LS-SVM) is a promising data-driven method with good fitting capability and fast computational speed [24]. In this subsection, an LS-SVM based predictor is designed to obtain the future values of the lumped disturbance adaptively based on the estimated lumped disturbance values at past and current moments. More precisely, the training inputs, d_in, contain several consecutive past values, and the training outputs, d_out, contain a value at the corresponding next moment. The test data, d_new, contains consecutive values from the current moment to the past. More precisely,

(9)$d_{in}=\left[\begin{array}{l}{d}_{in,1}^{\mathrm{T}}\\ d_{in,2}^{\mathrm{T}}\\ ⋮\\ d_{in,N}^{\mathrm{T}}\end{array}\right]=\left[\begin{array}{l}\Delta {\widehat{d}}_{j-n_p}, \Delta {\widehat{d}}_{j-n_p+1}, \dots , \Delta {\widehat{d}}_{j-1}\\ \Delta {\widehat{d}}_{j-n_p+1}, \Delta {\widehat{d}}_{j-n_p+2}, \dots , \Delta {\widehat{d}}_j\\ ⋮ ⋮ ⋮\\ \Delta {\widehat{d}}_{j-n_p+N-1} \dots \Delta {\widehat{d}}_{j+N-2}\end{array}\right], d_{out}=\left[\begin{array}{c}{d}_{out,1}^{\mathrm{T}}\\ d_{out,2}^{\mathrm{T}}\\ ⋮\\ d_{out,N}^{\mathrm{T}}\end{array}\right]=\left[\begin{array}{l}\Delta {\widehat{d}}_j\\ \Delta {\widehat{d}}_{j+1}\\ ⋮\\ \Delta {\widehat{d}}_{j+N-1}\end{array}\right]$

(10)$d_{new}=[\Delta {\widehat{d}}_{j-n_p+1},\Delta {\widehat{d}}_{j-n_p+2},\dots , \Delta {\widehat{d}}_j]$

where $\Delta {\widehat{d}}_j={\widehat{d}}_j-{\widehat{d}}_{j-1}$ is the increment, n_p is the length of training and test input data, and the superscript ‘T’ denotes the transpose.

The training and test data shown in (9) and (10) are all initialized to zero matrices at the start of the simulation. The training and test data are updated at each sampling instant adaptively. More specifically, at the kth sampling interval, one can obtain the ${\widehat{d}}_k$ by ESO (8) and calculate the $\Delta {\widehat{d}}_k$. Followed by updating the support vector as $d_{in,up}^{\mathrm{T}}=[\Delta {\widehat{d}}_{k-n_p},\Delta {\widehat{d}}_{k-n_p+1}\dots , \Delta {\widehat{d}}_{k-1}],$ $d_{out,up}=\Delta {\widehat{d}}_k.$ Next, update the training data by adding the (d_in,up, d_out,up) to its last row and removing its first row:

(11)$d_{in}^{\mathrm{T}}=[d_{in,2},\dots ,d_{in,N},d_{in,up}], d_{out}^{\mathrm{T}}=[d_{out,2},\dots ,d_{out,N},d_{out,up}]$

Then, the LS-SVM parameters α and b can be computed by solving the following optimization problem,

(12)$\begin{array}{ll}\min & {\displaystyle \frac{1}{2}}{w}^{\mathrm{T}}w+\lambda {\displaystyle \sum _{i=1}^N{e}_i^2}\\ \mathrm{s}.\mathrm{t}.& d_{out,i}=w^{\mathrm{T}}\cdot \phi (d_{\mathrm{in},i})+b+e_i, i=1,2,\cdots ,N\end{array}$

where φ(·) is the kernel function, λ is the regulation parameter, and e_i is the errors. One can construct the Lagrange equation as

(13)$L\left(w,b,e_i,{\alpha }_i\right)={\displaystyle \frac{1}{2}}{w}^{\mathrm{T}}w+\lambda {\displaystyle \sum _{i=1}^N{e}_i^2+}{\displaystyle \sum _{i=1}^N{\alpha }_i[d_{out,i}-w^{\mathrm{T}}\cdot \phi (d_{in,i})-b-e_i]}$

Based on the Karush-Kuhn-Tucker conditions, the partial derivatives of L(·) on w, b, e_i, and α_i should be zero, which are also the necessary conditions of optimal solutions. Then, the b and α can be obtained by further rearranging the necessary conditions as

(14)$\left\{\begin{array}{l}b={\displaystyle \frac{E^{\mathrm{T}}{H}^{-1}{d}_{out}}{E^{\mathrm{T}}{H}^{-1}E}}\\ \alpha =H^{-1}{d}_{out}-H^{-1}E\cdot b\end{array}\right.$

where E = [1, 1, …, 1]^T, α = [α₁, α₂, …, α_N] ^T, H is kernel function matrix.

Thus, the future behavior of the lumped disturbance can be predicted via (15) with the latest test input $d_{new}=[\Delta {\widehat{d}}_{k-n_p+1},\dots ,\Delta {\widehat{d}}_{k-1}, \Delta {\widehat{d}}_k]$ and (b, α).

(15)$\Delta {\widehat{d}}_{k+1}=f(d_{new})={\displaystyle \sum _{i=1}^N{\alpha }_{i}K\left(d_{in,i},d_{new}\right)+b}$

where K(·) is the Gaussian kernel function [25].

Update the test input as $d_{new}=[\Delta {\widehat{d}}_{k-n_p+2},\dots ,\Delta {\widehat{d}}_k, \Delta {\widehat{d}}_{k+1}]$ and one can obtain $\Delta {\widehat{d}}_{k+2}$ via (15). Repeat the above steps recursively, the $\Delta {\widehat{d}}_{k+3}$, …, $\Delta {\widehat{d}}_{k+n_a}$ can be obtained recursively.

Finally, convert the incremental values into their absolute form using

(16)${\widehat{d}}_{k+t}={\displaystyle {\sum }_{j=1}^t\Delta {\widehat{d}}_{k+j}}+{\widehat{d}}_k$

Assume that the first n_a-step dynamics of lumped disturbance are considered in the infinite prediction time domain of OMPC, namely

(17)${\widehat{d}}_{k+i}=\left\{\begin{array}{l}\mathrm{variable}, i=0, 1, \dots , n_a\\ {\widehat{d}}_{k+n_a}, i=n_a+1,\dots ,\infty \end{array}\right.$

The future disturbance dynamics at the kth interval can be described by ${\underrightarrow{\widehat{d}}}_{k+1}={[{\widehat{d}}_{k+1}, {\widehat{d}}_{k+2}, \dots , {\widehat{d}}_{k+n_a}]}^{\mathrm{T}}$.

3.2. Dynamic Correction of Steady-State Target Sequence

In order to integrate the obtained $[{\widehat{d}}_k,{\widehat{d}}_{k+1}, {\widehat{d}}_{k+2}, \dots , {\widehat{d}}_{k+n_a}]$ into the framework of the OMPC with the consideration of constraint satisfaction, a dynamic correction strategy of the steady-state target sequence, denoted as (x_ss|k+i, u_ss|k+i), is proposed. The correction strategy aims to correct the original equilibrium of the PTSF into a new one in the presence of disturbance dynamics and constraints imposed on the variables in real time.

Suppose the following constraints exist in the PTSF

(18)$\left\{\begin{array}{l}{u}_{\min }\le u_k\le u_{\max } \\ x_{\min }\le x_k\le x_{\max }\end{array}\right.$

where u_min (x_min) and u_max (x_max) are the upper and lower limits of the control input (state) variables.

Considering that the constraints of the variables are imposed for the safety purpose, which is the most important in engineering applications. Therefore, the principle of the dynamic correction strategy is to ensure constraint satisfaction first, followed by anti-disturbance and setpoint tracking capabilities.

Motivated by the above principle, the steady-state target sequence (x_ss|k+i, u_ss|k+i) can be modified by minimizing the Euclidean distance between the steady-state output sequence y_a|k+i and the setpoint y_r under the given constraints, namely

(19)$\begin{array}{l}\underset{{\underrightarrow{x}}_{ss|k}, {\underrightarrow{u}}_{ss|k-1}}{\min }{\displaystyle \sum _{i=0}^{n_a}{(y_{a|k+i}-y_r)}^T{Q}_s(y_{a|k+i}-y_r)+{(u_{ss|k+i-1}-u_r)}^T{R}_s(u_{ss|k+i-1}-u_r)}\\ \mathrm{(a)}\hfill & s.t. \left\{\begin{array}{l}{u}_{\min }+\epsilon \le u_{ss|k+i-1} \le u_{\max }-\epsilon \\ x_{\min }+\epsilon \le x_{ss|k+i}\le x_{\max }-\epsilon \end{array}\right.\hfill \mathrm{(b)}\hfill & \hspace{1em} \left\{\begin{array}{l}{x}_{ss|k+i}=A{x}_{ss|k+i}+B{u}_{ss|k+i-1}+B_d{\widehat{d}}_{k+i}\\ y_{a|k+i}=C{x}_{ss|k+i}\end{array}\right.\hfill \end{array}$

where ${\underrightarrow{x}}_{ss|k}^{\mathrm{T}}=[x_{ss|k}, x_{ss|k+1},\cdots , x_{ss|k+n_a}]$ and ${\underrightarrow{u}}_{ss|k-1}^{\mathrm{T}}=[u_{ss|k-1}, u_{ss|k},\cdots , u_{ss|k+n_a-1}]$ are the compact expression of the steady-state target sequence, and Q_s and R_s are the weight matrix.

The steady-state target sequence solver (19) can convert the effects of filled disturbances into the new equilibrium (${\underrightarrow{x}}_{\mathrm{ss}|k}$, ${\underrightarrow{u}}_{\mathrm{ss}|k-1}$) with constraints. As a result, the disturbance rejection and setpoint tracking problems under constraints are transformed into steady-state target sequence tracking problems.

It is noticeable that the obtained steady-state target sequence (${\underrightarrow{x}}_{\mathrm{ss}|k}$, ${\underrightarrow{u}}_{\mathrm{ss}|k-1}$) can hold constraint satisfaction at the cost of sacrificing certain anti-disturbance and setpoint tracking capabilities. Therefore, an infinite prediction time domain constraint optimization problem based on (${\underrightarrow{x}}_{\mathrm{ss}|k}$, ${\underrightarrow{u}}_{\mathrm{ss}|k-1}$) can be constructed as

(20)$\begin{array}{l}\underset{{\underrightarrow{u}}_k}{\underbrace{\min }}{J}_k^c={\displaystyle \sum _{i=0}^{\infty }{\overline{x}}_{k+i+1}^{\mathrm{T}}}Q{\overline{x}}_{k+i+1}+u_{k+i}^{\mathrm{T}}R{\overline{u}}_{k+i}\\ \hfill \mathrm{(a)}& \mathrm{s}.\mathrm{t}.\left\{\begin{array}{l}{x}_{k+i+1}=A{x}_{k+i}+B{u}_{k+i}+B_{\mathrm{d}}{d}_{k+i}\\ y_{k+i}=C{x}_{k+i}\end{array}\right.\hfill \hfill \mathrm{(b)}& \hspace{1em} \left\{\begin{array}{l}{u}_{\min }\le u_{k+i}\le u_{\max }\\ x_{\min }\le x_{k+i}\le x_{\max }\end{array}\right.\hfill \end{array}$

where ${\overline{x}}_{k+i}=x_{k+i}-x_{\mathrm{ss}|k+i}$, ${\overline{u}}_{k+i}=u_{k+i}-u_{\mathrm{ss}|k+i}$, ${\underrightarrow{u}}_k={[u_k,\dots ,u_{k+n_{\mathrm{c}}-1},u_{k+n_{\mathrm{c}}-1},\dots ]}^{\mathrm{T}}$.

The key to solving problem (20) lies in dealing with the infinite number of constraint inequalities (20b) derived from the infinite time domain prediction of OMPC.

3.3. The Construction of MCAS with Disturbance Prediction Compensation

A maximum controlled allowable set (MCAS) with disturbance prediction compensation is constructed to address the infinite constraint inequalities based on the obtained steady-state target sequence.

More precisely, the infinite prediction time domain of OMPC is divided into transient (mode 1) and asymptotical (mode 2) predictions and matched with a dual-mode control law in deviation form. When the first n_a step dynamics of the lumped disturbance are considered, the dual-mode control law can be described as

(21)$u_{k+i}-u_{ss|k+i-1}=\left\{\begin{array}{lll}-K(x_{k+i}-x_{ss|k+i})+c_{k+i},& i=0, \dots , n_c-1& (\mathrm{mode} 1) \\ -K(x_{k+i}-x_{ss|k+i}),& i=n_c, \dots , \infty & (\mathrm{mode} 2)\end{array}\right.$

where n_c (n_a ≤ n_c) is the step of mode 1, c_k+i is the control degree of freedom, and K signifies feedback control gain. Combine Equations (6), (19b) and (21) can obtain the one-step dynamic of the PTSF in deviation form:

(22)${\overline{x}}_{k+1}=\mathrm{\Phi }\cdot {\overline{x}}_k+B\cdot c_k+x_{\mathrm{ss}|k}-x_{\mathrm{ss}|k+1}$

where Φ = A − B·K. The dynamic characteristic shown in (22) is the basis for constructing the maximum controlled allowable set with disturbance prediction compensation.

Define new variables as $s_k^{\mathrm{T}}=\left[{\overline{x}}_k^{\mathrm{T}},{\underrightarrow{c}}_k^{\mathrm{T}},{\underrightarrow{x}}_{\mathrm{ss}|k}^{\mathrm{T}}\right]$, the free-response dynamic form of the PTSF can be derived as follows

(23)$\underset{s_{k+1}}{\underbrace{\left[\begin{array}{c}{\overline{x}}_{k+1}\\ {\underrightarrow{c}}_{k+1}\\ {\underrightarrow{x}}_{ss|k+1}\end{array}\right]}}=\underset{\mathrm{\Omega }}{\underbrace{\left[\begin{array}{ccc}\mathrm{\Phi },& [B,0,\dots ,0],& [I,-I,0,L,0]\\ 0,& I_E,& 0\\ 0,& 0,& I_s\end{array}\right]}}\underset{s_k}{\underbrace{\left[\begin{array}{c}{\overline{x}}_k\\ {\underrightarrow{c}}_k\\ {\underrightarrow{x}}_{ss|k}\end{array}\right]}}$

where $I_E=\left[\begin{array}{cccc}0,& I,& 0& \cdots \\ 0,& 0,& I& \cdots \\ ⋮& ⋮& ⋮& \ddots \\ 0,& 0,& \cdots ,& 0\end{array}\right], I_S=\left[\begin{array}{cccc}0,& I,& 0& \cdots \\ 0,& 0,& I& \cdots \\ ⋮& ⋮& ⋮& \ddots \\ 0,& 0,& \cdots ,& I\end{array}\right]$.

Furthermore, the inequality constraints (19a) in the infinite prediction time domain can be organized as

(24)$\underset{G}{\underbrace{\left[\begin{array}{c}-K_s\\ K_s\\ {\Gamma }_x\\ -{\Gamma }_x\end{array}\right]}}\cdot s_{k+i}\le \underset{h}{\underbrace{\left[\begin{array}{c}{u}_{\max }\\ -u_{\min }\\ x_{\max }\\ -x_{\min }\end{array}\right]+\left[\begin{array}{cc} & -u_{ss|k+i-1}\\ & u_{ss|k+i-1}\\ & -x_{ss|k+i}\\ & x_{ss|k+i}\end{array}\right]}}, i=0, 1, \cdots , \infty$

where ${\Gamma }_x=[I,0,0]$, $K_s=\left[K,[-I,0,\dots ,0],0\right]$.

The inequality set (19) represents that there exists an infinite number of inequality constraints over the whole prediction time domain. Combining (23) and (24) gives the following representation

(25)$\left[\begin{array}{c}G\\ G\mathrm{\Omega }\\ ⋮\\ G{\mathrm{\Omega }}^n\\ ⋮\end{array}\right]s_k\le \left[\begin{array}{c}h\\ h\\ ⋮\\ h\\ ⋮\end{array}\right]$

The constraint set shown in (25) can be replaced by its first (n + 1) constraint inequalities, where ‘n’ can be calculated by the following steps.

• ①. Define n = 0, F = G, t = h.

• ②. Solve the following linear programming problem

  (26)$\underset{s_k}{\max }G{\mathrm{\Omega }}^{n+1}{s}_k s.t. F{s}_k\le t$

• ③. If the solution in step ② violates constraint set (24), it means that the holding of the previous (n + 1) inequalities does not guarantee that (GΩⁿ⁺¹) s_k ≤ h also holds. In this case, let n = n + 1, $F=\left[\begin{array}{c} F\\ G{\mathrm{\Omega }}^n\end{array}\right]$, $t=\left[\begin{array}{c} t\\ h\end{array}\right]$. Then, to step ②, otherwise, to step ④.

• ④. If the solution in step ② satisfies the constraint set (24), it is proved that ‘n’ is large enough.

The final value of ‘n’ can be found using steps ①–④. The MCAS can be further expressed as an explicit form of the control degree of freedom, i.e.,

(27)$\underset{F_{\mathrm{s}}}{\underbrace{\left[M_{\mathrm{s}} N_{\mathrm{s}} V_{\mathrm{s}}\right]}}\underset{s_k}{\underbrace{\left[\begin{array}{c} {\overline{x}}_k\\ {\underrightarrow{c}}_k\\ {\underrightarrow{x}}_{\mathrm{ss}|k}\end{array}\right]}}\le t_{\mathrm{s}}$

Then, the maximum controlled allowable set with disturbance prediction compensation can be rewritten as

(28)$S_{\mathrm{MCAS}}=\{\overline{x}:\exists {\underrightarrow{c}}_k s.t. N_{\mathrm{s}}{\underrightarrow{c}}_k\le t_{\mathrm{s}}-M_{\mathrm{s}} {\overline{x}}_k-V_{\mathrm{s}} {\underrightarrow{x}}_{\mathrm{ss}|k}\}$

The maximum controlled allowable set (28) can convert infinite constraint inequalities into finite ones equivalently in the presence of multiple disturbances, which is essentially the feasible domain of problem (20).

3.4. Solving of Infinite Time Domain Constraint Optimization Problem

The infinite time domain constraint optimization problem (20) can be divided into transient zone (J_1k) and asymptotical zone (J_2k), namely

(29)$J_k^c=\underset{J_{1k}}{\underbrace{{\displaystyle \sum _{i=0}^{n_{\mathrm{c}}-1}Q{‖{\overline{x}}_{k+i+1}‖}_2^2}+R{‖{\overline{u}}_{k+i}‖}_2^2}} +\underset{J_{2k}}{\underbrace{{\displaystyle \sum _{i=0}^{\infty }Q{‖{\overline{x}}_{k+n_{\mathrm{c}}+i+1}‖}_2^2}+R{‖{\overline{u}}_{k+n_{\mathrm{c}}+i}‖}_2^2 }}$

In the transient zone, substituting (21) into (22) gives the following expression

(30)$\left\{\begin{array}{l}{\underrightarrow{\overline{x}}}_{k+1}=P_{xs}{x}_k+F_{xs}{\underrightarrow{x}}_{ss|k}+H_{xs}{\underrightarrow{c}}_k\\ {\underrightarrow{\overline{u}}}_k=P_{us}{x}_k+F_{us}{\underrightarrow{x}}_{ss|k}+T_{us}{\underrightarrow{u}}_{ss|k-1}+H_{us}{\underrightarrow{c}}_k\end{array}\right.$

where $P_{xs}=\left[\begin{array}{c} \mathrm{\Phi }\\ {\mathrm{\Phi }}^2\\ ⋮\\ {\mathrm{\Phi }}^{n_a+1}\\ ⋮\\ {\mathrm{\Phi }}^{n_c}\end{array}\right], F_{xs}=\left[\begin{array}{cccc}I-\mathrm{\Phi },& -I,& \dots ,& 0\\ \mathrm{\Phi }(I-\mathrm{\Phi }),& I-\mathrm{\Phi },& \dots ,& 0\\ ⋮& ⋮& & ⋮\\ {\mathrm{\Phi }}^{n_a}(I-\mathrm{\Phi }),& {\mathrm{\Phi }}^{n_a-1}(I-\mathrm{\Phi }),& \dots ,& -\mathrm{\Phi }\\ ⋮& ⋮& & ⋮\\ {\mathrm{\Phi }}^{n_c-1}(I-\mathrm{\Phi }),& {\mathrm{\Phi }}^{n_c-2}(I-\mathrm{\Phi }),& \dots ,& -{\mathrm{\Phi }}^{n_c-n_a}\end{array}\right],H_{xs}=\left[\begin{array}{ccc}B,& 0,\hfill & \dots ,0\\ \mathrm{\Phi }B,& B,\hfill & \dots ,0\\ ⋮& ⋮\hfill & ⋮\\ {\mathrm{\Phi }}^{n_a}B,& {\mathrm{\Phi }}^{n_a-1}B,\hfill & \dots ,0\\ ⋮& ⋮& ⋮\\ {\mathrm{\Phi }}^{n_c-1}B,& {\mathrm{\Phi }}^{n_c-2}B,& \dots ,B\end{array}\right],$ $P_{us}=\left[\begin{array}{c}-K\\ -K\mathrm{\Phi }\\ ⋮\\ -K{\mathrm{\Phi }}^{n_a}\\ ⋮\\ -K{\mathrm{\Phi }}^{n_c-1}\end{array}\right],F_{us}=\left[\begin{array}{cccc}K,& 0,& \dots ,& 0\\ -K(I-\mathrm{\Phi }),& K,& \dots ,& 0\\ ⋮& ⋮& & ⋮& \\ -K{\mathrm{\Phi }}^{n_a-1}(I-\mathrm{\Phi }),& -K{\mathrm{\Phi }}^{n_a-2}(I-\mathrm{\Phi }),& \dots ,& K\\ ⋮& ⋮& & ⋮\\ -K{\mathrm{\Phi }}^{n_c-2}(I-\mathrm{\Phi }),& -K{\mathrm{\Phi }}^{n_c-3}(I-\mathrm{\Phi }),& \dots ,& K{\mathrm{\Phi }}^{n_c-n_a-1}\end{array}\right],T_{us}=\left[\begin{array}{rrrrr}I,& -I,& 0,& \dots ,& 0\\ 0,& I,& -I,& \dots ,& 0\\ ⋮& ⋮& ⋮& & ⋮\\ 0,& 0,& 0,& \dots ,& 0\\ ⋮& ⋮& ⋮& & ⋮\\ 0,& 0,& 0,& \dots ,& 0\end{array}\right],$ $H_{us}=\left[\begin{array}{ccc}I,& 0,& \hfill \dots ,0\\ -KB,& I,& \hfill \dots ,0\\ ⋮& ⋮& \hfill ⋮\\ -K{\mathrm{\Phi }}^{n_a-1}B,& -K{\mathrm{\Phi }}^{n_a-2}B,& \hfill \dots ,0\\ ⋮& ⋮& \hfill ⋮\\ -K{\mathrm{\Phi }}^{n_c-2}B,& -K{\mathrm{\Phi }}^{n_c-3}B,& \hfill \dots , I\end{array}\right], {\underrightarrow{\overline{x}}}_{k+1}=\left[\begin{array}{c}{\overline{x}}_{k+1}\\ {\overline{x}}_{k+2}\\ ⋮\\ {\overline{x}}_{k+n_a}\\ ⋮\\ {\overline{x}}_{k+n_c}\end{array}\right], {\underrightarrow{\overline{u}}}_k=\left[\begin{array}{c}{\overline{u}}_k\\ {\overline{u}}_{k+1}\\ ⋮\\ {\overline{u}}_{k+n_a}\\ ⋮\\ {\overline{u}}_{k+n_c-1}\end{array}\right].$ One can obtain the compact form of J_1k coupled with (29) and (30) as

(31)$J_{1k}={\underrightarrow{\overline{x}}}_{k+1}^T\cdot Q_t\cdot {\underrightarrow{\overline{x}}}_{k+1}+{\underrightarrow{\overline{u}}}_k^T\cdot R_t\cdot {\underrightarrow{\overline{u}}}_k$

where Q_t = diag(Q), R_t = diag(R).

In the asymptotical zone, the lumped disturbance is regarded as a constant and $c_{k+n_c+i}$ = 0, (i = 0, 1, …, ∞). Therefore, J_2k in Equation (27) can be further arranged as

(32)$J_{2k}={\overline{x}}_{k+n_c}^{\mathrm{T}}\underset{P_s}{\underbrace{\{{\displaystyle \sum _{i=0}^{\infty }{({\mathrm{\Phi }}^i)}^{\mathrm{T}}[{\mathrm{\Phi }}^{\mathrm{T}}Q\mathrm{\Phi }+K^{\mathrm{T}}RK]} {\mathrm{\Phi }}^i\}}}{\overline{x}}_{k+n_c}$

where P_s is a Lyapunov function. ${\overline{x}}_{k+n_c}$ can be expressed as

(33)${\overline{x}}_{k+n_c}=P_{sn}{x}_k+F_{sn}{\underrightarrow{x}}_{ss|k}+H_{sn}{\underrightarrow{c}}_k$

where P_sn, F_sn, and H_sn are the last rows of P_s, F_s, and H_s, respectively.

Substituting (30)–(33) into (29), the index $J_k^c$ can be denoted as

(34)$\begin{array}{l}{J}_k^c={\underrightarrow{c}}_k^T\underset{S_{cs}}{\underbrace{[H_s^T\cdot Q_s\cdot H_s+H_{us}^T\cdot R_s\cdot H_{us}+H_{sn}^T\cdot P\cdot H_{sn}]}}{\underrightarrow{c}}_k+2{\underrightarrow{c}}_k^T\underset{S_{xs}}{\underbrace{[H_s^T\cdot Q_s\cdot P_s+H_{us}^T\cdot R_s\cdot P_{us}+H_{sn}^T\cdot P\cdot P_{sn}]}}{x}_k\\ +2{\underrightarrow{c}}_k^T\underset{S_{ss}}{\underbrace{[H_s^T\cdot Q_s\cdot F_s+H_{us}^T\cdot R_s\cdot F_{us}+H_{sn}^T\cdot P\cdot F_{sn}]}}{\underrightarrow{x}}_{ss|k}+2{\underrightarrow{c}}_k^T\underset{S_{us}}{\underbrace{H_{us}^T\cdot R_s\cdot T_{us}}}\cdot {\underrightarrow{u}}_{ss|k-1}+{\alpha }_t\end{array}$

where α_t is a term unrelated to ${\underrightarrow{c}}_k$.

To date, the original infinite time domain constraint optimization problem (20) has been transformed into an equivalent quadratic programming constraint optimization problem (35) by combining the plant (22), the maximum controlled allowable set (28), and the performance index $J_k^c$ (34).

(35)$\begin{array}{l}\underset{{\underrightarrow{c}}_k}{\underbrace{\min }}{J}_k^c={\underrightarrow{c}}_k^T\cdot S_{cs}\cdot {\underrightarrow{c}}_k+2{\underrightarrow{c}}_k^T(S_{xs}{x}_k+S_{ss}{\underrightarrow{x}}_{ss|k}+S_{us}{\underrightarrow{u}}_{ss|k-1})\\ s.t. \left\{\begin{array}{ccc}\hfill {\overline{x}}_{k+1}& =& \mathrm{\Phi }{\overline{x}}_k+B{c}_k\hfill \\ \hfill {\overline{y}}_k& =& C {\overline{x}}_k\hfill \end{array}\right.\\ N_s{\underrightarrow{c}}_k\le t_s-M_s {\overline{x}}_k- V_s {\underrightarrow{x}}_{ss|k}\end{array}$

The constraint optimization problem (35) can be solved by a quadratic programming solver, the resulting solution is denoted as ${\underrightarrow{c}}_k={\left[c_k,c_{k+1},\dots ,c_{k+n_{\mathrm{c}}-1}\right]}^{\mathrm{T}}$. Taking the first element c_k in the dual-mode control law (21) gives the control input

(36)$u_k=-K(x_k-x_{\mathrm{ss}|k})+u_{\mathrm{ss}|k-1}+c_k$

The scheme of the proposed C-ODPRC can be illustrated in Figure 2.

At each sampling instant k, the current and future dynamics of the lumped disturbance, denoted as ${\widehat{d}}_k$ and ${\underrightarrow{\widehat{d}}}_{k+1}$, are calculated via ESO (8) and LS-SVM based disturbance predictor (9)–(16). Then, the obtained disturbance series ${\underrightarrow{\widehat{d}}}_k={[{\widehat{d}}_k, {\underrightarrow{\widehat{d}}}_{k+1}]}^{\mathrm{T}}$ are sent to the steady-state target sequence solver (19) to obtain the (${\underrightarrow{x}}_{\mathrm{ss}|k}$, ${\underrightarrow{u}}_{\mathrm{ss}|k-1}$), followed by the construction of the maximum controlled allowable set according to steps ①–④. On this basis, the original optimization problem (20) is converted into an equivalent problem (35), and the control degree of freedom, ${\underrightarrow{c}}_k$, can be obtained by solving (35). By substituting the first element c_k into the dual-mode control law (21), we can obtain the feasible solution u_k. The above steps are repeated at the (k + 1)th sampling instant.

4. Results and Analysis

The control performance of the proposed C-ODPRC on the PTSF model (1) is demonstrated with two simulation cases. The control parameters are T_s = 1, n_c = 5, n_a = 5, Q = 0.1C^TC, R = 80, δ² = 10, G = 100, N = 45, n_p = 4.

Moreover, two other control methods are given as a comparison:

• Optimal disturbance predictive and rejection controller (ODPRC) proposed in [19], in which the future dynamics are considered in controller design without constraint handling mechanism.

• A special case of the C-ODPRC without considering the disturbance prediction (C-ODRC), namely n_a = 0.

All the above controllers are carried out in MATLAB2015b installed on a personal computer with a CPU of 2.5 GHz and RAM of 4.00 GB. In addition, the real running time is obtained by the ‘clock’ and ‘etime’ functions in MATLAB.

4.1. Simulation Research and Analysis Under Constant Operating Conditions

In this case, the field disturbances include the parabolic direct normal irradiation (I_n) of the sudden drop with an amplitude of −15 W/m² to imitate a passing cloud (3200 s–3700 s) and the fluctuation of the heat transfer fluid (HTF) with the form of d_o = 0.5cos(0.008t) (700 s–2700 s) to test the cosine form disturbance rejection ability. Both are drawn in the upper part of Figure 3a. Additionally, the constraints imposed on the heat transfer fluid are 6.5 kg/s ≤ q_H ≤ 8 kg/s.

The simulation is based on the typical condition #2 in Table 2. The reconstructed lumped disturbance and its estimate are shown in the lower part of Figure 3a. Figure 3 also exhibits the response curves of the field outlet temperature (b), the heat transfer fluid flow rate (c), and the control degree of freedom ${\underrightarrow{c}}_k$ (d) generated by the ODPRC, C-ODRC and C-ODPRC. The running time of all three control methods is given in Table 3.

Figure 3a shows that the ESO (8) can track the lumped disturbance with high estimate prediction accuracy. Moreover, it can be observable from Figure 3b that the overshoot corresponding to the field outlet temperature (T_H,out) generated by the ODPRC is the smallest (0.016%). In contrast, both the C-ODRC and the proposed C-ODPRC own larger overshoots of the field outlet temperature (the overshoot of both is 0.031% and 0.024%, respectively). However, the flow rate of heat transfer fluid (q_H) of the ODPRC violates the range of constraints, while the q_H of the other two controllers is within the limits of the constraints, implying that the ODPRC does not meet the control requirements.

Further comparisons are made between C-ODRC and C-ODPRC. It can be seen that the flow rate of heat transfer fluid of the two controllers is kept within the constraints despite the presence of multiple disturbances. Nevertheless, there are large differences in the flow curves shown in Figure 3c. In general, the heat transfer fluid flow rate of the C-ODPRC presents smoother changes, more timely control action, and a larger safety margin between the maximum fluid flow value and its upper limit than that of C-ODRC. The control performance of the two control methods is also compared. It is noticeable from Figure 3b that C-ODRC exhibits unsatisfactory dynamics with large oscillations and overshoot (0.031%), whereas the proposed C-ODPRC has better performance with smaller overshoot (0.024%). Finally, the actual running time of C-ODPRC at each sampling interval is slightly higher than that of C-ODRC (about 0.10 s and 0.12 s, respectively) because C-ODPRC contains an online prediction process for lumped disturbance. It is worth noting that the operating speed of C-ODPRC is also acceptable.

In summary, the proposed C-ODPRC achieves safer and more reasonable control curves due to the feedforward compensator of the future dynamics of the lumped disturbance. As a result, the goal of disturbance rejection can be better achieved with constraint satisfaction.

4.2. Simulation Research and Analysis Under Setpoint Tracking Conditions

This case aims to verify the control performance of each controller under setpoint tracking conditions in the presence of filed disturbances. The setpoint changes are set as follows:

• The PTSF operates at typical operating condition #2 from the beginning of the simulation to the 3000th second.

• The PTSF undergoes a step change from typical operating condition #2 to #1 at the 3000th second.

• The PTSF operates at typical condition #1 from the 3000th second to the 7000th second.

• The PTSF works from typical operating condition #1 to #3 (from the 7000th second to the 8000th second).

• The PTSF operates at typical condition #1 from the 8000th second to the 10,000th second.

The constraint on the heat transfer fluid flow rate is 6.5 kg/s ≤ q_H ≤ 8 kg/s. Field disturbances of the PTSF (see Figure 4a) contain parabolic direct normal irradiation (I_n), the heat transfer fluid (HTF) fluctuation with the form of d_o = 0.5cos(0.008t) (5500 s–6000 s), and the parameter perturbation d_in(j) = 0.5h_H·P_a,i[T_a(j) − T_H(j)]. The simulation results are shown in Figure 4b–d. The running time of all three control methods is exhibited in Table 3.

It is noticeable from Figure 4b,c that the proposed C-ODPRC can track the step and ramp changes in the setpoint in the presence of multiple disturbances. In addition, the overshoot and operating time of the outlet temperature curve generated by the C-ODPRC are higher than those of the ODPRC. This is mainly due to the idea of constraint satisfaction prior to disturbance rejection/setpoint tracking; thus, the heat transfer fluid is within the upper and lower limits at the cost of sacrificing certain disturbance rejection and setpoint tracking capabilities.

In contrast, the overshoot of the outlet temperature produced curve by C-ODPRC is comparable to that of C-ODRC, while the flow action of C-ODPRC is smoother, timelier, and with less oscillation. This benefits from the accurate prediction and feedforward compensation of the future dynamics of the lumped disturbance. Additionally, ODPRC owns the best tracking performance of the field outlet temperature with the smallest overshoot, the fastest running speed and the least oscillation, but at the price of constraint violation of the flow rate of the heat transfer fluid.

These all suggest that our proposed C-ODPRC possesses the best disturbance rejection and setpoint tracking capabilities under the premise of constraint satisfaction.

5. Conclusions

A constraint optimal model-based disturbance predictive and rejection controller (C-ODPRC) with disturbance prediction has been proposed to suppress disturbances of the parabolic trough solar field under constraint limitation. More precisely, a steady-state target sequence correction strategy is designed to integrate the current and future dynamics of lumped disturbance. This is followed by the construction of a maximum controlled allowable set based on the steady-state target sequence to convert infinite number constraint inequalities into finite ones. Finally, the transformed quadratic programming optimization problem can be solved. Simulation results suggest that the control performance of the proposed controller has been improved with constraint satisfaction at the cost of some anti-disturbance ability compared to ODPRC and C-ODRC. This is attributed to the combined action of the future dynamics of lumped disturbance and constraint handling mechanism. The proposed controller can be used to solve disturbance rejection problems under the constraint limitation for a class of thermal objects.

In future work, the time-delay processing mechanism in dealing with time-delay disturbances such as inlet temperature based on the proposed C-ODPRC is expected to be studied to provide better disturbance rejection performance.

Author Contributions

Conceptualization, S.W. and X.G.; methodology, S.W.; software, S.W. and X.G.; validation, S.W. and Y.L.; investigation, S.W. and X.G.; writing—original draft preparation, X.G.; writing—review and editing, S.W.; supervision, Y.L.; funding acquisition, S.W. and Y.L. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The original contributions presented in this study are included in the article. Further inquiries can be directed to the corresponding author.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. The diagram of the PTSF (a); the sectional drawing of A-A (b).

Figure 2. The scheme of the proposed C-ODPRC.

Figure 3. Response curves of different controllers with field disturbances under constant conditions. (a) Field disturbances and lumped disturbance; (b) Field outlet temperature; (c) The flow rate of heat transfer fluid; (d) Control degree of freedom of C-ODPRC.

Figure 4. Response curves of different controllers with field disturbances under setpoint tracking conditions. (a) Field disturbances and lumped disturbance; (b) Field outlet temperature; (c) The flow rate of heat transfer fluid; (d) The control of degree of freedom of C-ODPRC.

Table 1

Parameters of the PTSF.

Parameter	Value	Parameter	Value	Parameter	Value
η 1	0.94	Pa ,i	0.066	ρH	1487
η 2	0.96	Pa ,o	0.070	cH	2454
η 3	0.94	Pg ,i	0.109	ρa	7850
η 4	0.93	Pg ,o	0.115	ca	460
η 5	0.94	w	5.77	ρg	2400
η 6	0.94	l	600	cg	840
η 7	0.95	εa	0.18	Tatm	25
η 8	0.95	εg	0.88	Tsky	10

Table 2

Typical operating conditions of PTSF.

	TH,out,r/°C	qH,r/(kg/s)	In,r/(W/m2)	Tin,r/°C	Tatm,r/°C	Tsky,r/°C
#1	390.5	6.5	785	275	25	10
#2	390	7	851	275	25	10
#3	388	8	947	275	25	10

Note: the subscript r means the value at the typical operating condition.

Table 3

The real running time of the different control methods at each sampling interval (s).

Methods	ODPRC	C-ODRC	C-ODPRC
Case 1	0.104	0.100	0.120
Case 2	0.092	0.098	0.103

References

1. Pourasl, H.H.; Barenji, R.V.; Khojastehnezhad, V.M. Solar energy status in the world: A com-prehensive review. Energy Rep.; 2023; 10, pp. 3474-3493. [DOI: https://dx.doi.org/10.1016/j.egyr.2023.10.022]

2. Ahmad, N.; Liang, G.X.; Fan, P.; Zhou, H.Q. Anode interfacial modification for non-fullerene polymer solar cells: Recent advances and prospects. InfoMat; 2022; 4, e12370. [DOI: https://dx.doi.org/10.1002/inf2.12370]

3. Ferruzza, D.; Topel, M.; Laumert, B.; Haglind, F. Impact of steam generator start-up limitations on the performance of a parabolic trough solar power plant. Sol. Energy; 2018; 169, pp. 255-263. [DOI: https://dx.doi.org/10.1016/j.solener.2018.05.010]

4. Sanchez, A.J.; Gallego, A.J.; Escano, J.M.; Camacho, E.F. Adaptive incremental state space MPC for collector defocusing of a parabolic trough plant. Sol. Energy; 2019; 184, pp. 105-114. [DOI: https://dx.doi.org/10.1016/j.solener.2019.03.094]

5. Camacho, E.F.; Gallego, A.J.; Sánchez, A.J.; Berenguel, M. Incremental state-space model predictive control of a fresnel solar collector field. Energies; 2019; 12, 3. [DOI: https://dx.doi.org/10.3390/en12010003]

6. Sánchez, A.J.; Gallego, A.J.; Escaño, J.M.; Camacho, E.F. Event-based MPC for defocusing and power production of a parabolic trough plant under power limitation. Sol. Energy; 2018; 174, pp. 570-581. [DOI: https://dx.doi.org/10.1016/j.solener.2018.09.044]

7. Camacho, E.F.; Berenguel, M.; Gallego, A.J. Control of thermal solar energy plants. J. Process Contr.; 2014; 24, pp. 332-340. [DOI: https://dx.doi.org/10.1016/j.jprocont.2013.09.026]

8. de la Parte, M.P.; Cirre, C.M.; Camacho, E.F.; Berenguel, M. Application of predictive sliding mode controllers to a solar plant. IEEE Trans. Control Syst. Technol.; 2008; 16, pp. 819-825. [DOI: https://dx.doi.org/10.1109/TCST.2007.916298]

9. Camacho, E.F.; Gallego, A.J. Model predictive control in solar trough plants: A review. IFAC-PapersOnLine; 2015; 48, pp. 278-285. [DOI: https://dx.doi.org/10.1016/j.ifacol.2015.11.296]

10. Andrade, G.A.; Pagano, D.J.; Alvarez, J.D.; Berenguel, M. A practical NMPC with robustness of stability applied to distributed solar power plants. Sol. Energy; 2013; 92, pp. 106-122. [DOI: https://dx.doi.org/10.1016/j.solener.2013.02.013]

11. Reviriego, A.N.; Hernández-del-Olmo, F.; Alvarez-Barcia, L. Nonlinear adaptive control of heat transfer fluid temperature in a parabolic trough solar power plant. Energies; 2017; 10, 1155. [DOI: https://dx.doi.org/10.3390/en10081155]

12. Gallego, A.J.; Merello, G.M.; Berenguel, M.; Camacho, E.F. Gain-scheduling model predictive control of a Fresnel collector field. Control Eng. Pract.; 2019; 82, pp. 1-13. [DOI: https://dx.doi.org/10.1016/j.conengprac.2018.09.022]

13. Himour, Y.; Tadjine, M.; Boucherit, M.S. Nonlinear and infinite gain scheduling neural predictive control of the outlet temperature in a parabolic trough solar field: A comparative study. Eng. Appl. Artif. Intell.; 2023; 126, 106862. [DOI: https://dx.doi.org/10.1016/j.engappai.2023.106862]

14. Chanfreut, P.; Maestre, J.M.; Gallego, A.; Annaswamy, A.M.; Camacho, E.F. Clustering-based model predictive control of solar parabolic trough plants. Renew. Energy; 2023; 216, 118978. [DOI: https://dx.doi.org/10.1016/j.renene.2023.118978]

15. Rossiter, J.A. Model-Based Predictive Control: A Practical Approach; 1st ed. CRC Press: Boca Raton, FL, USA, 2003.

16. Alsharkawi, A.; Rossiter, J.A. Gain scheduling dual mode MPC for a solar thermal power plant. IFAC-PapersOnLine; 2016; 49, pp. 128-133. [DOI: https://dx.doi.org/10.1016/j.ifacol.2016.10.151]

17. Alsharkawi, A.; Rossiter, J.A. Towards an improved gain scheduling predictive control strategy for a solar thermal power plant. IET Control Theory Appl.; 2017; 11, pp. 1938-1947. [DOI: https://dx.doi.org/10.1049/iet-cta.2016.1319]

18. Gao, X.H.; Wei, S.S.; Su, Z.G. Optimal model predictive rejection control for nonlinear parabolic trough collector with lumped disturbances. Trans. Inst. Meas. Control; 2021; 43, pp. 1903-1914. [DOI: https://dx.doi.org/10.1177/0142331220983651]

19. Gao, X.H.; Wei, S.S.; Wang, M.L.; Su, Z.G. Optimal disturbance predictive and rejection control of a parabolic trough solar field. Int. J. Robust Nonlinear Control; 2022; early view [DOI: https://dx.doi.org/10.1002/rnc.6443]

20. Wei, S.S.; Liang, X.F.; Mohsin, T.; Wu, X.; Li, Y.G. A simplified dynamic model of integrated parabolic trough concentrating solar power plants: Modeling and validation. Appl. Therm. Eng.; 2020; 169, 114982. [DOI: https://dx.doi.org/10.1016/j.applthermaleng.2020.114982]

21. Pechlivanidou, G.; Karampetakis, N. Zero-order hold discretization of general state space systems with input delay. IMA J. Math. Control Inf.; 2022; 39, pp. 708-730. [DOI: https://dx.doi.org/10.1093/imamci/dnac005]

22. Su, Z.G.; Sun, Y.G.; Zhu, X.J.; Chen, Z.; Sun, L. Robust tuning of active disturbance rejection controller for time-delay systems with application to a factual electrostatic precipitator. IEEE Trans. Control Syst. Technol.; 2022; 30, pp. 2204-2211. [DOI: https://dx.doi.org/10.1109/TCST.2021.3127794]

23. Sun, L.; Dong, J.Y.; Li, D.H.; Xue, Y.L. Cascaded load-frequency digital electric-hydraulic regulation of turbine via extended state observer. Proc. CSEE; 2015; 35, pp. 1697-1703.

24. Zheng, Y.L.; He, F.; Shen, X.Z.; Jiang, X.S. Energy Control Strategy of Fuel Cell Hybrid Electric Vehicle Based on Working Conditions Identification by Least Square Support Vector Machine. Energies; 2020; 13, 426. [DOI: https://dx.doi.org/10.3390/en13020426]

25. Hossain, S.K.S.; Ayodele, B.V.; Almithn, A. Predictive modeling of bioenergy production from fountain grass using gaussian process regression: Effect of kernel functions. Energies; 2022; 15, 5570. [DOI: https://dx.doi.org/10.3390/en15155570]