Abstract
The finite element orthogonal collocation method is widely used in the discretization of differential algebraic equations (DAEs), while the discrete strategy significantly affects the accuracy and efficiency of the results. In this work, a finite element meshing method with error estimation on noncollocation point is proposed and several cases were studied. Firstly, the simultaneous strategy based on the finite element is used to transform the differential and algebraic optimization problems (DAOPs) into large scale nonlinear programming problems. Then, the state variables of the reaction process are obtained by simulating with fixed control variables. The noncollocation points are introduced to compute the error estimates of the state variables at noncollocation points. Finally, in order to improve the computational accuracy with less finite element, moving finite element strategy was used for dynamically adjusting the length of finite element appropriately to satisfy the set margin of error. The proposed strategy is applied to two classical control problems and a large scale reverse osmosis seawater desalination process. Computing result shows that the proposed strategy can effectively reduce the computing effort with satisfied accuracy for dynamic optimization problems.
1. Introduction
The direct transcription method is an important method to solve the problem of optimal control. By discretization of differential and algebraic optimization problems (DAOPs), the state variables and control variables are completely discretized. The discrete method uses the finite element orthogonal collocation, and generally the number of finite elements is empirically selected and the length of each finite element is equally divided. This results in low discretization accuracy for state and control variables, and to guarantee the satisfactory accuracy for some problems, the calculation time is too long to accept. Moving finite element strategy is a good idea for the solution. With the need of discrete differential algebraic equations, moving finite element is becoming the popular and practical technique for chemical process.
At present, we are concerned with calculation accuracy not only of the material, but also of the time in the process parameters for the chemical process attention. Modeling methods based on first principle and data-driven are used for practical control. With the development of solving technology, we can better understand the changes of state variables in the whole process of chemical reaction.
Betts and Kolmanovsky proposed a refinement procedure for nonlinear programming for discrete processes and estimating the discretization error for state variables [1]. Liu et al. used a novel penalty method to deal with nonlinear dynamic optimization problems with inequality path constraints [2]. Paiva and Fontes studied the adaptive mesh refinement algorithms which allow a nonuniform node collocation and apply a time mesh refinement strategy based on the local error into practical problems [3]. Zhao and Tsiotras introduced an efficient and simple method based on density (or monitor) functions, which have been used extensively for grid refinement. The accuracy and stability of the method were improved by choosing the appropriate density function. [4]. When using the discrete method to solve the nonlinear problem, iterative programming (IDP) algorithm is rather vulnerable to the stage of time in several aspects such as accuracy and the convergence rate. Li et al. introduce a self-adaptive variable-step IDP algorithm, taking account of the performance and control variables [5]. Zhang et al. presented an adaptive variable-step-size RKF method based on norm control. The method automatically adjusts the step size and the case where the local error norm value is 0 or the minimum value is discussed [6]. In view of the optimization of biochemical processes, some researchers have proposed adaptive parameterization methods to solve such problems [7, 8].
The above work is quite helpful to quick and stable solution of the dynamic optimization problem, but most of them put emphasis on method of choice and do not propose a specific operational process. As we know, the results of the optimization problem are often dependent on the specific algorithm.
This paper is a mesh refinement strategy based on the variable finite element mesh method [9–11]. Prior to optimization, a suitable initial mesh and length are obtained by GAMS platform simulation. Then, to optimize the objective function, the finite element is moved appropriately under the condition that the set error is satisfied. The method is applied to the chemical reaction of simple ODE equation [12, 13] and large scale reverse osmosis seawater desalination process, which proves its validity and feasibility.
2. Finite Element Discrete Method for Equation Model
According to the whole chemical reaction process, we want to know the changes in the material. Discrete solution of the model is the key part for computing process. Here, we introduce the orthogonal polynomials based on Lagrange’s use of the Radau collocation points on the finite element [14] and use of simultaneous method to the entire time domain . The simplified model can be written as the following equations:where is the scalar objective function, and as differential and algebraic variables, respectively; is the constraint equation for the state variables and the control variables, and as a control variable and time independent optimization variable. is a polynomial of order in collocation . denotes the length of element . is value of its first derivative in element at the collocation point . The advantage of the Lagrange interpolation polynomial over other interpolation methods is that the value of the variable at each collocation point is exactly equal to its coefficient
3. Error Calculation at Noncollocation Point
On the selection of collocation points on finite elements, if the mathematical expression is configured according to the Gauss point, the algebraic precision of the numerical integration is the highest. Gauss points can achieve the algebraic precision of order. When using Radau points to configure, the algebraic precision is lower than the Gauss point, . The discretization proposition based on Radau point has better stability [15]. Therefore, this paper uses Radau points to configure.
In order to improve the accuracy of the solution, Vasantharajan and Biegler [16] inserted the noncollocation points on the finite element to determine the increase and decrease of the finite element by the error of the noncollocation point.
Generally, in the solution of problem (1), it will cause compute error on state or control profiles due to finite mesh. In this study, the error of collocation finite elements can be estimated fromHere is constant depending on collocation points. is differential state equation at
From (3), will be zero at the collocations. So we need to select a noncollocation point ; that is, . Specifically, can be stated aswhere is the differential of . Then we add (3)-(4) to problem (1) and then resolve the OPC problem.
4. Finite Element Method
As the error estimation gotten with (3) and (4), we can divide the finite elements mesh. Generally, equally spaced method was widely used in many research, due to its simply operation. The flow chart of calculation is shown in Figure 1.
On the contrary, this approach will take more time and get larger after the mesh divided. Now we take an idea of adaptive mesh and obtain algorithm of mesh initialization.
Step 1. Divide the time domain equally with a coarse mesh and compute the error at every location of element by solving problem (1).
Step 2. When each finite element error , end the initialization of the finite element. If the requirement of less than 0.5 is not met, turn to Step 3.
Step 3. Once the error of each finite element , add finite element; here rounds up from . Then resolve problem (1).
However, the finite element mesh gotten above has not satisfied accuracy yet. So we first fix the mesh and suitable control profile and take simulate calculation by solving problem (1) without objective function. After that, relax and add constraint (7) as follows:
Through this way, seasonable initial values for optimization are provided. What is more, it also offers some freedom to search the breakpoint of control profiles.
From now on, we can summarize the approach described with program flow chart in Figure 2.
Through the above grid division and GAMS platform simulation and optimization, the optimized result is obtained finally. However, the optimality of the above result should be verified. According to the optimal control theory, the value of the Hamiltonian function comprising the objective function and the constraints should be kept as constant along the time axis. According to the optimal control theory, the discrete Hamiltonian function is denoted as
When the solution is optimal, the last two terms of the Hamiltonian are zero. According to the optimal control theory, when a system is optimal, the Hamiltonian must be a constant:
Therefore, if the finite element and the configuration point are not constant, the result is not optimal, and the meshing needs to readd the finite element. Here we use method of Tanartkit and Biegler [11], who put forward a Hamilton function criterion:
When is set to allowable errors, it is set as . If the meshing method proposed by us can not satisfy the above Hamilton function error requirement, we must readd the finite element mesh quantity.
5. Case Study
This numerical experiment is based on Intel (R) Core (TM) i3-2350M 2.3 GHz processor with 4 GB RAM; the nonlinear programming solver IPOPT [17, 18] is developed by Carnegie Mellon University.
To illustrate the above strategies, we consider two classical process dynamic optimization problems and a large scale reverse osmosis seawater (SWRO) desalination process. At the same time, we use the error analysis based on the equipartition finite element method for comparison. For the following example, we use the 3-order Radau collocation points; each finite element has a noncollocation point. Three of these collocation points are , , and , noncollocation point is , finite element tolerance error is , and default set point error estimation limit is .
5.1. Batch Temperature Profile with Parallel Reactions
This example is a chemical process, considering a nonisothermal reactor with first order , , where is reaction material, is target product, and is by-product. Now the goal is to find a feasible temperature profile that maximizes the final amount of product after one hour. The optimal control problem can be stated as
The problem has been solved in many pieces of literature. In [19], adaptive iterative dynamic programming leads to 50 time steps and needs more than 60 CPUs. Here we first adopt equally spaced elements method, which leads to 14 finite elements. Local errors for all finite elements are shown in Figure 3.
As shown in Figure 3, we can see that the first- and second-finite element errors are large, and other places have to meet the set error tolerance. This means that we can use fewer finite elements to calculate the problem. In Figure 4, when using our proposed meshing method, only 7 finite elements are needed. In the use of the equally finite element method, the posterior finite element error is very small, so that when we use a partitioning method, the error is close to the set error tolerance.
In Figure 5, fixed finite element length is gotten from the simulation. After adding the boundary constraint (7), a new finite element length is obtained (blue line). This can also help for locating breakpoints into specific problems.
Figures 6 and 7 show the optimal control image and the state variable change image, respectively. These changes are close to those obtained by previous researchers. This fully proves the feasibility of the proposed method.
Figure 8 shows that the Hamilton function is a constant and the derivative of Hamiltonian function satisfies the above error that the problem has been the optimal solution.
From Table 1, we can see that, by using an improved finite element method we propose, the finite element number can be reduced by half. The calculation time can be reduced by 4.5 s under the same precision.
5.2. Rayleigh Problem
Here we consider another optimal control problem which can be described as
In this example, the number of finite elements is often estimated by trial and error. Here, we use 104 equally spaced finite elements to make the two differential equations meet the error requirements, as shown in Figure 9.
As can be seen from Figure 9, most of the finite elements have been early to meet the error requirements. Here we use the new method we propose. Only 8 finite elements are required to satisfy the error requirements, as shown in Figure 10.
Figure 11 shows the length of each finite element position. It can be seen from the figure that the length of the 6, 7, and 8th finite element position is greater than the length of the simulation (red line), which indicates that the addition of the boundary constraint is effective when optimizing the objective function and finding the boundary breakpoint.
Figures 12 and 13 show the optimal control image and the state variable change image in process.
The Hamiltonian function image of Figure 14 indicates that the objective function has been optimized.
In Table 2, when using the finite element method, the total number of variables is 36, and the computer time is 67.277 s. But using our new method, the total number of discrete variables is 8, and the calculation time is 5.352 s. In terms of the number of finite elements and the computation time, the new method presented in this paper has some advantages.
5.3. Reverse Osmosis Sweater Desalination Process
The feasibility of the proposed method is illustrated by two classical chemical reaction control problems. The following is an operational optimization study of a large scale reverse osmosis seawater (SWRO) desalination process. The model equation can be expressed asAnd then according to the objective function of the cost of water production,
Here and are reverse osmosis process water intake and RO energy costs; is the cost of system chemicals; indicates other costs, including labor costs and maintenance costs.
The optimization proposition also needs to satisfy the RO process model, the reservoir model, the operation cost model, and the inequality constraint equation. The system cost model is shown as follows:
Here, is the unit price; is the motor energy efficiency; represents the load factor; is the outlet pressure of the pump.
From the optimization program mainly to the day, there are peaks and valleys of electricity prices, resulting in different periods of water production costs. The system can be done through the cistern reservoir volume, regulating the water system and the user’s demand for water to reduce the overall cost.
The reverse osmosis process model mainly includes three differential variables: membrane concentration, flow rate, and pressure. In this paper, the number of finite elements is divided into discrete errors based on the steady-state simulation process, so as to reduce the computational complexity of the optimization model.
Figure 15 shows three differential variables and the error of velocity, concentration, and pressure can only meet the error tolerance when the 22 finite elements are divided equally. Since the error of velocity and pressure is much smaller than that of channel concentration, the discrete error analysis was mainly focused on the channel concentration. Figure 16 is the finite element error analysis using the new method.
Figure 17 shows that after the initial finite element length is obtained by simulation, the constraint equation is added to obtain a more suitable finite element length.
Figures 18, 19, and 20 represent the flow rate, concentration, and pressure at the feed of the membrane module, respectively. It can be seen from these figures that adding boundary constraints to each curve is more smoothly and more responsive to changes in the variables within the actual membrane module.
Table 3 shows the use of a method proposed in this paper, both from the number of variables, the number of iterations, and computing time above a great upgrade.
It can be seen that the number of finite elements required for the model discretization is reduced to 14 when using the moving finite element method. Due to the finite element method, the finite element length can be adjusted to a small extent, which greatly increases the accuracy of the solution.
6. Conclusion
In this paper, we proposed a mesh-partitioning strategy based on the direct transcription method to solve the optimal control problem. This method discretizes the differential algebraic equation (DAE) using the Radau collocation point based on the variable finite element and finally transforms into a nonlinear programming problem. Here, the initialization variable finite element uses a valid termination criterion. At the same time, the computation time can be saved under the condition of satisfying certain precision.
Cases study of two classical control problems and a large scale reverse osmosis process optimization problem was carried with our proposed method and conventional method. The computing results show that the proposed method can effectively reduce the number of finite element and thus reduce the computing efforts with permitted discrete accuracy.
Nomenclature
| : | Terminal time (s) |
| : | The scalar objective function |
| : | Differential state variables |
| : | Algebraic state variables |
| : | Constraint equation for the state and control variables |
| : | Control variable |
| : | Time independent optimization variable |
| : | The length of element |
| : | Number of collocation points |
| : | First derivative in element at the collocation point |
| : | Constant depending on collocation points |
| : | Differential state equation at |
| : | Default set point error estimation limit |
| : | The setting error |
| : | Number of finite elements |
| Var: | Number of variables |
| Eq: | Number of equations |
| Iter: | The number of iterations calculated |
| ESM: | Equally spaced method |
| MFE: | Moving finite elements |
| : | Membrane water permeability (m·s−1·Pa−1) |
| : | Intrinsic membrane water permeability (m·s−1·Pa−1) |
| : | Membrane TDS permeability (m/s) |
| : | Intrinsic membrane TDS permeability (m/s) |
| : | Constant parameter |
| : | Constant parameter |
| : | Constant parameter |
| : | Feed (operational) temperature (K) |
| : | Feed pressure (bar) |
| : | Pressure drop along RO spiral wound module (bar) |
| : | Solvent flux (kg/m2·s) |
| : | Pressure loss of osmosis pressure (bar) |
| : | Solute flux (kg/m2·s) |
| : | Salt concentration of membrane surface (kg/m3) |
| : | Permeate concentration of RO unit (kg/m3) |
| : | Bulk concentration along feed channel (kg/m3) |
| : | Mass transfer coefficient (m/s) |
| : | Reynolds number (dimensionless) |
| : | Schmidt number (dimensionless) |
| : | Hydraulic diameter of the feed spacer channel (m) |
| : | Dynamic viscosity (m2/s) |
| : | Kinematic viscosity |
| : | Density of permeate water (kg/m3) |
| : | Friction factor |
| : | Empirical parameters |
| : | Axial velocity in feed channel (m/s) |
| : | Height of the feed spacer channel (m) |
| OC: | Operational cost of entire process |
| : | Energy cost of RO process |
| : | Energy cost for intake and pretreatment |
| : | Cost for chemicals |
| : | Other costs |
| SEC: | Specific energy consumption (kw·h/m3) |
| : | Permeate flow rate (m3/h) |
| : | Electricity price (CNY/kw·h) |
| : | Output pressure of the intake pump (bar) |
| PLF: | Load coefficient |
| : | Mechanical efficiency of intake pump. |
Conflicts of Interest
The authors declare no conflicts of interest.
Acknowledgments
This work is supported by National Natural Science Foundation (NNSF) of China (no. 61374142), the Natural Science Foundation of Zhejiang (LY16F030006), the Public Projects of Zhejiang Province, China (no. 2017C31065), Research and Innovation Fund of Hangzhou Dianzi University (CXJJ2016043), and Hangzhou Dianzi University Graduate Core Curriculum Construction Project (GK168800299024-012).