Dynamic Optimization of Variable Load Process for Combined Heat and Power Unit Based on Sequential Quadratic Programming and Interior Point Method Alternating Solution Method

Abstract

Aiming at the problem that the modeling and solving method of combined heat and power (CHP) unit variable load control process is challenging to meet the demand for efficient analysis of complex systems, this paper proposes a method based on sequential quadratic programming and interior point method (SQP-IPM) alternating solution for dynamic optimization of the CHP unit variable load process. Firstly, by constructing the CHP unit mechanism model, multi-variable coordination control constraints, and output variable process constraints, the dynamic optimization proposition of the CHP unit variable load control process is formed. Then, the large-scale nonlinear programming (NLP) problem formed by using the orthogonal configuration method to discrete the state and control variables is optimally solved using the IPM-SQP alternating solution method. Further, from the perspective of balancing the accuracy of the solution and computational efficiency, the flexible convergence depth control (CDC) strategy is introduced into the alternative solution method to improve the real-time performance of the algorithm. Finally, the variable load control process of 300MW extraction CHP unit is simulated to verify that the proposed method reduces the calculation time for 12 consecutive variable load scenarios by about 70%, effectively improving the real-time performance of scenario applications.

1. Introduction

In order to meet the demand for variable load capability of CHP units participating in integrated energy system frequency regulation [1], short-term dispatch, and other scenarios, optimal control technology modification solutions to improve the fast regulation capability of units are constantly proposed [2,3]. The dynamic regulation characteristics of the modified units are complex, showing characteristics such as multivariable, multi-process coupling, multiparameter superposition, and strong nonlinearity, which often lead to difficulties in accurate modeling and fast simulation calculation of the variable load response capability of the units. Differential algebraic equations (DAEs) are often used to describe the dynamic characteristics of the unit in the operating interval, which makes NLP with the constraints of the DAEs model of CHP units a dynamic optimization problem. Therefore, how to efficiently and reliably solve the optimal control problem of the modified CHP units becomes an urgent problem for integrated energy utilization.

The widely used simultaneous method [4,5,6] for dynamic process control optimization is often used for dynamic optimization problem solving, and the method has advantages in dealing with complex constrained CHP unit variable load control optimization problems. However, its practice of discretizing both state and control variables often makes the size of the discretized NLP proposition exponentially larger. The discrete operation will lead to more difficult optimization calculations and increased solution time, which often makes it difficult to meet the computational performance requirements of optimization solutions for different application scenarios of integrated energy systems. The existing research on NLP problem solving is mainly in the following areas. Full-space methods and simplex space methods are often used to deal with problems with equation constraints. For large-scale optimization problems with a large number of equation constraints and limited degrees of freedom, the simplex space method can be used to reduce the dimensionality of the SQP subproblem through a spatial decomposition strategy, thereby reducing the computational effort of the Hessian matrix and improving the solution efficiency of the algorithm [7,8,9]. However, the spatial decomposition process adds additional computational cost and its convergence is not as good as the full-space method based on exact second-order derivative information. The positive set method and the IPM are commonly used for problems with inequality constraints. IPM converts inequality constraints into obstacle terms added to the objective function, and the original problem solution is converted into a series of equationally constrained obstacle subproblem solutions, thus avoiding the determination of the positive set by the solution process [10,11]. Nevertheless, algorithms that use linear search strategies based on evaluation functions and trust strategies [12,13,14] to ensure global convergence, although less dependent on the initial point, cannot balance convergence and convergence speed when solving large-scale complex process optimization problems where the initial point is far from the optimal point. The algorithm has better local convergence and can converge quickly only when the initial point is closer to the optimal point; otherwise, it may fail to converge.

Meanwhile, another reason limiting the convergence of the above algorithm is the waste of computational resources caused by the convergence criterion. The convergence criterion significantly impacts the computational time, and the computational delay will lead to a decrease in controller performance. The traditional Himmelblau termination criterion (H-criterion for short) [15] is commonly used to control all or part of the convergence of optimization algorithms. The H-criterion is rigid and can only give a simple conclusion of “convergence” or “non-convergence”. When considering the degree of change of variables and the degree of reduction (or increase) of the objective function, it often makes the optimization calculation index enter a stage where the computational cost is much greater than the improvement of the convergence effect (“over convergence”) after reaching a certain level of convergence. Nonlinear optimization solution algorithms based on traditional criteria often need to satisfy pre-given tolerance values when solving dynamic optimization problems. Without significantly improving the numerical solution, over-iterations consume much computational time, making the control action not act on the controlled system in time. Moreover, there is no general termination criterion applicable to all optimization problems and various optimization algorithms in existing studies.

To address the above problems, this paper proposes the dynamic optimization of the variable load process of the CHP unit based on the IPM-SQP alternating solution method. The paper is structured as follows: Section 2 constructs the CHP unit mechanism model, designs the multivariable coordinated control model and the control process output variable constraints, and forms the CHP unit variable load control process dynamic optimization proposition. In Section 3, the finite element orthogonal collocation method is used for discretization to generate large-scale NLP problems. In Section 4, an alternating IPM-SQP solution strategy based on CDC is designed to solve the NLP problems to improve the algorithm’s real-time performance. Section 5 describes a case study to prove the effectiveness of the proposed method. Section 6 presents our conclusions.

2. Proposition Construction for Dynamic Optimization of CHP Unit Variable Load Process

