Two-Layer Predictive Control of a Continuous Biodiesel Transesterification Reactor

Abstract

A novel two-layer predictive control scheme for a continuous biodiesel transesterification reactor is presented. Based on a validated mechanistic model, the least squares (LS) algorithm is used to identify the finite step response (FSR) process model adapted in the controller. The two-layer predictive control method achieves the steady-state optimal setpoints and resolves the multivariable dynamic control problems synchronously. Simulation results show that the two-layer predictive control strategy leads to a significant improvement of control performance in terms of the optimal set-points tracking and disturbances rejection, as compared to conventional PID controller within a multiloop framework.

1. Introduction

With the depletion of fossil fuels and global environmental degradation, the development of alternative fuels from renewable resources has received considerable attention. Biodiesel has become the foremost alternative fuel to those refined from petroleum products. It can be produced from renewable sources, such as vegetable and animal oils, as well as from wastes, such as used cooking oil. Transesterification is the primary method of converting these oils to biodiesel [1–3]. A block diagram for a biodiesel production process by transesterification is shown in Figure 1.

Figure 1 

Biodiesel production by transesterification.

A modern transesterification plant is continuous instead of batch. A continuous biodiesel production leads to better heat economization, better product purity from phase separation by removing only the portion of the layer furthest from the interface, better recovery of excess methanol in order to save on methanol cost and regulatory issues, minimal operator interference in adjusting plant parameters, and lower capital costs per unit of biodiesel produced. The same technology can also be applied to other biofuels production [4–6].

Biodiesel transesterification reactor is the most crucial operation unit to be controlled because any drift in standard operating condition may lead to significant changes in process variable and production quality specification [4–7]. These reactors have complicated dynamics and heat transfer characteristics. Moreover, they are inherently concerned with nonlinearity which arises from fluctuations of reactant concentration, reactant temperature, coolant temperature, and instrumentation noise or complex microbial interactions. The complicated nonlinear, multivariable, and coupling in nature are the fundamental control problems involved in biodiesel reactor [8, 9].

Recently, a number of reports have appeared on the controller design and dynamic optimization in continuous and batch biodiesel reactors. Mjalli et al. developed a rigorous mechanistic model of a continuous biodiesel reactor and proposed a multimodel adaptive control strategy which realized the set-point tracking and disturbance rejection [4]. Ho et al. further adopted adaptive generalized predictive control strategy to handle multivariable problems of a biodiesel reactor [8]. Wali et al. proposed an artificial intelligence technique to design online genetic-ANFIS temperature control based on LabVIEW for a novel continuous microwave biodiesel reactor [10]. Benavides and Diwekar realized the optimal control of a batch biodiesel reactor involved optimization of the concentration based on maximum principle [11].

This work considers the advanced control strategy of biodiesel continuous transesterification reactor. Model predictive control (MPC) is one of the most popular advanced control strategies. It is a class of model-based control algorithm, which has become a complex standard process industry solving complicated constrained multivariable control problems, and widely used in the chemical and petrochemical processes [12]. The main technical characteristics of MPC, include using mathematical models and history input and output data to predict future output, combined with the established control objectives, to calculate the optimal feedback rate. Compared with the traditional multiloop PID controllers, MPC takes into account simultaneously the effects of all manipulated variables to all controlled variables. Usually successfully put into operation, MPC can significantly reduce the standard deviation of the controlled variable and then through the card edge operations, improve the overall efficiency of the control system.

In recent years, there has been an integrated steady-state optimization of the two-layer predictive control strategy in MPC industry technology [13–15]. Two-layer predictive control is divided into upper steady-state optimization (SSO) layer and lower dynamic control layer. SSO can achieve real time optimization (RTO) objectives tracking asymptotically, independently complete local economic optimization of the corresponding MPC procedure. Specifically, the upper SSO uses steady-state gain of MPC dynamic mathematical model as the mathematical model and searches the optimum value within the constraints space of MPC. Part steady-state values of the operating or output variables will be in the position of “card edge”. The calculation results of the SSO layer will be as the set- points to the lower MPC layer.

Although two-layer predictive control strategy has been widely used in many applications of chemical reactors, hardly any work was done on the biodiesel transesterification reactor. In this paper, a two-layer predictive control strategy is designed, tested, and simulated on a continuous biodiesel transesterification reactor. The scheme can amplify the advantages of both technologies in terms of process stability, and optimal and improved performances. Section 2 discusses the transesterification mechanism, which uses a validated mechanistic model of Mjalli et al. [4]. Then the two-layer predictive control strategy is developed in Section 3. Section 4 gives the control system design based on two-layer predictive control theory. Section 5 discusses model identification results and the performances of the control strategy.

2. Mathematical Models

The modeling of transesterification reactors starts with understanding the complex reaction kinetic mechanism. The stoichiometry of vegetable oil methanolysis reaction requires three mol of methanol (A) and one mol of triglyceride (TG) to give three mol of fatty acid methyl ester (E) and one mol of glycerol (G) [16]. The overall reaction scheme for this reaction is TG+3A⟷3E+G. (1)

The methanolysis, in turn, consists of three consecutive reversible reactions, where a mole of fatty acid methyl ester is released in each step, and monoglycerides (MG) and diglycerides (DG) are intermediate products. The stepwise reactions are(2)

