An Improved Finite Element Meshing Strategy for Dynamic Optimization Problems

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.


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][10][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.

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 [0,   ].The simplified model can be written as the following equations: min where  is the scalar objective (2)

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 2 − 1 order.When using  Radau points to configure, the algebraic precision is lower than the Gauss point, 2 − 2. 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 from Here  is constant depending on collocation points. nc is differential state equation at From (3),   will be zero at the collocations.So we need to select a noncollocation point  nc ; that is,  nc =   + ℎ   nc .Specifically, ( nc )/ can be stated as where Ω( nc ) is the differential of Ω  (  ).Then we add (3)-( 4) to problem (1) and then resolve the OPC problem.

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 [ 0 ,   ] 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 radio  = log(|error  / |) < 0.5, 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 radio  ≥ 0.5, add  finite element; here  rounds up from radio  .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 tol is set to allowable errors, it is set as max , (0.001, 0.001‖ , ‖).If the meshing method proposed by us can not satisfy the above Hamilton function error requirement, we must readd the finite element mesh quantity.

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  1 = 0.155051,  2 = 0.644949, and  3 = 1, noncollocation point is  nc = 0.5, finite element tolerance error is  = 1.0 − 5, and default set point error estimation limit is  = 1.0 − 5.

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 min −  (1)   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.

Rayleigh Problem.
Here we consider another optimal control problem which can be described as min ∫ 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.seawater (SWRO) desalination process.The model equation can be expressed as

Reverse Osmosis Sweater Desalination
And then according to the objective function of the cost of water production, min ,,, Here OC IP and OC EN are reverse osmosis process water intake and RO energy costs; OC CH is the cost of system chemicals; OC OTR 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,  elc is the unit price;  IP is the motor energy efficiency; LF represents the load factor;  in 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 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.

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    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.The number of iterations calculated ESM: Equally spaced method MFE: Moving finite elements   : Membrane water permeability (m⋅s −1 ⋅Pa −1 )  0 : Intrinsic membrane water permeability (m⋅s −1 ⋅Pa −1 )   : Membrane TDS permeability (m/s)  0 : Intrinsic membrane TDS permeability (m/s)  1 : Constant parameter  2 : Constant parameter  1 : Constant parameter : Feed (operational) temperature (K)   : Feed pressure (bar)   : Pressure drop along RO spiral wound module (bar) V: Solvent flux (kg/m 2 ⋅s) Δ: Pressure loss of osmosis pressure (bar) : Solute flux (kg/m 2 ⋅s) : Salt concentration of membrane surface (kg/m 3 ) : Permeate concentration of RO unit (kg/m 3 ) : Bulk concentration along feed channel (kg/m 3 )   : Mass transfer coefficient (m/s) Re: Reynolds number (dimensionless) Sc: Schmidt number (dimensionless)   : Hydraulic diameter of the feed spacer channel (m)   : Dynamicviscosity(m 2 /s) : Kinematic viscosity : Density of permeate water (kg/m 3 ) : Friction factor   : Empirical parameters : Axial velocity in feed channel (m/s) ℎ sp : Height of the feed spacer channel (m) OC: Operational cost of entire process OC EN : Energy cost of RO process OC IP : Energy cost for intake and pretreatment OC CH : Cost for chemicals OC OTR : Other costs SEC: Specific energy consumption (kw⋅h/m 3 )   : Permeateflowrate(m 3 /h)  elc : Electricityprice(CNY/kw⋅h)  in : Output pressure of the intake pump (bar) PLF: Load coefficient  IP : Mechanicalefficiencyofintakepump.

Figure 1 :
Figure 1: Equally spaced method for finite element mesh.

Figure 3 :
Figure 3: Error of differential profiles at finite elements with equally spaced method.

Figure 4 :
Figure 4: Error of differential profiles at finite elements with new method.

Figure 5 :Figure 6 :
Figure 5: Size of finite element in simulation and optimization.

Figure 7 :
Figure 7: State profiles for the batch reaction problem.

Figure 9 :
Figure 9: Error of differential profiles at finite elements with equally spaced method.

Figure 10 :
Figure 10: Error of differential profiles at finite elements with new method.

Figure 11 :Figure 12 :
Figure 11: Size of finite element in simulation and optimization.

Figure 13 :
Figure 13: State profiles for Rayleigh problem.

Figure 15 :
Figure 15: Error of differential profiles at finite elements with equally spaced method.The black dotted line represents the maximum error of each finite element at the noncollocation point.

Figure 16 :
Figure 16: Error of differential profiles at finite elements with new method.The black dotted line represents the maximum error of each finite element at the noncollocation point.

Figure 17 :
Figure 17: Size of finite element in simulation and optimization.

Figure 18 :
Figure 18: Membrane channel within one hour of the velocity changes.

3 )Figure 19 :
Figure 19: Membrane channel within one hour of the concentration changes.

Figure 20 :
Figure 20: Membrane channel within one hour of the pressure changes.
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   +1 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   =  −1 + (  −  −1 )   ,   (  ) =   ,   (  ) =   .

Table 1 :
Batch reaction numerical results.