CHP unit variable load control is the process of regulating the operating unit’s electrical output and thermal output, mainly in response to changes in the electrical and thermal energy demand of the energy supply system. As shown in Figure 1, the variable load process controls the amount of fuel entering the boiler to produce high-temperature steam by adjusting the coal feed mass flow rate V_B through the controller. The amount of steam entering the high-pressure cylinder is controlled by adjusting the turbine regulating valve opening V_T. The amount of steam entering the low-pressure cylinder to do work is controlled by adjusting the opening of the pumping regulating butterfly valve V_H to output electrical power, while the amount of steam in the medium-pressure cylinder is controlled to supply the heat source. DAEs describe the CHP unit mechanism relationship. The unit control variables are constrained by the control method to meet the variable load demand. To ensure safe and stable operation, the unit control process is subject to state and output constraints. This process control optimization seeks the feasible regulation path of the unit electric output and thermal output change process from initial to final value under satisfying the performance requirements. Therefore, the CHP unit variable load problem can be described as an optimal control problem.

2.1. Objective Function of the Optimization Proposition

To ensure that the control process of the CHP unit is accurate and fast, and smooth to achieve variable load regulation, the mathematical form of the combination of control variables and output variables is chosen as the optimal control performance objective function [16], which is described as

$$\underset{u,y}{min}J={{\displaystyle \int }}_{t_0}^{t_f}{\left(\parallel \mathit{u}\left(t\right)-\mathit{u}\left(t-1\right){\parallel }_R^2+\parallel \mathit{y}\left(t\right)-{\mathit{y}}^{\ast }\parallel {}_S^2\right)}^{2}dt$$

(1)

where ${\mathit{y}}^{\ast }$ are the set values of output variables; $\mathit{u}\left(t\right)$ is the control variable at the t-th moment, $\parallel \mathit{u}\left(t\right)-\mathit{u}\left(t-1\right){\parallel }_{\mathit{R}}^2$ minimizes the smooth variation of the control variable; $\mathit{y}\left(t\right)$ is the measured output of the control object at the t−1-th moment as the feedback of the closed-loop system at the t-th moment; R and S are positive definite weighting matrix to determine the initial values by normalization.

In the model, $\mathit{u}$ is the control variable, $\mathit{y}$ is the algebraic state variable, and $\mathit{x}$ is the differential state variable, where:

$$\left\{\begin{array}{c}\mathit{x}={\left[{\psi }_d,P_H,{\psi }_z,V_f\right]}^T\\ \mathit{u}={\left[V_T,V_B,V_H\right]}^T\\ \mathit{y}={\left[{\psi }_t,P_H,m_H\right]}^T\end{array}\right.$$

(2)

2.2. CHP Unit Model Constraints

2.2.1. Mechanism Model

By analyzing the thermal component characteristics of the heating unit and the steam work process in the turbine [17], a simplified nonlinear mechanism model description of the CHP unit can be obtained as follows:

$$V_m\left(t\right)=V_B\left(t-t_B\right)$$

(3)

$$\frac{d{V}_f\left(t\right)}{dt}=-\frac{V_f\left(t\right)}{T_f}+\frac{V_B\left(t-t_B\right)}{T_f}$$

(4)

$$\frac{d{\psi }_d\left(t\right)}{dt}=-\frac{K_3{\psi }_t\left(t\right)V_T\left(t\right)}{C_d}+\frac{K_1{V}_f\left(t\right)}{C_d}$$

(5)

$${\psi }_d\left(t\right)-{\psi }_t\left(t\right)=K_2{\left(K_1{V}_f\left(t\right)\right)}^{1.5}$$

(6)

$$\frac{d{P}_H\left(t\right)}{dt}=\frac{K_3{K}_4{\psi }_t\left(t\right)V_T\left(t\right)}{T_t}+\frac{K_5{\psi }_z\left(t\right)V_H\left(t\right)}{T_t}-\frac{P_H\left(t\right)}{T_t}$$

(7)

$${\psi }_1\left(t\right)=0.01{\psi }_t\left(t\right)V_T\left(t\right)$$

(8)

$$\frac{d{\psi }_z\left(t\right)}{dt}=-\frac{K_6{m}_r\left(t\right)\left(96{\psi }_z\left(t\right)-{\epsilon }_r\left(t\right)+103\right)}{C_z}+\frac{K_3\left(1-K_4\right){\psi }_t\left(t\right)V_T\left(t\right)}{C_z}-\frac{K_5{\psi }_z\left(t\right)V_H\left(t\right)}{C_z}$$

(9)

$$m_H\left(t\right)=K_7{K}_6{m}_r\left(t\right)\left(96{\psi }_z\left(t\right)-{\epsilon }_r\left(t\right)+103\right)$$

(10)

Equations (3) and (4) characterize the dynamic relationship of the coal-fired pulverizing system. Equations (5) and (6) characterize the dynamic relationship of the steam package boiler system. Equations (7)–(9) characterize the dynamic relationship of the turbine system. Equation (10) characterizes the dynamic relationship of the heating system. K₁, K₂, K₃, K₄, K₅, K₆, and K₇ are static parameters.

2.2.2. Multivariate Coordinated Control Model

The CHP unit controller cooperatively regulates the direction and magnitude of $V_B\left(t\right)$, $V_T\left(t\right)$, and $V_H\left(t\right)$ through the designed control algorithm, which makes the output variables ${\psi }_t\left(t\right)$, $m_H\left(t\right)$, and $P_H\left(t\right)$ change. For the above mechanism model, a coordinated thermal-electric control algorithm is designed to achieve the three critical functions of accurate electric power tracking, fast thermal power recovery, and system pressure safety [2,3]. Its mathematical model is as follows:

$$\left\{\begin{array}{c}{V}_T\left(t\right)=K_{PT}{E}_T\left(t\right)+K_{IT}{{\displaystyle \int }}_0^t{E}_T\left(t\right)dt+K_{DT}\frac{d{E}_T\left(t\right)}{dt}\\ V_B\left(t\right)=K_{PB}{E}_B\left(t\right)+K_{IB}{{\displaystyle \int }}_0^t{E}_B\left(t\right)dt+K_{DB}\frac{d{E}_B\left(t\right)}{dt}\\ V_H\left(t\right)=K_{PH}{E}_H\left(t\right)+K_{IH}{{\displaystyle \int }}_0^t{E}_H\left(t\right)dt+K_{DH}\frac{d{E}_H\left(t\right)}{dt}\\ E_T\left(t\right)=S{P}_{\psi }-{\psi }_t\left(t\right)\\ E_B\left(t\right)=S{P}_P-\left(m_H\left(t\right)-\frac{K{e}^{-\frac{t}{T_c}}}{T_c}{P}_H\left(t\right)\right)\\ E_H\left(t\right)=S{P}_m-m_H\left(t\right)\end{array}\right.$$

(11)

where K_PT, K_IT, K_DT, K_PB, K_IB, K_DB, K_PH, K_IH, K_DH, K, and Tc are the control system regulation parameters.

2.2.3. Output Variable Constraints for Control Processes

Considering the safety and stability of CHP unit operation, the fluctuation range of the process value of the output variable and the error range of the steady-state value during the control process should not exceed the allowed limits, and the constraints should be satisfied as follows:

$$\left\{\begin{array}{c}-{\mathsf{\Delta }}_{\psi }\le S{P}_{\psi }-{\psi }_t\left(t\right)\le {\mathsf{\Delta }}_{\psi }\\ -{\mathsf{\Delta }}_P\le S{P}_P-P_H\left(t\right)\le {\mathsf{\Delta }}_P\\ -{\mathsf{\Delta }}_m\le S{P}_m-m_H\left(t\right)\le {\mathsf{\Delta }}_m\end{array}\right.$$

(12)

$$\left\{\begin{array}{c}-{\delta }_{\psi }\le {\psi }_t\left(t_f\right)-S{P}_{\psi }\le {\delta }_{\psi }\\ -{\delta }_P\le P_H\left(t_f\right)-S{P}_P\le {\delta }_P\\ -{\delta }_m\le m_H\left(t_f\right)-S{P}_m\le {\delta }_m\end{array}\right.$$

(13)

2.3. Boundary Constraints

The CHP unit completes variable load regulation in the time period $\left[t_0,t_f\right]$, and the initial and termination states of each control variable and output variable are shown below:

$$\left\{\begin{array}{c}{\psi }_t\left(t_0\right)={\psi }_t^0\\ V_T\left(t_0\right)=V_T^0\\ P_H\left(t_0\right)=P_H^0\\ V_B\left(t_0\right)=V_B^0\\ m_H\left(t_0\right)=m_H^0\\ V_H\left(t_0\right)=V_H^0\end{array}\right.$$

(14)

$$\left\{\begin{array}{c}{\psi }_t\left(t_f\right)={\psi }_t^f\\ V_T\left(t_f\right)=V_T^f\\ P_H\left(t_f\right)=P_H^f\\ V_B\left(t_f\right)=V_B^f\\ m_H\left(t_f\right)=m_H^f\\ V_H\left(t_f\right)=V_H^f\end{array}\right.$$

(15)

where $\left({\psi }_t^0,P_H^0,m_H^0\right)$ and $\left({\psi }_t^f,P_H^f,m_H^f\right)$ are the output values at the start and end times, respectively. $\left(V_T^0,V_B^0,V_H^0\right)$ and $\left(V_T^f,V_B^f,V_H^f\right)$ are the control values at the start and end times, respectively.

2.4. Path Constraints

Due to the limits of the CHP unit’s own component properties, V_T, V_B, and V_H need to meet the following constraints during variable load regulation.

$$\left\{\begin{array}{c}{V}_T^{min}\le V_T\left(t\right)\le V_T^{max}\\ V_B^{min}\le V_B\left(t\right)\le V_B^{max}\\ V_H^{min}\le V_H\left(t\right)\le V_H^{max}\end{array}\right.$$

(16)

where $V_T^{min}$ and $V_T^{max}$ are the minimum and maximum values of the variation of the turbine regulator opening, respectively; V_B and V_H are the same.

3. Discretization of Dynamic Optimization Problems

Without loss of generality, the unit variable load optimal control problem obtained in Chapter 1 can be expressed as a dynamic optimization problem in the form of DAEs as follows:

$$\begin{gathered} \min J\left(\mathit{x}\left(t_f\right)\right) \\ \mathrm{s.t.}\left\{\begin{array}{c}\dot{\mathit{x}}\left(t\right)=F\left(\mathit{x}\left(t\right),\mathit{y}\left(t\right),\mathit{u}\left(t\right)\right),\\ g_E\left(\mathit{x}\left(t\right),\mathit{y}\left(t\right),\mathit{u}\left(t\right)\right)=0,\\ h_E\left(\mathit{x}\left(t_f\right),\mathit{y}\left(t_f\right),\mathit{u}\left(t_f\right)\right)=0,\\ g_{IE}\left(\mathit{x}\left(t\right),\mathit{y}\left(t\right),\mathit{u}\left(t\right)\right)<0,\\ h_{IE}\left(\mathit{x}\left(t_f\right),\mathit{y}\left(t_f\right),\mathit{u}\left(t_f\right)\right)<0,\\ {\mathit{u}}_L\le \mathit{u}\left(t\right)\le {\mathit{u}}_U,\\ \mathit{x}\left(t_0\right)={\mathit{x}}^0,\mathit{x}\left(t_f\right)={\mathit{x}}^f,\\ t\in \left[t_0,t_f\right]\end{array}\right. \end{gathered}$$

(17)

where F is the dynamic model of the CHP unit DAEs, corresponding to Equations (3)–(10); $g_E$ and $g_{IE}$ are the equal and unequal path constraints, corresponding to Equations (11), (12) and (16); h_E and h_IE are the final value constraints at moment $t_f$, corresponding to Equation (13); the initial and termination boundary values of the state variable $\mathit{x}$ are calculated by Equations (14) and (15).

The whole time domain [$t_0,t_f$] is divided into 𝑁-segment finite elements, and on each finite element $\left[t_{i-1},t_i\right]\left(i=1,2,\dots ,N\right)$, the Lagrange interpolation function based on Radau orthogonal configuration points approximates the state variables, control variables and algebraic output variables of the original functions. The interpolation equation is as follows:

$$\left\{\begin{array}{c}x\left(t\right)={\displaystyle {\displaystyle \sum }_{j=0}^K}{l}_j\left(\tau \right)x_{ij}\\ u\left(t\right)={\displaystyle {\displaystyle \sum }_{j=1}^K}{\bar{l}}_j\left(\tau \right)u_{ij}\\ y\left(t\right)={\displaystyle {\displaystyle \sum }_{j=1}^K}{\bar{l}}_j\left(\tau \right)y_{ij}\end{array}\right.$$

(18)

where $l_j\left(\tau \right)={\prod }_{k=0,\ne j}^K\left[\left(\tau -{\tau }_k\right)/\left({\tau }_j-{\tau }_k\right)\right]$ and ${\bar{l}}_j\left(\tau \right)={\prod }_{k=1,\ne j}^K\left[\left(\tau -{\tau }_k\right)/\left({\tau }_j-{\tau }_k\right)\right]$ are the Lagrangian interpolation polynomials of the corresponding variables, respectively. K is the interpolation order.

The initial and final value conditions of the state variables are:

$$x_{1,0}=x^0, x_f=x_{N,K}$$

(19)

Since the state variables are derivable, the values of the state variables at the nodes of adjacent finite element connections should be continuous, so there is a continuity condition as:

$$x_{i+1,0}={\displaystyle {\displaystyle \sum }_{j=0}^K}{l}_j\left(1\right)x_{ij},i=1,2,\cdots ,N-1$$

(20)

The advantage of the Lagrange interpolation polynomial is that the value of the variable at each configuration point is exactly equal to its coefficient.

$$t_{ij}=t_{i-1}+\left(t_i-t_{i-1}\right){\tau }_j$$

(21)

$$x\left(t_{ij}\right)=x_{ij}, u\left(t_{ij}\right)=u_{ij}, y\left(t_{ij}\right)=y_{ij}$$

(22)

Substituting Equations (18), (21), and (22) into the set of differential equations of proposition (17), we obtain the configuration equation:

$$R_i\left({\tau }_j\right)={\displaystyle {\displaystyle \sum }_{k=0}^K}{\dot{l}}_k\left({\tau }_j\right)x_{ik}-h_{i}f\left(x_{ij},u_{ij}\right)=0,i=1,2,\cdots ,N,j=1,2,\cdots ,K.$$

(23)

$$t=t_{i-1}+h_i\tau , h_i=t_i-t_{i-1}$$

(24)

So far, the NLP proposition form of the original optimal proposition (17) after discretization is obtained:

$$\begin{gathered} \min J\left(x\left(t_f\right)\right) \\ \mathrm{s.t.}\begin{array}{c}\left\{\begin{array}{c}{\displaystyle \sum _{k=0}^K}{\dot{l}}_k\left({\tau }_j\right)x_{ik}-h_{i}F\left(x_{ij},y_{ij},u_{ij}\right)=0,\\ g_E\left(x_{ij},y_{ij},u_{ij}\right)=0,\\ h_E\left(x_{ij},y_{ij},u_{ij}\right)=0,\\ g_{IE}\left(x_{ij},y_{ij},u_{ij}\right)<0,\\ h_{IE}\left(x_{ij},y_{ij},u_{ij}\right)<0,\\ u_L\le u_{ij}\le u_U,\\ i=1,2,\cdots ,N;j=1,2,\cdots ,K;\\ x\left(0\right)=x^0,x\left(t_f\right)=x^f\\ x_{i+1,0}={\displaystyle {\displaystyle \sum }_{j=0}^K}{l}_j\left(1\right)x_{ij}, i=1,2,\cdots ,N-1.\end{array}\right.\end{array} \end{gathered}$$

(25)

For the above NLP problems containing many inequality constraints after discretization, SQP and IPM [18] are often used for optimization solutions. However, they still do not provide a good balance of convergence and convergence speed when solving large-scale complex process optimization problems. The SQP algorithm with good local convergence often converges quickly only when the initial point is close to the optimal point; otherwise, it may fail to converge. The IPM algorithm with good global convergence does not depend strongly on the initial point, but the local convergence and convergence speed are difficult to guarantee. Currently, the solution of dynamic optimization problems faces the challenge that the convergence rate of the optimization search and the solution accuracy are difficult to be balanced at the same time. In the CHP unit nonlinear optimal control of variable load problems, the real-time requirement is very high, and the solution time is vital. The computational delay prevents the control action from being applied to the controlled object promptly, resulting in a degradation of the system control performance.

4. IPM-SQP Alternate Solution Based on Convergence Depth Control

In this paper, we design an alternate solution method using the advantages of each IPM and SQP algorithm. The two nonlinear planning algorithms are regarded as the elemental algorithms of the alternate solution method, and the CDC is used to realize the alternate solution of the two elemental algorithms in collaboration with each other. This method can improve the algorithm’s overall robustness and computational efficiency to complete the solution of large-scale complex problems.

4.1. Convergence Depth Control

Due to the requirements of the CHP unit variable load process on the real-time performance of the solution algorithm, the design of the dynamic optimization convergence criterion should take into account the computational accuracy and operational efficiency, and a balance should be achieved between them. The traditional convergence criterion has only convergence success and failure states when terminating the algorithm, which is too strict in practical applications [19] and is not conducive to real-time optimization calculation. In order to better control the termination condition of the algorithm, this paper introduces the convergence depth and progress degree to describe the convergence state at the current iteration point quantitatively. Furthermore, the Sigmoid function is used to soften the traditional rigid convergence criterion, to control the termination of the algorithm.

The nonlinear programming problem shown in the original proposition (25) can be abbreviated as:

$$\begin{array}{c}\underset{x\in R^{\mathrm{n}}}{\min } f\left(x\right) \\ \mathrm{s}.\mathrm{t}.\left\{\begin{array}{c}c\left(x\right)=0\\ x\ge 0\end{array}\right.\end{array}$$

(26)

where $f$ is a scalar objective function, $x\in R^n$ is an n-dimensional real vector, and $c\in R^m$ is an m-dimensional equation. The solution to the problem (26) must satisfy the following first-order optimality condition:

$$\left\{\begin{array}{c}\nabla f\left(x\right)+\nabla c\left(x\right)\lambda -v=0\\ c\left(x\right)=0\\ Xv=0,v\ge 0\end{array}\right.$$

(27)

where $\nabla f\left(x\right)$ is the gradient vector of the objective function; $\nabla c\left(x\right)$ is the Jacobi matrix of the constraint equation; 𝜆 is the Lagrangian multiplier corresponding to the constraint equation; $\nu$ is the pairwise variable corresponding to $x$, and $X=\mathrm{diag}\left(x\right)$.

Assuming that the sequence of iterations generated by the optimization algorithm is {x^k}, define the feasibility error ${\delta }_{\mathit{feasErr} }^k$ at the current iteration point x^k, and the amount of improvement predicted by the objective function ${\delta }_{objErr}^k$ [20]:

$${\delta }_{feasErr}^k=\max \left\{\parallel c\left(x^k\right)\parallel ;\max \left\{-x_{\left(i\right)}^k\right\}\right\}$$

(28)

$${\delta }_{objErr}^k=\left|\nabla f^T\left(x^k\right)d^k\right|,d^k=x^k-x^{k-1}$$

(29)

where $x_{\left(i\right)}^k$ is the i-th component of the vector.

S is introduced as a deformed Sigmoid function for softening the traditional rigid convergence criterion defined as:

$$S\left({\delta }^k,{\epsilon }_0\right)=\left[\tanh \left({\xi }^{\ast }\left(\log {\delta }^k/\log {\epsilon }_0\right)\right)\right]/\left[\tanh \left(\xi \right)\right]$$

(30)

where the function $S\left({\delta }^k,{\epsilon }_0\right)$ smoothly connects the interval $\left[{\epsilon }_0,1/{\epsilon }_0\right]$; $\xi$ measures the degree of smoothness of the function S in this interval, and $\xi$ takes 1.5.

Then the depth of convergence ${\theta }_{\mathit{convg}}^k$ is:

$${\theta }_{\mathit{convg}}^k=S\left(\max \left\{{\delta }_{\mathit{feasErr}}^k,{\delta }_{\mathit{objErr}}^k\right\},{\epsilon }_0\right)$$

(31)

The softening of the Sigmoid function to the rigid convergence criterion is shown in Figure 2. Given different tolerances ${\epsilon }_0$ (1 × 10⁻¹⁰ for curve 1, 1 × 10⁻⁵ for curve 2, 1 × 10⁻³ for curve 3, and rigid criterion for curve 4). The degree of deformation of curve 2 is more appropriate, while curve 1 is too flat and curve 3 is too steep.

The depth of convergence based on the flexible criterion is more intuitive than the first-order optimality condition (Equation (27)). The criterion uniformly characterizes the quality (feasibility and optimality) of each iteration point generated by the optimization algorithm and can quantitatively describe the distance between it and the optimal point [19,21]. When the depth of convergence reaches 1.0, the feasibility error and the amount of improvement in the objective function prediction also reach ɛ₀, at which point the iteration point is already a good approximation of the optimal point.

4.2. IPM-SQP Alternate Solution Algorithm

The CDC-based IPM-SQP alternate solution method uses the globally convergent IPM algorithm to advance the initial point to the optimal point and uses the depth of convergence to measure the distance between the current iteration point and the optimal point. By controlling the degree of advancement, and then switching to the SQP algorithm with better local convergence to make the iteration point converge to the optimal point quickly. This paper introduces the progress degree index for switching between IPM-SQP alternate solution methods.

The amount of feasible change ${\delta }_{\mathit{feasChg} }^k$ and the amount of change in the objective function ${\delta }_{\mathit{objChg}}^k$ are:

$${\delta }_{\mathit{feasChg}}^k=\left|{\delta }_{\mathit{feasErr}}^k-{\delta }_{\mathit{feasErr}}^{k-1}\right|$$

(32)

$${\delta }_{\mathit{objChg}}^k=\left|f\left(x^k\right)-f\left(x^{k-1}\right)\right|$$

(33)

From the above equation, the degree of progress ${\theta }_{\mathit{prog} }^k$ can be defined as [22]:

$${\theta }_{\mathit{prog} }^k=S\left(\max \left\{{\delta }_{\mathit{feasChg} }^k,{\delta }_{objChg}^k\right\},{\epsilon }_0\right)$$

(34)

When the degree of progress reaches 1.0, the algorithm has reached ${\epsilon }_0$ in terms of the amount of feasibility change and the amount of improvement in the objective function, which means that the algorithm has almost no more progress in terms of the objective function and feasibility improvement.

The CDC-IPM-SQP alternate solution method framework is designed based on the depth of convergence and the degree of progress introduced above. The framework for solving the dynamic optimization problem in the form of Equation (25) using the simultaneous method is shown in Figure 3. First, the finite element Orthogonal collocation method is used for discretization, and then the IPM-SQP optimization algorithm is used for optimization.

(1)

Establishing dynamic optimization problems

Design input for dynamic optimization problems based on the dynamic process scenario of CHP unit variable load optimization control, which includes: the CHP unit mechanism model, multiple control variable coordination control model, control process output variable constraints, boundary constraints, and path constraints.

(2)

Discrete processing of original problems

The variables of the original optimization problem are fully discretized by the finite element Orthogonal collocation method. Lagrange interpolation function based on Radau Orthogonal collocation points approximates the original functions of state variables and control variables. The original dynamic optimization problem after discretization becomes an NLP problem.

(3)

Real-time planning NLP model

Follow the process of inputting “initial information → generating initial values → optimizing the solver → obtaining the optimal solution” for real-time planning. During the optimization process, input NLP model information and calculated residuals into the solver. The variable values obtained from each optimization iteration are fed back to the NLP model. When using the SQP-IPM optimization algorithm for calculation and solution, it is necessary to select configuration points, calculate derivatives, approximate Jacobian matrices for the NLP model, and input the sparse Jacobian matrix into the optimization algorithm.

In the optimization algorithm solving section, IPM is first selected to solve the NLP problem, and the depth of convergence at each step is calculated using Equation (31) during the iteration. When the convergence depth satisfies the specified threshold ${\theta }_c$, the information of the iteration points is saved; Terminate the current algorithm and switch to SQP, and the last saved iteration point is used as the initial point to continue the iteration. Since different algorithms have different local convergence radii, this paper first gives the threshold ${\theta }_c$ an initial value of ${\theta }_c^0\in \left(0,1\right)$. If the SQP convergence is unsuccessful, the iteration point obtained from IPM enters the local convergence radius; the convergence depth threshold is increased by $d_{\theta }\in \left(0,1\right)$, and returns to IPM for iteration. The specific steps of the alternate solution method are as follows:

(1)

Let k = 0, given the initial point x⁰; specify the convergence depth threshold ${\theta }_c$, the progress degree threshold ${\theta }_p$, the convergence depth advance step $d_{\theta }$, and the maximum convergence depth value ${\theta }_{max}$; and set the flag bit IPM to 0.

(2)

Select the IPM and set its flag bit to 1 for solving the optimization problem.

(3)

Iterate one step, k = k + 1, if the termination criterion of IPM is satisfied, the solution is successful, and the whole algorithm is terminated; otherwise, go to step 4).

(4)

Calculate the depth of convergence, check whether ${\theta }_{conv\mathrm{g}}^k\ge {\theta }_c$ is satisfied; if it is satisfied, assign the current iteration point x^k to x⁰ and clear the flag bit of SQP, go to step 6); otherwise, go to step 5).

(5)

If the termination criterion (solution failure) of IPM is not satisfied, turn to step 3); otherwise, the solution fails, and the whole algorithm terminates.

(6)

Select SQP, set its flag bit to 1, and solve the optimization problem with x⁰ as the initial point.

(7)

Iterate one step, k = k + 1; if the termination criterion of SQP is satisfied, the solution succeeds, and the whole algorithm terminates; otherwise, go to step 8).

(8)

Calculate the degree of progress; if both ${\theta }_{conv\mathrm{g}}^k\ge {\theta }_p$ and the termination criterion of SQP are not satisfied, turn to step 7); otherwise, go to step 9).

(9)

Assign the current iteration point x^k to x⁰ and go to step 10).

(10)

Let ${\theta }_c={\theta }_c+d_{\theta }$; if ${\theta }_c\ge {\theta }_{max}$, the solution fails, and the whole algorithm terminates; otherwise, the flag bit of IPM is cleared to 0; turn to step 2).

The simultaneous method framework and the specific implementation flow of the CDC-based IPM-SQP alternate solution to dynamic optimization problems in this paper are shown in Figure 3.

Figure 3. Simultaneous method framework for alternating SQP-IPM solutions based on CDC.

Figure 3. Simultaneous method framework for alternating SQP-IPM solutions based on CDC.

5. Simulation Case and Analysis

The numerical experiments in this paper are conducted on a computer-based Intel(R) Core(TM) i7-8550U processor with 32 GB RAM. The NLP solvers used are IPM-based IPOPT and SQP-based SNOPT. The initial value of the convergence depth threshold ${\theta }_c^0$ is set to 0.6; the convergence depth advancement step $d_{\theta }$ is 0.1; the degree of progress threshold ${\theta }_p$ is set to 1.0 (the optimization algorithm has made no progress in the objective function and feasibility improvement); the maximum convergence depth ${\theta }_{max}$ is 1.0 (the iteration point is already a good approximation of the optimal point). The configuration points are of order 3 (k = 3), the finite elements are 200 (N = 200), and the discrete points are 600 in total.

5.1. Simulation Scenario of CHP Unit Variable Load Control

In this paper, simulation analysis is performed using the parameters of the extracted CHP unit in the literature [17], and the control system parameters of the CHP unit are shown in

Table 1. CHP unit control system parameters.

VT Control Parameters	VB Control Parameters	VH Control Parameters	Coordination Parameters
KPT = −0.9	KPB = 0.1	KPH = 0.01	K = 0.5
KIT = −1	KIB = 0.01	KIH = 10	TC = 15
KDT = 0	KDB = 0	KDH = 0

. The data of three cases of common variable load requirements for CHP units in integrated energy system applications are shown in

Table 2. Variable load requirements for CHP units under different operating conditions.

Case	$\mathbf{Initial} \mathbf{Value} ({\mathit{P}}_{\mathit{H}}^0,{\mathit{\psi }}_{\mathit{t}}^0,{\mathit{m}}_{\mathit{H}}^0)$	$\mathbf{Terminal} \mathbf{Value} ({\mathit{P}}_{\mathit{H}}^{\mathit{f}},{\mathit{\psi }}_{\mathit{t}}^{\mathit{f}},{\mathit{m}}_{\mathit{H}}^{\mathit{f}})$
Case 1 (Heating condition)	(220, 16.67, 380)	(260, 16.67, 390)
Case 2 (Pure condensing condition)	(260, 16.67, 0)	(235, 16.67, 0)

and

Table 3. Continuous variable load demand for CHP units under heating conditions (Case 3).

Time Period	$({\mathit{P}}_{\mathit{H}}^0,{\mathit{\psi }}_{\mathit{t}}^0,{\mathit{m}}_{\mathit{H}}^0)$	$({\mathit{P}}_{\mathit{H}}^{\mathit{f}},{\mathit{\psi }}_{\mathit{t}}^{\mathit{f}},{\mathit{m}}_{\mathit{H}}^{\mathit{f}})$
300–600 s	(250, 16.67, 390)	(220, 16.67, 380)
600–900 s	(220, 16.67, 380)	(200, 16.67, 390)
900–1200 s	(200, 16.67, 390)	(230, 16.67, 395)
1200–1500 s	(230, 16.67, 395)	(210, 16.67, 385)
1500–1800 s	(210, 16.67, 385)	(240, 16.67, 395)
1800–2100 s	(240, 16.67, 395)	(250, 16.67, 385)
2100–2400 s	(250, 16.67, 385)	(220, 16.67, 380)
2400–2700 s	(220, 16.67, 380)	(200, 16.67, 390)
2700–3000 s	(200, 16.67, 390)	(230, 16.67, 395)
3000–3300 s	(230, 16.67, 395)	(210, 16.67, 385)
3300–3600 s	(210, 16.67, 385)	(240, 16.67, 395)
3600–3900 s	(240, 16.67, 395)	(260, 16.67, 385)

. The adjustment time of the control process from the initial value to the terminal value is 300 s. Three different optimization solution algorithms are set for the simulation calculation. Method 1 is IPM; Method 2 is CDC-based IPM; Method 3 is CDC-based IPM-SQP.

5.2. Simulation Results Analysis

5.2.1. Comparative Analysis of the Effect of Variable Load Simulation

Simulation calculations were carried out by three methods for Case 1 and Case 2, and the results are shown in Figure 4. Different methods are used to solve the CHP unit variable load process under heating and pure condensing conditions. It can be seen that the output variables (${\psi }_t,P_{H,}{m}_H$) curves are consistent, and the dynamic process’s control variables (V_T, V_B, V_H) have the same trajectory. All of them achieve the optimal control effect.

The mechanistic model-related parameters of Equations (3)–(10) are taken at different values for CHP units under different operating conditions, resulting in the variation of DAEs of the unit model in the dynamic optimization problem. The simultaneous solution framework and algorithm improvement method proposed in this paper can converge and achieve the demanded effect when solving the unit variable load problem under different operating conditions. The simulation results prove that the unified description framework established and the optimization performance objectives designed in this paper for the dynamic optimization problems of CHP units with different operating conditions and different variable load demand scenarios have the generality to be applied to the integrated energy system calculation scenarios.

For the demand of continuously changing load in the application scenario of an integrated energy system, the control process of continuously changing electric and thermal loads for 12 time periods (300 s each) under the heating condition of the CHP unit is selected for simulation. The method proposed in this paper is used for the continuous dynamic optimization solution, and the output quantity curve and control quantity trajectory in this process are shown in Figure 5. The results show that the terminal values of electric power P_H and thermal power (related variable m_H) of the CHP unit can accurately realize the variable load demand. The unit’s main steam pressure ${\psi }_t$ and heating steam flow m_H curves continuously portray the pressure and heat fluctuation changes caused by the variable load process, which is conducive to the analysis of the integrated energy system calculation scenario for safe and stable operation.

Figure 5b shows that the control variables (V_T, V_B, V_H) trajectories enable coherent planning for optimal control of the variable load dynamic process in CHP units. These trajectories provide precise quantitative indicators for repeated adjustment of flow rate and valve opening during unit operation and guide actual CHP units’ continuous variable load regulation.

5.2.2. Method Performance Comparison Analysis

Different solution methods are used to calculate the optimal control problem of the CHP unit for each of the three case scenarios, and a set of optimal indicators are solved to obtain variable load control process paths. Their computational performance indexes are shown in

Table 4. Comparison of computational performance of different methods.

Case	Performance Indicators	CDC-IPM-SQP (IPOTP-SNOPT)	CDC-IPM (IPOTP)	IPM (IPOTP)
Case 1	Solution time	1.34 s	2.27 s	3.35 s
	Iteration number	33	47	61
	Feasibility error	2.03 × 10−5	1.13 × 10−5	5.13 × 10−8
	Convergence	Successful	Successful	Successful
	Terminal value	(259.9, 16.67, 390)	(259.9, 16.67, 390)	(260, 16.67, 390)
Case 2	Solution time	1.11 s	2.32 s	3.07 s
	Iteration number	23	41	61
	Feasibility error	2.71 × 10−5	1.25 × 10−5	5.93 × 10−8
	Convergence	Successful	Successful	Successful
	Terminal value	(234.9, 16.67, 0)	(234.9, 16.67, 0)	(235, 16.67, 0)
Case 3	Solution time	10.17 s	19.92 s	38.68 s
	Convergence	Successful in all periods	Successful in all periods	Successful in all periods
	Terminal value	All periods are satisfied	All periods are satisfied	All periods are satisfied

. It can be seen that the optimal control terminal values of the three methods do not differ much under different working conditions and variable load demands, and the feasibility error also meets the actual requirements. However, the number of iterations and the solution time of the discrete model have significant differences.

When the CDC-based improved method is used to solve the problem, the time interval between two calculations can be shorter, and the deviation between the initial value and the optimal solution becomes smaller, thus reducing the number of iterations in the calculation process. This feature can effectively improve the computational efficiency of the algorithm and greatly reduce the computing time. In case 1, compared with the IPM solution alone, the CDC criterion used in the CDC-IPM solution makes the constraint violation slightly larger than the result of the traditional H criterion. In addition, the objective function value is slightly smaller than the optimal value. However, the computation time of the CDC-IPM solution is reduced by about 30%, and the number of iterations is reduced by about 20%. Under the traditional H-criterion, the termination position of the optimization process (objective function and maximum constraint violation) is the 61st iteration, while the CDC criterion terminates at the 47th iteration. The simulation results show that the CDC criterion can effectively determine the degree of convergence of the solution process, the degree of the constraint violation, and the degree of improvement by continuing the iteration and can terminate the iteration in time to avoid over-convergence.

When the CDC-IPM-SQP method is used to solve the optimization problem of case 1, IPM reaches a convergence depth of 0.67 after 19 iterations, which exceeds the convergence depth threshold (0.6); the algorithm considers that the current iteration point is close to the optimal point and switches to the SQP algorithm; after continuing to iterate for six steps, SQP finds the optimal solution, and the algorithm terminates; the maximum value of the progress degree of SQP during the iteration process is 0.26, which is far from the threshold (1.0), so no switching is performed. Compared with the CDC-IPM solution, its computation time is shortened by about 40%, and the number of iterations is reduced by about 30%.

The improved IPM-SQP alternating solution method based on the CDC criterion can reduce the calculation time of variable load for the 12 time periods in case 3 by about 70%, which effectively improves the real-time performance of the scenario application.

6. Conclusions

In this paper, a method based on alternating SQP-IPM to solve the optimization of the variable-load dynamic process of CHP units is proposed. The computational performance of the CHP unit variable-load dynamic optimization under the three different solution methods is compared through numerical simulations of different arithmetic cases and the following conclusions are drawn:

(1)

Compared with the IPM solution alone, the CDC-IPM solution method uses the CDC criterion so that the constraint violation is slightly larger than the result of the conventional H criterion, while the objective function value is slightly smaller than the optimal value. However, the computation time of CDC-IPM is reduced by about 30%, and the number of iterations is reduced by about 20%.

(2)

Compared with the CDC-IPM solution, the computation time is reduced by about 40%, and the number of iterations is reduced by about 30% using the CDC-IPM-SQP.

(3)

The use of CDC criterion improvement and IPM-SQP alternate solution reduces the calculation time of variable load for 12 consecutive periods by about 70%, which effectively improves the real-time performance of energy system optimization scenario applications.

In summary, the dynamic optimization model and solution method for a variable load of CHP units proposed in this paper can improve the calculation speed while guaranteeing the solution accuracy and provide the possibility of incorporating the finely modeled CHP units into the energy scenario optimization application. This paper compares three solution methods for variable-load process optimization problems on specific CHP units, and further validation is needed for other types of units. It does not consider the technical difficulty of unit optimization and control modification, maintenance cost, and time scale of scenario application. Further discussion and research are needed in practical applications.