The stepwise reactions can be termed as pseudo-homogeneous catalyzed reactions, following second-order kinetics. The second-order kinetic model can be explained through the following set of differential equations [17]: 𝑑𝐶_TG𝑑𝑡=−𝑘^′₁𝐶_TG𝐶_A+𝑘^′₂𝐶_DG𝐶_E,𝑑𝐶_DG𝑑𝑡=𝑘^′₁𝐶_TG𝐶_A−𝑘^′₂𝐶_DG𝐶_E−𝑘^′₃𝐶_DG𝐶_A+𝑘^′₄𝐶_MG𝐶_E,𝑑𝐶_MG𝑑𝑡=𝑘^′₃𝐶_DG𝐶_A−𝑘^′₄𝐶_MG𝐶_E−𝑘^′₅𝐶_MG𝐶_A+𝑘^′₆𝐶_GL𝐶_E,𝑑𝐶_E𝑑𝑡=𝑘^′₁𝐶_TG𝐶_A−𝑘^′₂𝐶_DG𝐶_E+𝑘^′₃𝐶_DG𝐶_A−𝑘^′₄𝐶_MG𝐶_E+𝑘^′₅𝐶_MG𝐶_A−𝑘^′₆𝐶_GL𝐶_E,𝑑𝐶_A𝑑𝑡=−𝑑𝐶_E,𝑑𝑡𝑑𝐶_GL𝑑𝑡=𝑘^′₅𝐶_MG𝐶_A−𝑘^′₆𝐶_GL𝐶_E, (3) where 𝐶_TG, 𝐶_DG, 𝐶_MG, 𝐶_E, 𝐶_A, and 𝐶_GL are concentrations of triglyceride, diglyceride, monoglyceride, methyl ester, methanol, and glycerol, respectively. 𝑘^′₁, 𝑘^′₃, and 𝑘^′₅ are the effective rate constants for the forward reactions, and 𝑘^′₂, 𝑘^′₄, and 𝑘^′₆ are the effective rate constants for the reverse reactions.

The previously selected kinetic model can be formulated in terms of a general reaction equation 𝑟_𝑗=𝑘^′_𝑗𝐶_𝑖². (4)

The catalyst concentration remained constant because the sidereactions that consume the catalyst were supposed to be negligible. Therefore, each effective rate constant includes the catalyst concentration (𝐶_cat) and the corresponding rate constant for the catalyzed reaction [18]: 𝑘^′_𝑗=𝑘_𝑗𝐶_cat. (5)

The temperature influence on the reaction rate was studied from the Arrhenius equation (6) that shows the temperature dependency of the reaction rate constant 𝑘_𝑗=𝑘₀𝑒^{(−𝐸𝑎/𝑅𝑇)}, (6) where 𝑘₀ is a constant called the preexponential factor, 𝐸_𝑎 is the activation energy of the reaction, and 𝑅 is the gas constant.

In order to realize the optimization and control of continuous biodiesel production process, the model used in the paper on the basis of the second-order kinetic model jointing the material and energy balance equations as well as the dynamic equation of the coolant temperature. The material balance for each component is expressed as follows [4]: 𝑉𝑑𝐶_𝑖𝑑𝑡=𝐹_𝑖0𝐶_𝑖0−𝐹_𝑖𝐶_𝑖−^𝑛_𝑗=1𝑟_𝑗𝑉. (7)

The reactor energy balance is expressed as 𝑉^𝑛𝑖_𝑖=1𝐶_𝑖𝐶_𝑝𝑖𝑑𝑇=𝐹𝑑𝑡^TG0𝐶_TG0𝐶_𝑃TG+𝐹_A0𝐶_A0𝐶_𝑃A𝑇₀−𝑉−𝑇_{𝑛𝑖=𝑗}𝑟_𝑗Δ𝐻_𝑗−𝑈𝐴_𝐻Δ𝑇. (8)

The coolant fluid energy balance is expressed as 𝑑𝑇_𝐶=𝐹𝑑𝑡𝑉^𝐶0_𝐶𝑇_𝐶0−𝑇_𝐶+𝑈𝐴_𝐻Δ𝑇𝜌_𝐶𝑉_𝐶𝐶_𝑃𝐶. (9)

The function equation of heat transfer coefficient is approximately expressed as 𝑈=𝛼𝐹_𝐶^𝛽𝑁^𝛾=735.5𝐹_𝐶^1.095𝑁^0.405. (10)

3. Theory of Two-Layer Predictive Control

In modern process industries, the MPC controller is part of a multilevel hierarchy of optimization and control functions. Typically it is three-layer structure; that is, an RTO block is at the top layer, a MPC block is at the middle, and a PID block is at the bottom [19]. Therefore, under this multilevel hierarchy control system structure, the primary task of the MPC is to dynamic track the computational target calculated by the RTO. RTO layer should be optimized for the whole device.

Reference [20] proposed the framework of two-layer predictive control shown in Figure 2. SSO is added between RTO and MPC. Left branch, the SSO layer is used for recalculating the results of RTO layer, make the output steady-state target be located in the steady state gain matrix column space, so as to meet the compatibility and consistency conditions of steady state solution. Right branch, the role of SSO is to conduct local optimization to further improve the MPC steady-state performance, which can effectively resolve the nonparty system setpoints in the given problem.

Figure 2 

Framework of two-layer predictive control of industrial processes.

Mathematical description of the two-layer predictive control include establishing steady-state mathematical model, steady-state target calculation, and a dynamic controller design [21].

3.1. Establish Steady-State Mathematical Model

Assume an MIMO plant with 𝑚 control input and 𝑝 controlled output and the coefficients of the corresponding step response model between control input 𝑢_𝑗 and output 𝑦_𝑖 are given; the model vector is 𝑎_𝑖𝑗(𝑡)=𝑎_𝑖𝑗(1),…,𝑎_𝑖𝑗(𝑁)^𝑇, (11) where 𝑖=1,…,𝑝;  𝑗=1,…,𝑚. 𝑁 in (11) denotes modeling horizon of step response model. Thus, a multistep predictive model can be obtained: 𝑦(𝑘+1)=𝑦(𝑘)+𝐴₁Δ𝑢(𝑘), (12) where 𝑦()=𝑦𝑘+1¹()⋮𝑦𝑘+1_𝑝()(𝑘)=𝑦𝑘+1;𝑦¹(𝑘)⋮𝑦_𝑝(𝑘);Δ𝑢(𝑘)=Δ𝑢₁(𝑘)⋮Δ𝑢_𝑚(𝑘);𝐴₁=𝑎₁₁(1)⋯𝑎_1𝑚(1)𝑎⋮⋮⋮_𝑝1(1)⋯𝑎_𝑝𝑚(1). (13)

Under the control increment Δ𝑢(𝑘),…,Δ𝑢(𝑘+𝑀−1) action, the output predictive value of the system is 𝑦(𝑘+1)=𝑦(𝑘)+𝐴₁(𝑘),𝑦()(𝑘)Δ𝑢𝑘+2=𝑦+𝐴₂Δ𝑢(𝑘)+𝐴₁(),⋮𝑦()(𝑘)Δ𝑢𝑘+1𝑘+𝑁=𝑦+𝐴_𝑁Δ𝑢(𝑘)+⋯+𝐴_{𝑁−𝑀+1}Δ𝑢(𝑘+𝑀−1), (14)

abbreviated as 𝜕𝑦(𝑘)=𝐴Δ𝑢_𝑀(𝑘), (15) where (𝑘)=𝑦()(𝑘)⋮𝑦()(𝑘)𝜕𝑦,𝑘+1−𝑦𝑘+𝑁−𝑦Δ𝑢_𝑀(𝑘)=(𝑘)⋮(),𝐴Δ𝑢Δ𝑢𝑘+𝑀−1𝐴=¹0𝐴⋮⋱_𝑀⋯𝐴₁𝐴⋮⋮_𝑁⋯𝐴_{𝑁−𝑀+1}. (16)

The system can be written at the steady-state time Δ𝑦(∞)=𝐴_𝑁Δ𝑢(∞), (17) where Δ𝑦(∞)=[Δ𝑦₁(∞),Δ𝑦₂(∞),…,Δ𝑦_𝑝(∞)]^𝑇, Δ𝑢(∞)=[Δ𝑢₁(∞),Δ𝑢₂(∞),…,Δ𝑢_𝑚(∞)]^𝑇 are the steady-state output increment and input increment, respectively, and 𝐴_𝑁 is the steady-state step response coefficients matrix 𝐴_𝑁=𝑎₁₁(𝑁)⋯𝑎_1𝑚(𝑁)𝑎⋮⋮⋮_𝑝1(𝑁)⋯𝑎_𝑝𝑚(𝑁). (18)

To meet the requirements of steady-state target calculation, model (17) can also be written as Δ𝑦_∞(𝑘)=𝐴_𝑁Δ𝑢_∞(𝑘). (19)

3.2. Steady-State Target Calculation

3.2.1. Basic Problem Description

Steady-state target calculation is to maximize economic benefits for the purpose of self-optimization under MPC existing configuration mode according to the process conditions. According to the production process characteristics and objectives, the basic problem of steady-state target calculation is the optimization process, which controlled input as cost variables, controlled output as steady-state variables. A common description of the objective function is as follows [21]: min_{Δ𝑢∞(𝑘),Δ𝑦∞(𝑘)}𝐽=𝛼^𝑇Δ𝑢_∞(𝑘)+𝛽^𝑇Δ𝑦_∞(𝑘). (20)

Since Δ𝑢_∞ and Δ𝑦_∞ are linearly related, the input output variation of objective function can be unified to control the input change. The formula (20) can be unified as min_{Δ𝑢∞(𝑘)}𝐽=𝑐^𝑇Δ𝑢_∞(𝑘), (21) where 𝑐^𝑇=[𝑐₁,…,𝑐_𝑚] is the cost coefficient vector, constructed by the normalized benefit, or cost of each input variable. Δ𝑢_∞(𝑘)=[Δ𝑢¹_∞,…,Δ𝑢^𝑚_∞]^𝑇 is the steady-state change value of every input at time 𝑘.

Given the steady-state constraints of input and output variables, global-optimization problem of steady-state target calculation can be described as the following linear program (LP) problem: min_{Δ𝑢∞(𝑘)}𝐽=𝑐^𝑇Δ𝑢_∞(𝑘)s.t.Δ𝑦_∞(𝑘)=𝐺_𝑢Δ𝑢_∞(𝑘)+𝐺_𝑓Δ𝑓_∞(𝑘)𝑢+𝑒,_min≤𝑢_∞(𝑘)+Δ𝑢_∞(𝑘)≤𝑢_max,𝑦_min≤𝑦_∞(𝑘)+Δ𝑦_∞(𝑘)≤𝑦_max, (22) where 𝐺_𝑢, 𝐺_𝑓 are the steady-state gain matrices of control input and disturbance variables; and 𝑒 is the model bias. 𝑢_min, 𝑢_max are low limit and upper limit of steady-state input variables 𝑦_min, 𝑦_max are low limit and upper limit of steady state output variables.

The global-optimization problem of steady-state target calculation can be described as the following quadratic program (QP) problem: min_{Δ𝑢∞(𝑘)}𝐽=𝑐^𝑇Δ𝑢_∞(𝑘)−Maxprofit²s.t.Δ𝑦_∞(𝑘)=𝐺_𝑢Δ𝑢_∞(𝑘)+𝐺_𝑓Δ𝑓_∞(𝑘)𝑢+𝑒,_min≤𝑢_∞(𝑘)+Δ𝑢_∞(𝑘)≤𝑢_max,𝑦_min≤𝑦_∞(𝑘)+Δ𝑦_∞(𝑘)≤𝑦_max, (23) where Maxprofit is the potential maximum economic profit.

3.2.2. Feasibility Judgment and Soft Constraint Adjustment

Mathematically, optimization feasibility is the existence problem of the optimal solution. Feasibility of steady-state target calculation means that optimal steady state of input-output should meet their operating constraints; if feasible solution does not exist, the optimization calculation has no solution. The solving process is as follows: first, judge the existence of space domain formed by the constraints and if there is in it for optimization, if does not exist, then through the soft constraints adjustment to obtain the feasible space domain, and then to solve.

Soft constraints adjustment is an effective way to solve infeasible optimization [22, 23]. By relaxing the output constraints within the hard constraints, increasing the optimization problem feasible region that feasible solution to be optimized. Hard constraints refer to unalterable constraints limited by the actual industrial process.

Engineering standards of the priority strategy of soft constraints adjustment are the following: give priority to meet the highly important operating constraints, and allow less important operating constraints to be violated appropriately under the premise of satisfying the engineering constraints.

Considering the following constraints (24), constituted by steady-state model input constraints and output constraints containing slack variables, the priority rank is “𝑁”, where Δ𝑦_∞(𝑘)=𝐺_𝑢Δ𝑢_∞(𝑘)+𝐺_𝑓Δ𝑓_∞(𝑘)𝑢+𝑒,_𝐿𝐿≤𝑢_∞(𝑘)+Δ𝑢_∞(𝑘)≤𝑢_𝐻𝐿,𝑦^𝑗_𝐿𝐿−𝜀^𝑗₂≤𝑦_∞(𝑘)+Δ𝑦_∞(𝑘)≤𝑦^𝑗_𝐻𝐿+𝜀^𝑗₁,𝜀^𝑗₁≥0,𝜀^𝑗₂𝜀≥0,_𝑗1≤𝑦_𝐻𝐻𝐿−𝑦_𝐻𝐿,𝜀^𝑗₂≤𝑦_𝐿𝐿−𝑦_𝐿𝐿𝐿,𝑗=1,…,𝑁. (24)

The algorithm steps of feasibility judgment and soft constraint adjustment based on the priority strategy are as follows.

Step 1. Initialization: according to the characteristics of the output variables and process conditions, set the upper and lower output constraints priority ranks, the same priority rank setting adjustments according to actual situation constraint weights.

Step 2. According to the priority ranks, judge the feasibility and adjust the soft constraints in accordance with the ranks from large to small. Under a larger priority rank if cannot find a feasible solution, the constraints of the rank will be relaxed to hard constraints, and then consider less priority rank constraints, until we find a feasible solution.

Step 3. Then the steady-state target calculation entered the stage of economy optimization or target tracking.

For Step 2, constraints of the highest priority rank 𝑁 are adjusted first by solving the following optimization problem: min_𝜀𝑁𝐽=𝑊^𝑁𝑇𝜀^𝑁,𝑊^𝑁𝑇=𝑊^𝑁₁,…,𝑊^𝑁_2×𝑛𝑁s.t.Θ^𝑁𝑍^𝑁=𝑏^𝑁,Ω^𝑁𝑍^𝑁≤Ψ^𝑁, (25) where𝑍^𝑁=𝑋^𝑇₁,𝑋^𝑇₂,𝑋¹₃^𝑇,…,𝑋^𝑁₃^𝑇,𝑋¹₄^𝑇𝑋,…,_𝑁4𝑇,𝜀^𝑁₁^𝑇,𝜀^𝑁₂^𝑇,𝜀^𝑁₁^𝑇,𝜀^𝑁₂^𝑇𝑇,Ω^𝑁=block-diag−𝐼_𝑚,−𝐼_𝑚,−𝐼_𝑛1,…,−𝐼_𝑛𝑁,−𝐼_𝑛1,…,−𝐼_𝑛𝑁,−𝐼_𝑛𝑁,−𝐼_𝑛𝑁,𝐼_𝑛𝑁,𝐼_𝑛𝑁,Ψ^𝑁=0_𝑚×1^𝑇,0_𝑚×1^𝑇,0_𝑛1×1^𝑇,…,0_𝑛𝑁×1^𝑇,0_𝑛1×1^𝑇,…,0_𝑛𝑁×1^𝑇,0_𝑛𝑁×1^𝑇,0_𝑛𝑁×1^𝑇,𝑦^𝑁_𝐻𝐻𝐿−𝑦^𝑁_𝐻𝐿^𝑇,𝑦^𝑁_𝐿𝐿−𝑦^𝑁_𝐿𝐿𝐿^𝑇𝑇,𝑏^𝑁=𝑢_𝐻𝐿−𝑢_𝐿𝐿𝐺¹_𝑢𝑢_∞(𝑘)−𝐺¹_𝑢𝑢_𝐿𝐿(𝑘)+𝑦¹_𝐻𝐿−𝑦¹_∞(𝑘)−𝐺¹_𝑓Δ𝑓_∞(𝑘)−𝑒¹⋮𝐺^𝑁_𝑢𝑢_∞(𝑘)−𝐺^𝑁_𝑢𝑢_𝐿𝐿(𝑘)+𝑦^𝑁_𝐻𝐿−𝑦^𝑁_∞(𝑘)−𝐺^𝑁_𝑓Δ𝑓_∞(𝑘)−𝑒^𝑁𝐺¹_𝑢𝑢_𝐻𝐿−𝐺¹_𝑢𝑈_∞(𝑘)+𝑌¹_∞(𝑘)+𝐺¹_𝑓Δ𝑓_∞(𝑘)−𝑦¹_𝐿𝐿+𝑒¹⋮𝐺^𝑁_𝑢𝑢_𝐻𝐿−𝐺^𝑁_𝑢𝑈_∞(𝑘)+𝑦^𝑁_∞(𝑘)+𝐺^𝑁_𝑓Δ𝑓_∞(𝑘)−𝑦^𝑁_𝐿𝐿+𝑒^𝑁00,Θ^𝑁=𝐼_𝑚𝐼_𝑚𝐺0⋯00⋯00000_1𝑢00⋯0𝐼_𝑛1𝐺000000⋮⋮0⋯00⋱00000_𝑁𝑢00⋯000𝐼_𝑛𝑁−𝐼_𝑛𝑁0000𝐺¹_𝑢𝐼_𝑛1000⋯00000⋮⋮0⋱00⋯000000𝐺^𝑁_𝑢00𝐼_𝑛𝑁0⋯00−𝐼_𝑛𝑁00000⋯00⋯0𝐼_𝑛𝑁0−𝐼_𝑛𝑁0000⋯00⋯00𝐼_𝑛𝑁0−𝐼_𝑛𝑁. (26)

Solving (25) may appear in three different cases, respectively: if (25) is feasible, and the optimum solution is 𝜀^𝑁=0, subject to 𝐽=0, that is, no need for soft constraints adjustment, directly solve the original problem (22); if (25) is feasible, but 𝜀^𝑁≠0, just need to relax constraints of priority ranks 𝑁, and further optimization solution; if (25) is infeasible, not get a feasible solution to soft constraints adjustment of the priority rank 𝑁, relaxing the constraints of the priority rank 𝑁 to hard constraints; that is, 𝜀^𝑁₁=𝑦^𝑁_𝐻𝐻𝐿−𝑦^𝑁_𝐻𝐿,𝜀^𝑁₂=𝑦^𝑁_𝐿𝐿−𝑦^𝑁_𝐿𝐿𝐿. (27)

Go to the procedure of judging rank 𝑁−1 constraints min_𝜀𝑁−1𝐽=𝑊^𝑁−1𝑇𝜀^𝑁−1,𝑊^𝑁−1𝑇=𝑊₁^𝑁−1,…,𝑊^𝑁−1_{2×𝑛𝑁−1}s.t.Θ^𝑁−1𝑍^𝑁−1=𝑏^𝑁−1,Ω^𝑁−1𝑍^𝑁−1≤Ψ^𝑁−1. (28)

For (28), the matrix form is the same with priority rank 𝑁, only in the corresponding position of 𝜀^𝑁−1 to replace 𝜀^𝑁, 𝑏^𝑁−1 matrix is adjusted𝑏^𝑁−1=𝑢_𝐻𝐿−𝑢_𝐿𝐿𝐺¹_𝑢𝑢_∞(𝑘)−𝐺¹_𝑢𝑢_𝐿𝐿(𝑘)+𝑦¹_𝐻𝐿−𝑦¹_∞(𝑘)−𝐺¹_𝑓Δ𝑓_∞(𝑘)−𝑒¹⋮𝐺^𝑁_𝑢𝑢_∞(𝑘)−𝐺^𝑁_𝑢𝑢_𝐿𝐿(𝑘)+𝑦^𝑁_𝐻𝐿−𝑦^𝑁_∞(𝑘)−𝐺^𝑁_𝑓Δ𝑓_∞(𝑘)−𝑒^𝑁+𝑦^𝑁_𝐻𝐻𝐿−𝑦^𝑁_𝐻𝐿^𝑇𝐺¹_𝑢𝑢_𝐻𝐿−𝐺¹_𝑢𝑢_∞(𝑘)+𝑦¹_∞(𝑘)+𝐺¹_𝑓Δ𝑓_∞(𝑘)−𝑦¹_𝐿𝐿+𝑒¹⋮𝐺^𝑁_𝑢𝑢_𝐻𝐿−𝐺^𝑁_𝑢𝑢_∞(𝑘)+𝑦^𝑁_∞(𝑘)+𝐺^𝑁_𝑓Δ𝑓_∞(𝑘)−𝑦^𝑁_𝐿𝐿+𝑒^𝑁+𝑦^𝑁_𝐿𝐿−𝑦^𝑁_𝐿𝐿𝐿^𝑇00. (29)

𝑁−1 rank and 𝑁 rank are the same for the soft constraints adjustment processing, until the end of constraint adjustment of the priority rank 1. If all ranks of constraints are relaxed to the hard constrain and a feasible solution still can't be found, then the original problem of soft constraints adjustment is infeasible and needs to be redesigned.

3.3. Dynamic Controller Design

In engineering applications, dynamic matrix control (DMC) algorithm is one of the most widely used MPC algorithms based on the step response model of the plant. This paper adopts DMC and steady-state target calculation integration strategy.

The difference is that the general DMC algorithms have no requirements on the steady-state position of the control input, and they only require the controlled output as close as possible to arrive at its set point. However, the integration strategy DMC requires both input and output variables to approach their steady-state targets (𝐮_𝑠,𝐲_𝑠) as far as possible. The algorithm has three basic characteristics: predictive model, receding horizon optimization, and feedback correction [24].

3.3.1. Predictive Model

Based on system process step response model, at the current time 𝑘, the future 𝑃-step prediction output can be written as follows: 𝐲_𝑃𝑀(𝑘)=𝐲_𝑃0(𝑘)+𝐀Δ𝐮_𝑀(𝑘), (30) where 𝑃 denotes the prediction horizon, 𝑀 is the control horizon, 𝐀 is the prediction matrix composed by the corresponding step response coefficients, 𝐲_𝑃0 is the initial output prediction value when control action starting from the present time does not change, Δ𝐮_𝑀(𝑘) is the prediction incremental in 𝑀 control horizon, and 𝐲_𝑃𝑀(𝑘) is the future 𝑃-step prediction output under 𝑀-step control action change. Among them 𝐲_𝑃𝑀(𝑘)=𝑦_1,𝑃𝑀(𝑘)⋮𝑦_{𝑝,𝑃𝑀}(𝑘),𝐲_𝑃0(𝑘)=𝑦_1,𝑃0(𝑘)⋮𝑦_𝑝,𝑃0(𝑘),Δ𝐮_𝑀(𝑘)=Δ𝑢_1,𝑀(𝑘)⋮Δ𝑢_𝑚,𝑀(𝑘)𝐴,𝐀=¹¹⋯𝐴_1𝑚𝐴⋮⋱⋮_𝑝1⋯𝐴_𝑝𝑚. (31)

3.3.2. Receding Horizon Optimization

In the receding horizon optimization process, control increment can be obtained in every execution cycle by minimizing the following performance index: min_{Δ𝐮𝑀(𝑘)}𝐽(𝑘)=𝐰(𝑘)−𝐲_𝑃𝑀(𝑘)²_𝐐+‖𝜺(𝑘)‖²_𝐒+𝐮_𝑀(𝑘)−𝐮_∞²_𝐓+Δ𝐮_𝑀(𝑘)²_𝐑. (32)

Subject to the model 𝐲_𝑃𝑀(𝑘)=𝐲_𝑃0(𝑘)+𝐀Δ𝐮_𝑀(𝑘). (33)

Subject to bound constraints 𝐲_min−𝜺≤𝐲_𝑃𝑀(𝑘)≤𝐲_max𝐮+𝜺,_min≤𝐮_𝑀≤𝐮_max,Δ𝐮_min≤Δ𝐮_𝑀(𝑘)≤Δ𝐮_max, (34) where 𝜺 denotes the slack variables, guaranteeing the feasibility of the DMC optimization, and 𝑤(𝑘)=[𝑤₁(𝑘),…,𝑤_𝑝(𝑘)]^𝑇 is the setpoint of controlled output obtained from upper SSO layer. 𝐐, 𝐑 are the weight coefficient matrix 𝐐=block-diag𝑄₁,…,𝑄_𝑝,𝐐_𝑖=diag𝑞_𝑖(1),…,𝑞_𝑖(𝑃)𝑅,𝑖=1,…,𝑝,𝐑=block-diag₁,…𝑅_𝑚,𝐑_𝑗=diag𝑟_𝑖(1),…,𝑟_𝑖(𝑀),𝑗=1,…,𝑚. (35)

Through the necessary conditions of extreme value 𝜕𝐽/𝜕Δ𝑢_𝑀(𝑘)=0, the optimal increment of control input can be obtained: Δ𝐮_𝑀(𝑘)=𝐀^𝑇𝐐𝐀+𝐑⁻¹𝐀^𝑇𝐐𝐰(𝑘)−𝐲_𝑃0(𝑘). (36)

The instant increment can be calculated as follows: Δ𝐮(𝑘)=𝐋𝐃𝐰(𝑘)−𝐲_𝑃0(𝑘), (37) where 𝐃=(𝐀^𝑇𝐐𝐀+𝐑)⁻¹𝐀^𝑇𝐐; remark the operation of only the first element with 𝐿=⋱10⋯00010⋯0. (38)

3.3.3. Feedback Correction

The difference between the process sample values by the present moment 𝑘 and prediction values of (30) is 𝑒()=𝑒𝑘+1¹()⋮𝑒𝑘+1_𝑝()=𝑦𝑘+1¹(𝑘+1)−𝑦_1,1()⋮𝑦𝑘+1∣𝑘_𝑝(𝑘+1)−𝑦_𝑝,1(𝑘+1∣𝑘), (39) where 𝑦_𝑖,1(𝑘+1∣𝑘) is the first element of 𝑦_{𝑖,𝑃𝑀}(𝑘+1∣𝑘), and the corrected output prediction value can be obtained using the error vector; that is, 𝐲_cor(𝑘+1)=𝐲_𝑁1(𝑘)+𝐇𝑒(𝑘+1), (40) where 𝐲_𝑁1(𝑘)=𝐲_𝑁0(𝑘)+𝐀_𝑁Δ𝐮, 𝐲_𝑁0(𝑘) is the future 𝑁 moment initial prediction value when all of the input remained unchanged at the time 𝑘; 𝐲_𝑁1(𝑘) is the future 𝑁 moment output prediction value under one-step control input action; 𝐻 is the error correct matrix. Then using a shift matrix 𝑆, next time the initial prediction value can be obtained, which is 𝑦_𝑁0(𝑘+1)=𝑆𝑦_cor(𝑘+1), (41) where 𝑆=01001⋱⋱0101_𝑁∗𝑁. (42)

4. Control System Design

In the biodiesel reactor control, multiloops are necessary to stabilize the plant. One loop is needed to maintain the set point of specifying the product purity, and another loop is needed to ensure an optimal yield of biodiesel and to minimize the generation of unwanted by-products even in the presence of disturbances.

To achieve these goals, the control loop configurations analysis is meaningful. Based on the analysis of Mjalli et al. [4], the favorable pairings are as follows: the biodiesel concentration (𝐶_𝐸) is maintained by manipulating reactant flow rate (𝐹_𝑜), the reactor temperature (𝑇) is maintained by manipulating coolant flow rate (𝐹_𝑐), respectively, and the effect of stirred rotational speed on the reactor output is insignificant, and it would be regarded as one of disturbances to the control system. The relative gain array (RGA) shows that there are some interactions among the controlled and manipulated variables which make two-layer predictive controller better qualified.

Consequently, the two-layer predictive controller is designed to handle a 2 × 2 system of inputs and outputs. The controlled output variables include biodiesel concentration (𝐶_𝐸) and reactor temperature (𝑇); the manipulated variables include reactant flow rate (𝐹_𝑜) and coolant flow rate (𝐹_𝑐). It is very important for a reactor to handle the disturbances in the feed concentration and initial temperatures, as these disturbances heavily change the system performance.

The design of the control loop based on the two-layer predictive control strategy for the biodiesel reactor is shown in Figure 3. The SSO layer searches the optimal output set-points 𝐶_𝐸𝑠𝑠 and 𝑇_𝑠𝑠 according to the economic optimization goal of the actual production process. The MPC layer selects the real-time control actions Δ𝑢 to complete the dynamic tracking control.

Figure 3 

Two-layer predictive framework of biodiesel process.

5. Simulation Results and Analysis

5.1. Model Identification

For the two-layer predictive control scheme to be successful, process modeling plays a key role in capturing the varying dynamics of the system. Section 4 shows that the biodiesel process is a two-input two-output multivariable process. The process nonlinear model was programmed and simulated in Matlab as a function. Simulation results show system is open stable process.

Firstly, generalized binary noise (GBN) signal is selected as the excitation signal. GBN signals switch between 𝑎 and −𝑎 according to the following rules: 𝑃[𝑢(𝑡)=−𝑢(𝑡−1)]=𝑝_𝑠𝑤,𝑃[𝑢(𝑡)=𝑢(𝑡−1)]=1−𝑝_𝑠𝑤, (43) where 𝑝_𝑠𝑤 is transition probability; 𝑇_min is defined as the sampling time of the signal held constant; 𝑇_𝑠𝑤 is time interval of twice conversion. The average conversion time and power spectrum are, respectively, 𝐸𝑇_𝑠𝑤=𝑇_min𝑝_𝑠𝑤,Φ_𝑢(𝜔)=1−𝑞²𝑇_min1−2𝑞cos𝑇_min𝜔+𝑞²,𝑞=1−2𝑝_𝑠𝑤. (44)

Next, least squares (LS) identification method is used to estimate the process model parameters. Suppose an MIMO plant with 𝑚 input 𝑝 output, for the 𝑖th output of the finite impulse response (FIR) model, is described as 𝑦_𝑖(𝑘)=^𝑚𝑁_𝑗=1𝑙=1ℎ_𝑖𝑗𝑙𝑢_𝑗(𝑘−𝑙). (45)

Consider experimental tests of collecting input sequence 𝑢₁(1)𝑢₁(2)⋯𝑢₁(𝐿)𝑢⋮⋮_𝑚(1)𝑢_𝑚(2)⋯𝑢_𝑚(𝐿) (46) and output sequence 𝑦₁(1)𝑦₁(2)⋯𝑦₁(𝐿)𝑦⋮⋮_𝑝(1)𝑦_𝑝(2)⋯𝑦_𝑝(𝐿). (47)

Consider matching between data and models; the introduction of residuals for each output can be independently expressed as follows: 𝑦_𝑖(𝑘)=𝜑(𝑘)𝜃_𝑖+𝑒(𝑘). (48)

Matrix form is written as 𝑦_𝑖=Φ𝜃_𝑖+𝑒, (49) where𝑦_𝑖=𝑦_𝑖()𝑦𝑁+1_𝑖()⋮𝑦𝑁+2_𝑖(𝐿)𝑒()𝑒()⋮𝑒(𝐿),𝑢,𝑒=𝑁+1𝑁+2Φ=¹(𝑁)𝑢₁(𝑁−1)⋯𝑢₁(1)𝑢_𝑚(𝑁)𝑢_𝑚(𝑁−1)⋯𝑢_𝑚(1)𝑢₁(𝑁+1)𝑢₁(𝑁)⋯𝑢₁(2)𝑢_𝑚(𝑁+1)𝑢_𝑚(𝑁)⋯𝑢_𝑚(2)𝑢⋮⋮⋮⋯⋮⋮⋮₁(𝐿−1)𝑢₁(𝐿−2)𝑢₁(𝐿−𝑁)𝑢_𝑚(𝐿−1)𝑢_𝑚(𝐿−2)𝑢_𝑚().𝐿−𝑁 (50)

Minimize the squared residuals min𝐽=𝑒^𝑇𝑒=𝑦−Φ𝜃^𝑇𝑦−Φ𝜃. (51)

Obtain the optimal estimate 𝜃=Φ^𝑇Φ⁻¹Φ^𝑇𝑦. (52)

For the model predictive controller design, the FIR model of system identification needs to be further converted into finite step response (FSR) model. The relationship between FSR coefficients and FIR coefficients is as follows: 𝑔_𝑗=^𝑗_𝑖=1ℎ_𝑗. (53)

Coefficients matrix of FSR is 𝐺^𝑢_𝑙=𝑠_11𝑙𝑠_12𝑙⋯𝑠_1𝑚𝑙𝑠_21𝑙𝑠_22𝑙⋯𝑠_2𝑚𝑙𝑠⋮⋮⋱⋮_𝑝1𝑙𝑠_𝑝2𝑙⋯𝑠_𝑝𝑚𝑙. (54)

Finally, (11)–(19) are used to create a steady-state mathematical model of two-layer prediction control. The concrete simulation process is as follows.

In the work, GBN as the excitation signal was added to the model input to produce output data. The parameters of GBN signal applied to the first input are 𝑇_𝑠𝑤=65, 𝑎𝑚𝑝=0.1 the parameters of GBN applied to the second input are 𝑇_𝑠𝑤=65, 𝑎𝑚𝑝=0.005, both the conversion probabilities are taken to be 𝑃_𝑠𝑤=1/𝑇_𝑠𝑤. Simulation time 𝑡=2000s, and sample time equals 2 s, under each input excitation, corresponding to two sets of output data each set of data capacity is 1000. Among them, the former 500 data as model identification, the remaining data are used as model validations, and FSR model length value is taken as 200.

Under the action of two inputs, reactant flow rate 𝐹_𝑜 and coolant flow rate 𝐹_𝑐, respectively, predicted value, actual value, and the relative error of two outputs biodiesel concentration 𝐶_𝐸 and reactor temperature 𝑇 were shown in Figures 4, 5, 6, and 7. Figures 4–7 show that relative error is small enough, and the model can describe 𝐶_𝐸 and 𝑇 change trends under 𝐹_𝑜 and 𝐹_𝑐.

(a)

(b)

Figure 4 

Biodiesel concentration prediction result and relative error under reactor flow rate 𝐹𝑜 action.

(a)

(b)

Figure 5 

Reactor temperature prediction result and relative error under reactor flow rate 𝐹𝑜 action.

(a)

(b)

Figure 6 

Biodiesel concentration prediction result and relative error under reactor flow rate 𝐹𝑐 action.

(a)

(b)

Figure 7 

Reactor temperature prediction result and relative error under reactor flow rate 𝐹𝑐 action.

Figures 8 and 9 give the two output step response curves under two input 𝐹_𝑜, 𝐹_𝑐 action, respectively, further shows the multiple-input multiple-output system is open-loop stable and the step response model has been identified successfully. The FSR model will be utilized to represent the actual process in latter optimization and controller design.

(a)

(b)

Figure 8 

Step response curve of biodiesel concentration and reactor temperature, respectively, under 𝐹𝑜 action.

(a)

(b)

Figure 9 

Step response curve of biodiesel concentration and reactor temperature, respectively, under 𝐹𝑐 action.

5.2. Dynamic Simulation

To validate the effectiveness and immunity in two-layer predictive control, the models obtained in Section 5.1 are used in the simulations.

The reaction rate constants come from [18] under the common industrial conditions of 6 : 1 methanol/oil mole ratio, 1.0 wt% catalyst KOH, and 600 rpm stirrer rotational speed. These kinetics parameters can be considered as constants. The initial operating conditions refer to the literature [4] the validated data. According to these parameters and reaction conditions, the simulation of biodiesel transesterification reactor can be carried out.

The economic optimization method described in (22) is adopted as SSO whose main parameters are selected as follows: the cost coefficients of control input in steady-state optimization are set to [1;−1], the input 𝐹_𝑜 is constrained between 0 and 0.2 m³/s, the input 𝐹_𝑐 is constrained between 0 and 0.1 m³/s, and the output 𝐶_𝐸 is constrained between 3.0536 kmol/m³ and 3.196 kmol/m³, the output 𝑇 is constrained between 337.77 K and 338.25 K.

The parameters of the dynamic control layer adopted the unconstrained DMC algorithm: the modeling time domain 𝑁=200, prediction horizon 𝑃=200, control horizon 𝑀=20. The weight coefficient values of weight matrix 𝑄 and 𝑅 equal to 10 and 1000, respectively.

Conventional PID controller has also been designed in this simulation for comparison of performance to two-layer predictive controller. The parameters of PID controller for 𝐶_𝐸 with 𝐹_𝑜 control loop are 𝑘𝑝=−6𝑒−5, 𝑘𝑖=−0.05, and 𝑘𝑑=0; the parameters for 𝑇 with 𝐹_𝑐 control loop are 𝑘𝑝=−0.02, 𝑘𝑖=−0.001, and 𝑘𝑑=0. The simulations of general PID controller and two-layer predictive controller are compared to validate the performance of the latter algorithm, whose results are shown in Figures 10 and 11.

(a)

(b)

Figure 10 

Biodiesel concentration and controller moves of two-layer predictive controller and PID controller.

(a)

(b)

Figure 11 

Reactor temperature and controller moves of two-layer predictive controller and PID controller.

As Figures 10 and 11 show, the two-layer predictive controller starts running at the time 𝑡=0. The results of steady state optimization are 𝑦_𝑠𝑠=[3.196,337.77],𝑢_𝑠𝑠=[0.073,0.0062]. (55)

The optimized values as the setpoints were send to the lower layer DMC. In the beginning, the closed loop response of the two-layer predictive controller was a little sluggish in bringing the biodiesel concentration back the optimum steady-state values, this is because that the algorithm enter the constraint adjustment stage based on the priority strategy which adjusting the upper limit and lower limit to be handled. About At the time 𝑡=400, the response gradually becomes stable. It can be seen that the two-layer predictive controller preceded the PID controller in terms of the ability to attain lower overshoot, smaller oscillation, and faster response time.

Considering the actual application, the control input is also an important indicator of good or bad controller. From Figures 10 and 11, the two-layer predictive controller has much more stable controller moves than does PID that meets the practical implementation constrains.

To challenge the stability of two-layer predictive controller, some disturbances were exerted alone and at the same time. The chosen disturbance variables include coolant input temperature (𝑇_𝑐0), feed temperature (𝑇₀), triglyceride initial concentration (𝐶_TG0), and stirrer rotational speed (𝑁). After the system has attained the steady state, The nominal values of 𝑇_𝑐0, 𝑇₀ were increased 3 K, respectively, and 𝐶_TG0, 𝑁 were increased 5%, respectively, at the time 𝑡=1000s. Figures 12 and 13 show the biodiesel concentration and reactor temperature profiles when these disturbance variables were introduced.

(a)

(b)

Figure 12 

Biodiesel concentration and controller moves of four individual disturbance variables effects.

(a)

(b)

Figure 13 

Reactor temperature and controller moves of four individual disturbance variables effects.

Figures 12 and 13 showed satisfactory rejection of all disturbances. Two-layer predictive controller was able to bring back the controlled variables to their setpoints in less than 1000 s, and overshoot was within the acceptable range. For the biodiesel concentration loop, the initial concentration 𝐶_TG0 has the highest effect, with an overshoot of less than 0.01 kmol/m³. For the reactor temperature loop, the feed temperature 𝑇₀ has the largest effect, with an overshoot of less than 0.33 K. For the two loops, the stirrer rotational speed almost has no effect on the controlled variables.

6. Conclusions

Biodiesel transesterification reactor control has become very important in recent years due to the difficulty in controlling the complex and highly nonlinear dynamic behavior. In this paper, a novel two-layer predictive control scheme for a continuous biodiesel transesterification reactor has been proposed. The SSO layer achieved optimal output setpoints according to the local economic optimization goal of the actual production process, and the MPC layer realized the dynamic tracking control. The main aim was to optimize and control the biodiesel concentration and reactor temperature in order to obtain the product of the highest quality at the lower cost. With steady-state optimum target calculation and DMC algorithm implement, the performance of the two-layer predictive controller was superior to that of a conventional PID controller. The two-layer predictive control is not only stable but also tracks set points more efficiently with minimal overshoots and shorter settling times. Moreover, it exhibits good disturbance rejection characteristics.

Acknowledgments

This work is supported by the National Natural Science Foundation of China (61034008) and the Science Research Foundation of Liaoning Provincial Department of Education (L2012145).