# Two-Layer Predictive Control of a Continuous Biodiesel Transesterification Reactor

Authors: Hongyan Shi; Dingding Wang; Decheng Yuan; Tianran Wang
Journal: Journal of Applied Mathematics
 (2013)

---

## 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.


---

## Body

## 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. Section2 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
(1)TG+3A⟷3E+G.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]:
(3)dCTGdt=-k1′CTGCA+k2′CDGCE,dCDGdt=k1′CTGCA-k2′CDGCE-k3′CDGCA+k4′CMGCE,dCMGdt=k3′CDGCA-k4′CMGCE-k5′CMGCA+k6′CGLCE,dCEdt=k1′CTGCA-k2′CDGCE+k3′CDGCA-k4′CMGCE+k5′CMGCA-k6′CGLCE,dCAdt=-dCEdt,dCGLdt=k5′CMGCA-k6′CGLCE,
where CTG, CDG, CMG, CE, CA, and CGL are concentrations of triglyceride, diglyceride, monoglyceride, methyl ester, methanol, and glycerol, respectively. k1′, k3′, and k5′ are the effective rate constants for the forward reactions, and k2′, k4′, and k6′ 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)rj=kj′[Ci]2.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 (Ccat) and the corresponding rate constant for the catalyzed reaction [18]:
(5)kj′=kjCcat.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)kj=k0e(-Ea/RT),
where k0 is a constant called the preexponential factor, Ea is the activation energy of the reaction, and R 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]:
(7)VdCidt=Fi0Ci0-FiCi-∑j=1nrjV.The reactor energy balance is expressed as(8)V∑i=1niCiCpidTdt=(FTG0CTG0CPTG+FA0CA0CPA)(T0-T)-(V∑i=jnrjΔHj)-(UAHΔT).The coolant fluid energy balance is expressed as(9)dTCdt=FC0VC(TC0-TC)+UAHΔTρCVCCPC.The function equation of heat transfer coefficient is approximately expressed as(10)U=αFCβNγ=735.5FC1.095N0.405.


## 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 withm control input and p controlled output and the coefficients of the corresponding step response model between control input uj and output yi are given; the model vector is
(11)aij(t)=[aij(1),…,aij(N)]T,
where i=1,…,p;  j=1,…,m. N in (11) denotes modeling horizon of step response model. Thus, a multistep predictive model can be obtained:
(12)y(k+1)=y(k)+A1Δu(k),
where
(13)y(k+1)=[y1(k+1)⋮yp(k+1)];y(k)=[y1(k)⋮yp(k)];Δu(k)=[Δu1(k)⋮Δum(k)];A1=[a11(1)⋯a1m(1)⋮⋮⋮ap1(1)⋯apm(1)].Under the control incrementΔu(k),…,Δu(k+M-1) action, the output predictive value of the system is
(14)y(k+1)=y(k)+A1Δu(k),y(k+2)=y(k)+A2Δu(k)+A1Δu(k+1),⋮y(k+N)=y(k)+ANΔu(k)+⋯+AN-M+1Δu(k+M-1),abbreviated as(15)∂y(k)=AΔuM(k),
where
(16)∂y(k)=[y(k+1)-y(k)⋮y(k+N)-y(k)],ΔuM(k)=[Δu(k)⋮Δu(k+M-1)],A=[A10⋮⋱AM⋯A1⋮⋮AN⋯AN-M+1].The system can be written at the steady-state time(17)Δy(∞)=ANΔu(∞),
where Δy(∞)=[Δy1(∞),Δy2(∞),…,Δyp(∞)]T, Δu(∞)=[Δu1(∞),Δu2(∞),…,Δum(∞)]T are the steady-state output increment and input increment, respectively, and AN is the steady-state step response coefficients matrix
(18)AN=[a11(N)⋯a1m(N)⋮⋮⋮ap1(N)⋯apm(N)].To meet the requirements of steady-state target calculation, model (17) can also be written as
(19)Δy∞(k)=ANΔu∞(k).


### 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]:
(20)minΔu∞(k),Δy∞(k)J=αTΔu∞(k)+βTΔy∞(k).SinceΔu∞ and Δy∞ 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
(21)minΔu∞(k)J=cTΔu∞(k),
where cT=[c1,…,cm] is the cost coefficient vector, constructed by the normalized benefit, or cost of each input variable. Δu∞(k)=[Δu∞1,…,Δu∞m]T is the steady-state change value of every input at time k.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:
(22)minΔu∞(k)J=cTΔu∞(k)s.t.Δy∞(k)=GuΔu∞(k)+GfΔf∞(k)+e,s.t.umin≤u∞(k)+Δu∞(k)≤umax,s.t.ymin≤y∞(k)+Δy∞(k)≤ymax,
where Gu, Gf are the steady-state gain matrices of control input and disturbance variables; and e is the model bias. umin, umax are low limit and upper limit of steady-state input variables ymin, ymax 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:
(23)minΔu∞(k)J=cT(Δu∞(k)-Maxprofit)2s.t.Δy∞(k)=GuΔu∞(k)+GfΔf∞(k)+e,s.t.umin≤u∞(k)+Δu∞(k)≤umax,s.t.ymin≤y∞(k)+Δy∞(k)≤ymax,
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 “N”, where
(24)Δy∞(k)=GuΔu∞(k)+GfΔf∞(k)+e,uLL≤u∞(k)+Δu∞(k)≤uHL,yLLj-ε2j≤y∞(k)+Δy∞(k)≤yHLj+ε1j,ε1j≥0,ε2j≥0,ε1j≤yHHL-yHL,ε2j≤yLL-yLLL,j=1,…,N.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 Step2, constraints of the highest priority rank N are adjusted first by solving the following optimization problem:
(25)minεNJ=(WN)TεN,(WN)T=[W1N,…,W2×nNN]s.t.ΘNZN=bN,s.t.ΩNZN≤ΨN,
where(26)ZN=[X4N)T,(ε1N)T,(ε2N)T,(ε1N)T,(ε2N)T]X1T,X2T,(X31)T,…,(X3N)T,(X41)T,…,(X4N)T,(ε1N)T,(ε2N)T,(ε1N)T,(ε2N)T]T,ΩN=block-diag(-Im,-Im,-In1,…,-InN,=block-diaghh-In1,…,-InN,-InN,-InN,InN,InN),ΨN=[(0n1×1)T(0m×1)T,(0m×1)T,(0n1×1)T,…,(0nN×1)T,(0n1×1)T,…,(0nN×1)T,(0nN×1)T,(0nN×1)T,(yHHLN-yHLN)T,(yLLN-yLLLN)T(0n1×1)T]T,bN=[uHL-uLLGu1u∞(k)-Gu1uLL(k)+yHL1-y∞1(k)-Gf1Δf∞(k)-e1⋮GuNu∞(k)-GuNuLL(k)+yHLN-y∞N(k)-GfNΔf∞(k)-eNGu1uHL-Gu1U∞(k)+Y∞1(k)+Gf1Δf∞(k)-yLL1+e1⋮GuNuHL-GuNU∞(k)+y∞N(k)+GfNΔf∞(k)-yLLN+eN00],ΘN=[ImIm0⋯00⋯00000Gu100⋯0In1000000⋮⋮0⋯00⋱00000GuN00⋯000InN-InN0000Gu1In1000⋯00000⋮⋮0⋱00⋯000000GuN00InN0⋯00-InN00000⋯00⋯0InN0-InN0000⋯00⋯00InN0-InN].Solving (25) may appear in three different cases, respectively: if (25) is feasible, and the optimum solution is εN=0, subject to J=0, that is, no need for soft constraints adjustment, directly solve the original problem (22); if (25) is feasible, but εN≠0, just need to relax constraints of priority ranks N, and further optimization solution; if (25) is infeasible, not get a feasible solution to soft constraints adjustment of the priority rank N, relaxing the constraints of the priority rank N to hard constraints; that is,
(27)ε1N=yHHLN-yHLN,ε2N=yLLN-yLLLN.Go to the procedure of judging rankN-1 constraints
(28)minεN-1J=(WN-1)TεN-1,(WN-1)T=[W1N-1,…,W2×nN-1N-1]s.t.ΘN-1ZN-1=bN-1,ΩN-1ZN-1≤ΨN-1.For (28), the matrix form is the same with priority rank N, only in the corresponding position of εN-1   to replace εN, bN-1 matrix is adjusted(29)bN-1=[uHL-uLLGu1u∞(k)-Gu1uLL(k)+yHL1-y∞1(k)-Gf1Δf∞(k)-e1⋮GuNu∞(k)-GuNuLL(k)+yHLN-y∞N(k)-GfNΔf∞(k)-eN+(yHHLN-yHLN)TGu1uHL-Gu1u∞(k)+y∞1(k)+Gf1Δf∞(k)-yLL1+e1⋮GuNuHL-GuNu∞(k)+y∞N(k)+GfNΔf∞(k)-yLLN+eN+(yLLN-yLLLN)T00].N
-
1 rank and N 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 (us,ys) 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 timek, the future P-step prediction output can be written as follows:
(30)y~PM(k)=y~P0(k)+AΔuM(k),
where P denotes the prediction horizon, M is the control horizon, A is the prediction matrix composed by the corresponding step response coefficients, y~P0 is the initial output prediction value when control action starting from the present time does not change, ΔuM(k) is the prediction incremental in M control horizon, and y~PM(k) is the future P-step prediction output under M-step control action change. Among them
(31)y~PM(k)=[y~1,PM(k)⋮y~p,PM(k)],y~P0(k)=[y~1,P0(k)⋮y~p,P0(k)],ΔuM(k)=[Δu1,M(k)⋮Δum,M(k)],A=[A11⋯A1m⋮⋱⋮Ap1⋯Apm].


#### 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:(32)minΔuM(k)J(k)=∥w(k)-y~PM(k)∥Q2+∥ε(k)∥S2+∥uM(k)-u∞∥T2+∥ΔuM(k)∥R2.Subject to the model(33)y~PM(k)=y~P0(k)+AΔuM(k).Subject to bound constraints(34)ymin-ε≤y~PM(k)≤ymax+ε,umin≤uM≤umax,Δumin≤ΔuM(k)≤Δumax,
where ε denotes the slack variables, guaranteeing the feasibility of the DMC optimization, and w(k)=[w1(k),…,wp(k)]T is the setpoint of controlled output obtained from upper SSO layer. Q, R are the weight coefficient matrix
(35)Q=block-diag(Q1,…,Qp),Qi=diag(qi(1),…,qi(P)),i=1,…,p,R=block-diag(R1,…Rm),Rj=diag(ri(1),…,ri(M)),j=1,…,m.Through the necessary conditions of extreme value∂J/∂ΔuM(k)=0, the optimal increment of control input can be obtained:
(36)ΔuM(k)=(ATQA+R)-1ATQ[w(k)-y~P0(k)].The instant increment can be calculated as follows:(37)Δu(k)=LD[w(k)-y~P0(k)],
where D=(ATQA+R)-1ATQ; remark the operation of only the first element with
(38)L=[10⋯00⋱010⋯0].


#### 3.3.3. Feedback Correction

The difference between the process sample values by the present momentk and prediction values of (30) is
(39)e(k+1)=[e1(k+1)⋮ep(k+1)]=[y1(k+1)-y~1,1(k+1∣k)⋮yp(k+1)-y~p,1(k+1∣k)],
where y~i,1(k+1∣k) is the first element of y~i,PM(k+1∣k), and the corrected output prediction value can be obtained using the error vector; that is,
(40)y~cor(k+1)=y~N1(k)+He(k+1),
where y~N1(k)=y~N0(k)+ANΔu, y~N0(k) is the future N moment initial prediction value when all of the input remained unchanged at the time k; y~N1(k) is the future N moment output prediction value under one-step control input action; H is the error correct matrix. Then using a shift matrix S, next time the initial prediction value can be obtained, which is
(41)y~N0(k+1)=Sy~cor(k+1),
where
(42)S=[01001⋱⋱0101]N*N.


## 3.1. Establish Steady-State Mathematical Model

Assume an MIMO plant withm control input and p controlled output and the coefficients of the corresponding step response model between control input uj and output yi are given; the model vector is
(11)aij(t)=[aij(1),…,aij(N)]T,
where i=1,…,p;  j=1,…,m. N in (11) denotes modeling horizon of step response model. Thus, a multistep predictive model can be obtained:
(12)y(k+1)=y(k)+A1Δu(k),
where
(13)y(k+1)=[y1(k+1)⋮yp(k+1)];y(k)=[y1(k)⋮yp(k)];Δu(k)=[Δu1(k)⋮Δum(k)];A1=[a11(1)⋯a1m(1)⋮⋮⋮ap1(1)⋯apm(1)].Under the control incrementΔu(k),…,Δu(k+M-1) action, the output predictive value of the system is
(14)y(k+1)=y(k)+A1Δu(k),y(k+2)=y(k)+A2Δu(k)+A1Δu(k+1),⋮y(k+N)=y(k)+ANΔu(k)+⋯+AN-M+1Δu(k+M-1),abbreviated as(15)∂y(k)=AΔuM(k),
where
(16)∂y(k)=[y(k+1)-y(k)⋮y(k+N)-y(k)],ΔuM(k)=[Δu(k)⋮Δu(k+M-1)],A=[A10⋮⋱AM⋯A1⋮⋮AN⋯AN-M+1].The system can be written at the steady-state time(17)Δy(∞)=ANΔu(∞),
where Δy(∞)=[Δy1(∞),Δy2(∞),…,Δyp(∞)]T, Δu(∞)=[Δu1(∞),Δu2(∞),…,Δum(∞)]T are the steady-state output increment and input increment, respectively, and AN is the steady-state step response coefficients matrix
(18)AN=[a11(N)⋯a1m(N)⋮⋮⋮ap1(N)⋯apm(N)].To meet the requirements of steady-state target calculation, model (17) can also be written as
(19)Δy∞(k)=ANΔu∞(k).


## 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]:
(20)minΔu∞(k),Δy∞(k)J=αTΔu∞(k)+βTΔy∞(k).SinceΔu∞ and Δy∞ 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
(21)minΔu∞(k)J=cTΔu∞(k),
where cT=[c1,…,cm] is the cost coefficient vector, constructed by the normalized benefit, or cost of each input variable. Δu∞(k)=[Δu∞1,…,Δu∞m]T is the steady-state change value of every input at time k.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:
(22)minΔu∞(k)J=cTΔu∞(k)s.t.Δy∞(k)=GuΔu∞(k)+GfΔf∞(k)+e,s.t.umin≤u∞(k)+Δu∞(k)≤umax,s.t.ymin≤y∞(k)+Δy∞(k)≤ymax,
where Gu, Gf are the steady-state gain matrices of control input and disturbance variables; and e is the model bias. umin, umax are low limit and upper limit of steady-state input variables ymin, ymax 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:
(23)minΔu∞(k)J=cT(Δu∞(k)-Maxprofit)2s.t.Δy∞(k)=GuΔu∞(k)+GfΔf∞(k)+e,s.t.umin≤u∞(k)+Δu∞(k)≤umax,s.t.ymin≤y∞(k)+Δy∞(k)≤ymax,
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 “N”, where
(24)Δy∞(k)=GuΔu∞(k)+GfΔf∞(k)+e,uLL≤u∞(k)+Δu∞(k)≤uHL,yLLj-ε2j≤y∞(k)+Δy∞(k)≤yHLj+ε1j,ε1j≥0,ε2j≥0,ε1j≤yHHL-yHL,ε2j≤yLL-yLLL,j=1,…,N.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 Step2, constraints of the highest priority rank N are adjusted first by solving the following optimization problem:
(25)minεNJ=(WN)TεN,(WN)T=[W1N,…,W2×nNN]s.t.ΘNZN=bN,s.t.ΩNZN≤ΨN,
where(26)ZN=[X4N)T,(ε1N)T,(ε2N)T,(ε1N)T,(ε2N)T]X1T,X2T,(X31)T,…,(X3N)T,(X41)T,…,(X4N)T,(ε1N)T,(ε2N)T,(ε1N)T,(ε2N)T]T,ΩN=block-diag(-Im,-Im,-In1,…,-InN,=block-diaghh-In1,…,-InN,-InN,-InN,InN,InN),ΨN=[(0n1×1)T(0m×1)T,(0m×1)T,(0n1×1)T,…,(0nN×1)T,(0n1×1)T,…,(0nN×1)T,(0nN×1)T,(0nN×1)T,(yHHLN-yHLN)T,(yLLN-yLLLN)T(0n1×1)T]T,bN=[uHL-uLLGu1u∞(k)-Gu1uLL(k)+yHL1-y∞1(k)-Gf1Δf∞(k)-e1⋮GuNu∞(k)-GuNuLL(k)+yHLN-y∞N(k)-GfNΔf∞(k)-eNGu1uHL-Gu1U∞(k)+Y∞1(k)+Gf1Δf∞(k)-yLL1+e1⋮GuNuHL-GuNU∞(k)+y∞N(k)+GfNΔf∞(k)-yLLN+eN00],ΘN=[ImIm0⋯00⋯00000Gu100⋯0In1000000⋮⋮0⋯00⋱00000GuN00⋯000InN-InN0000Gu1In1000⋯00000⋮⋮0⋱00⋯000000GuN00InN0⋯00-InN00000⋯00⋯0InN0-InN0000⋯00⋯00InN0-InN].Solving (25) may appear in three different cases, respectively: if (25) is feasible, and the optimum solution is εN=0, subject to J=0, that is, no need for soft constraints adjustment, directly solve the original problem (22); if (25) is feasible, but εN≠0, just need to relax constraints of priority ranks N, and further optimization solution; if (25) is infeasible, not get a feasible solution to soft constraints adjustment of the priority rank N, relaxing the constraints of the priority rank N to hard constraints; that is,
(27)ε1N=yHHLN-yHLN,ε2N=yLLN-yLLLN.Go to the procedure of judging rankN-1 constraints
(28)minεN-1J=(WN-1)TεN-1,(WN-1)T=[W1N-1,…,W2×nN-1N-1]s.t.ΘN-1ZN-1=bN-1,ΩN-1ZN-1≤ΨN-1.For (28), the matrix form is the same with priority rank N, only in the corresponding position of εN-1   to replace εN, bN-1 matrix is adjusted(29)bN-1=[uHL-uLLGu1u∞(k)-Gu1uLL(k)+yHL1-y∞1(k)-Gf1Δf∞(k)-e1⋮GuNu∞(k)-GuNuLL(k)+yHLN-y∞N(k)-GfNΔf∞(k)-eN+(yHHLN-yHLN)TGu1uHL-Gu1u∞(k)+y∞1(k)+Gf1Δf∞(k)-yLL1+e1⋮GuNuHL-GuNu∞(k)+y∞N(k)+GfNΔf∞(k)-yLLN+eN+(yLLN-yLLLN)T00].N
-
1 rank and N 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.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]:
(20)minΔu∞(k),Δy∞(k)J=αTΔu∞(k)+βTΔy∞(k).SinceΔu∞ and Δy∞ 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
(21)minΔu∞(k)J=cTΔu∞(k),
where cT=[c1,…,cm] is the cost coefficient vector, constructed by the normalized benefit, or cost of each input variable. Δu∞(k)=[Δu∞1,…,Δu∞m]T is the steady-state change value of every input at time k.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:
(22)minΔu∞(k)J=cTΔu∞(k)s.t.Δy∞(k)=GuΔu∞(k)+GfΔf∞(k)+e,s.t.umin≤u∞(k)+Δu∞(k)≤umax,s.t.ymin≤y∞(k)+Δy∞(k)≤ymax,
where Gu, Gf are the steady-state gain matrices of control input and disturbance variables; and e is the model bias. umin, umax are low limit and upper limit of steady-state input variables ymin, ymax 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:
(23)minΔu∞(k)J=cT(Δu∞(k)-Maxprofit)2s.t.Δy∞(k)=GuΔu∞(k)+GfΔf∞(k)+e,s.t.umin≤u∞(k)+Δu∞(k)≤umax,s.t.ymin≤y∞(k)+Δy∞(k)≤ymax,
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 “N”, where
(24)Δy∞(k)=GuΔu∞(k)+GfΔf∞(k)+e,uLL≤u∞(k)+Δu∞(k)≤uHL,yLLj-ε2j≤y∞(k)+Δy∞(k)≤yHLj+ε1j,ε1j≥0,ε2j≥0,ε1j≤yHHL-yHL,ε2j≤yLL-yLLL,j=1,…,N.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 Step2, constraints of the highest priority rank N are adjusted first by solving the following optimization problem:
(25)minεNJ=(WN)TεN,(WN)T=[W1N,…,W2×nNN]s.t.ΘNZN=bN,s.t.ΩNZN≤ΨN,
where(26)ZN=[X4N)T,(ε1N)T,(ε2N)T,(ε1N)T,(ε2N)T]X1T,X2T,(X31)T,…,(X3N)T,(X41)T,…,(X4N)T,(ε1N)T,(ε2N)T,(ε1N)T,(ε2N)T]T,ΩN=block-diag(-Im,-Im,-In1,…,-InN,=block-diaghh-In1,…,-InN,-InN,-InN,InN,InN),ΨN=[(0n1×1)T(0m×1)T,(0m×1)T,(0n1×1)T,…,(0nN×1)T,(0n1×1)T,…,(0nN×1)T,(0nN×1)T,(0nN×1)T,(yHHLN-yHLN)T,(yLLN-yLLLN)T(0n1×1)T]T,bN=[uHL-uLLGu1u∞(k)-Gu1uLL(k)+yHL1-y∞1(k)-Gf1Δf∞(k)-e1⋮GuNu∞(k)-GuNuLL(k)+yHLN-y∞N(k)-GfNΔf∞(k)-eNGu1uHL-Gu1U∞(k)+Y∞1(k)+Gf1Δf∞(k)-yLL1+e1⋮GuNuHL-GuNU∞(k)+y∞N(k)+GfNΔf∞(k)-yLLN+eN00],ΘN=[ImIm0⋯00⋯00000Gu100⋯0In1000000⋮⋮0⋯00⋱00000GuN00⋯000InN-InN0000Gu1In1000⋯00000⋮⋮0⋱00⋯000000GuN00InN0⋯00-InN00000⋯00⋯0InN0-InN0000⋯00⋯00InN0-InN].Solving (25) may appear in three different cases, respectively: if (25) is feasible, and the optimum solution is εN=0, subject to J=0, that is, no need for soft constraints adjustment, directly solve the original problem (22); if (25) is feasible, but εN≠0, just need to relax constraints of priority ranks N, and further optimization solution; if (25) is infeasible, not get a feasible solution to soft constraints adjustment of the priority rank N, relaxing the constraints of the priority rank N to hard constraints; that is,
(27)ε1N=yHHLN-yHLN,ε2N=yLLN-yLLLN.Go to the procedure of judging rankN-1 constraints
(28)minεN-1J=(WN-1)TεN-1,(WN-1)T=[W1N-1,…,W2×nN-1N-1]s.t.ΘN-1ZN-1=bN-1,ΩN-1ZN-1≤ΨN-1.For (28), the matrix form is the same with priority rank N, only in the corresponding position of εN-1   to replace εN, bN-1 matrix is adjusted(29)bN-1=[uHL-uLLGu1u∞(k)-Gu1uLL(k)+yHL1-y∞1(k)-Gf1Δf∞(k)-e1⋮GuNu∞(k)-GuNuLL(k)+yHLN-y∞N(k)-GfNΔf∞(k)-eN+(yHHLN-yHLN)TGu1uHL-Gu1u∞(k)+y∞1(k)+Gf1Δf∞(k)-yLL1+e1⋮GuNuHL-GuNu∞(k)+y∞N(k)+GfNΔf∞(k)-yLLN+eN+(yLLN-yLLLN)T00].N
-
1 rank and N 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 (us,ys) 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 timek, the future P-step prediction output can be written as follows:
(30)y~PM(k)=y~P0(k)+AΔuM(k),
where P denotes the prediction horizon, M is the control horizon, A is the prediction matrix composed by the corresponding step response coefficients, y~P0 is the initial output prediction value when control action starting from the present time does not change, ΔuM(k) is the prediction incremental in M control horizon, and y~PM(k) is the future P-step prediction output under M-step control action change. Among them
(31)y~PM(k)=[y~1,PM(k)⋮y~p,PM(k)],y~P0(k)=[y~1,P0(k)⋮y~p,P0(k)],ΔuM(k)=[Δu1,M(k)⋮Δum,M(k)],A=[A11⋯A1m⋮⋱⋮Ap1⋯Apm].


### 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:(32)minΔuM(k)J(k)=∥w(k)-y~PM(k)∥Q2+∥ε(k)∥S2+∥uM(k)-u∞∥T2+∥ΔuM(k)∥R2.Subject to the model(33)y~PM(k)=y~P0(k)+AΔuM(k).Subject to bound constraints(34)ymin-ε≤y~PM(k)≤ymax+ε,umin≤uM≤umax,Δumin≤ΔuM(k)≤Δumax,
where ε denotes the slack variables, guaranteeing the feasibility of the DMC optimization, and w(k)=[w1(k),…,wp(k)]T is the setpoint of controlled output obtained from upper SSO layer. Q, R are the weight coefficient matrix
(35)Q=block-diag(Q1,…,Qp),Qi=diag(qi(1),…,qi(P)),i=1,…,p,R=block-diag(R1,…Rm),Rj=diag(ri(1),…,ri(M)),j=1,…,m.Through the necessary conditions of extreme value∂J/∂ΔuM(k)=0, the optimal increment of control input can be obtained:
(36)ΔuM(k)=(ATQA+R)-1ATQ[w(k)-y~P0(k)].The instant increment can be calculated as follows:(37)Δu(k)=LD[w(k)-y~P0(k)],
where D=(ATQA+R)-1ATQ; remark the operation of only the first element with
(38)L=[10⋯00⋱010⋯0].


### 3.3.3. Feedback Correction

The difference between the process sample values by the present momentk and prediction values of (30) is
(39)e(k+1)=[e1(k+1)⋮ep(k+1)]=[y1(k+1)-y~1,1(k+1∣k)⋮yp(k+1)-y~p,1(k+1∣k)],
where y~i,1(k+1∣k) is the first element of y~i,PM(k+1∣k), and the corrected output prediction value can be obtained using the error vector; that is,
(40)y~cor(k+1)=y~N1(k)+He(k+1),
where y~N1(k)=y~N0(k)+ANΔu, y~N0(k) is the future N moment initial prediction value when all of the input remained unchanged at the time k; y~N1(k) is the future N moment output prediction value under one-step control input action; H is the error correct matrix. Then using a shift matrix S, next time the initial prediction value can be obtained, which is
(41)y~N0(k+1)=Sy~cor(k+1),
where
(42)S=[01001⋱⋱0101]N*N.


## 3.3.1. Predictive Model

Based on system process step response model, at the current timek, the future P-step prediction output can be written as follows:
(30)y~PM(k)=y~P0(k)+AΔuM(k),
where P denotes the prediction horizon, M is the control horizon, A is the prediction matrix composed by the corresponding step response coefficients, y~P0 is the initial output prediction value when control action starting from the present time does not change, ΔuM(k) is the prediction incremental in M control horizon, and y~PM(k) is the future P-step prediction output under M-step control action change. Among them
(31)y~PM(k)=[y~1,PM(k)⋮y~p,PM(k)],y~P0(k)=[y~1,P0(k)⋮y~p,P0(k)],ΔuM(k)=[Δu1,M(k)⋮Δum,M(k)],A=[A11⋯A1m⋮⋱⋮Ap1⋯Apm].


## 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:(32)minΔuM(k)J(k)=∥w(k)-y~PM(k)∥Q2+∥ε(k)∥S2+∥uM(k)-u∞∥T2+∥ΔuM(k)∥R2.Subject to the model(33)y~PM(k)=y~P0(k)+AΔuM(k).Subject to bound constraints(34)ymin-ε≤y~PM(k)≤ymax+ε,umin≤uM≤umax,Δumin≤ΔuM(k)≤Δumax,
where ε denotes the slack variables, guaranteeing the feasibility of the DMC optimization, and w(k)=[w1(k),…,wp(k)]T is the setpoint of controlled output obtained from upper SSO layer. Q, R are the weight coefficient matrix
(35)Q=block-diag(Q1,…,Qp),Qi=diag(qi(1),…,qi(P)),i=1,…,p,R=block-diag(R1,…Rm),Rj=diag(ri(1),…,ri(M)),j=1,…,m.Through the necessary conditions of extreme value∂J/∂ΔuM(k)=0, the optimal increment of control input can be obtained:
(36)ΔuM(k)=(ATQA+R)-1ATQ[w(k)-y~P0(k)].The instant increment can be calculated as follows:(37)Δu(k)=LD[w(k)-y~P0(k)],
where D=(ATQA+R)-1ATQ; remark the operation of only the first element with
(38)L=[10⋯00⋱010⋯0].


## 3.3.3. Feedback Correction

The difference between the process sample values by the present momentk and prediction values of (30) is
(39)e(k+1)=[e1(k+1)⋮ep(k+1)]=[y1(k+1)-y~1,1(k+1∣k)⋮yp(k+1)-y~p,1(k+1∣k)],
where y~i,1(k+1∣k) is the first element of y~i,PM(k+1∣k), and the corrected output prediction value can be obtained using the error vector; that is,
(40)y~cor(k+1)=y~N1(k)+He(k+1),
where y~N1(k)=y~N0(k)+ANΔu, y~N0(k) is the future N moment initial prediction value when all of the input remained unchanged at the time k; y~N1(k) is the future N moment output prediction value under one-step control input action; H is the error correct matrix. Then using a shift matrix S, next time the initial prediction value can be obtained, which is
(41)y~N0(k+1)=Sy~cor(k+1),
where
(42)S=[01001⋱⋱0101]N*N.


## 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 (CE) is maintained by manipulating reactant flow rate (Fo), the reactor temperature (T) is maintained by manipulating coolant flow rate (Fc), 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 (CE) and reactor temperature (T); the manipulated variables include reactant flow rate (Fo) and coolant flow rate (Fc). 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 Figure3. The SSO layer searches the optimal output set-points CEss and Tss according to the economic optimization goal of the actual production process. The MPC layer selects the real-time control actions Δu 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. Section4 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 betweena and -a according to the following rules:
(43)P[u(t)=-u(t-1)]=psw,P[u(t)=u(t-1)]=1-psw,
where psw is transition probability; Tmin is defined as the sampling time of the signal held constant; Tsw is time interval of twice conversion. The average conversion time and power spectrum are, respectively,
(44)ETsw=Tminpsw,Φu(ω)=(1-q2)Tmin1-2qcosTminω+q2,q=1-2psw.Next, least squares (LS) identification method is used to estimate the process model parameters. Suppose an MIMO plant with m input p output, for the ith output of the finite impulse response (FIR) model, is described as
(45)yi(k)=∑j=1m∑l=1Nhijluj(k-l).Consider experimental tests of collecting input sequence(46)u1(1)u1(2)⋯u1(L)⋮⋮um(1)um(2)⋯um(L)
and output sequence
(47)y1(1)y1(2)⋯y1(L)⋮⋮yp(1)yp(2)⋯yp(L).Consider matching between data and models; the introduction of residuals for each output can be independently expressed as follows:(48)yi(k)=φ(k)θi+e(k).Matrix form is written as(49)yi=Φθi+e,
where(50)yi=[yi(N+1)yi(N+2)⋮yi(L)],e=[e(N+1)e(N+2)⋮e(L)],Φ=[u1(N)u1(N-1)⋯u1(1)um(N)um(N-1)⋯um(1)u1(N+1)u1(N)⋯u1(2)um(N+1)um(N)⋯um(2)⋮⋮⋮⋯⋮⋮⋮u1(L-1)u1(L-2)u1(L-N)um(L-1)um(L-2)um(L-N)].Minimize the squared residuals(51)minJ=eTe=[y-Φθ]T[y-Φθ].Obtain the optimal estimate(52)θ^=[ΦTΦ]-1ΦTy.For the model predictive controller design, theFIR 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:
(53)gj=∑i=1jhj.Coefficients matrix ofFSR is
(54)Glu=[s11ls12l⋯s1mls21ls22l⋯s2ml⋮⋮⋱⋮sp1lsp2l⋯spml].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 areTsw=65, amp=0.1 the parameters of GBN applied to the second input are Tsw=65, amp=0.005, both the conversion probabilities are taken to be Psw=1/Tsw. Simulation time t=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 rateFo and coolant flow rate Fc, respectively, predicted value, actual value, and the relative error of two outputs biodiesel concentration CE and reactor temperature T were shown in Figures 4, 5, 6, and 7. Figures 4–7 show that relative error is small enough, and the model can describe CE and T change trends under Fo and Fc.Biodiesel concentration prediction result and relative error under reactor flow rateFo action.
(a)
(b)Reactor temperature prediction result and relative error under reactor flow rateFo action.
(a)
(b)Biodiesel concentration prediction result and relative error under reactor flow rateFc action.
(a)
(b)Reactor temperature prediction result and relative error under reactor flow rateFc action.
(a)
(b)Figures8 and 9 give the two output step response curves under two input Fo, Fc 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.Step response curve of biodiesel concentration and reactor temperature, respectively, underFo action.
(a)
(b)Step response curve of biodiesel concentration and reactor temperature, respectively, underFc action.
(a)
(b)


### 5.2. Dynamic Simulation

To validate the effectiveness and immunity in two-layer predictive control, the models obtained in Section5.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 Fo is constrained between 0 and 0.2 m3/s, the input Fc is constrained between 0 and 0.1 m3/s, and the output CE is constrained between 3.0536 kmol/m3 and 3.196 kmol/m3, the output T 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 domainN=200, prediction horizon P=200, control horizon M=20. The weight coefficient values of weight matrix Q and R 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 forCE with Fo control loop are kp=-6e-5, ki=-0.05, and kd=0; the parameters for T with Fc control loop are kp=-0.02, ki=-0.001, and kd=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.Biodiesel concentration and controller moves of two-layer predictive controller and PID controller.
(a)
(b)Reactor temperature and controller moves of two-layer predictive controller and PID controller.
(a)
(b)As Figures10 and 11 show, the two-layer predictive controller starts running at the time t=0. The results of steady state optimization are
(55)yss=[3.196,337.77],uss=[0.073,0.0062].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 timet=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 Figures10 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 (Tc0), feed temperature (T0), triglyceride initial concentration (CTG0), and stirrer rotational speed (N). After the system has attained the steady state, The nominal values of Tc0, T0 were increased 3 K, respectively, and CTG0, N were increased 5%, respectively, at the time t=1000s. Figures 12 and 13 show the biodiesel concentration and reactor temperature profiles when these disturbance variables were introduced.Biodiesel concentration and controller moves of four individual disturbance variables effects.
(a)
(b)Reactor temperature and controller moves of four individual disturbance variables effects.
(a)
(b)Figures12 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 CTG0 has the highest effect, with an overshoot of less than 0.01 kmol/m3. For the reactor temperature loop, the feed temperature T0 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.


## 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. Section4 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 betweena and -a according to the following rules:
(43)P[u(t)=-u(t-1)]=psw,P[u(t)=u(t-1)]=1-psw,
where psw is transition probability; Tmin is defined as the sampling time of the signal held constant; Tsw is time interval of twice conversion. The average conversion time and power spectrum are, respectively,
(44)ETsw=Tminpsw,Φu(ω)=(1-q2)Tmin1-2qcosTminω+q2,q=1-2psw.Next, least squares (LS) identification method is used to estimate the process model parameters. Suppose an MIMO plant with m input p output, for the ith output of the finite impulse response (FIR) model, is described as
(45)yi(k)=∑j=1m∑l=1Nhijluj(k-l).Consider experimental tests of collecting input sequence(46)u1(1)u1(2)⋯u1(L)⋮⋮um(1)um(2)⋯um(L)
and output sequence
(47)y1(1)y1(2)⋯y1(L)⋮⋮yp(1)yp(2)⋯yp(L).Consider matching between data and models; the introduction of residuals for each output can be independently expressed as follows:(48)yi(k)=φ(k)θi+e(k).Matrix form is written as(49)yi=Φθi+e,
where(50)yi=[yi(N+1)yi(N+2)⋮yi(L)],e=[e(N+1)e(N+2)⋮e(L)],Φ=[u1(N)u1(N-1)⋯u1(1)um(N)um(N-1)⋯um(1)u1(N+1)u1(N)⋯u1(2)um(N+1)um(N)⋯um(2)⋮⋮⋮⋯⋮⋮⋮u1(L-1)u1(L-2)u1(L-N)um(L-1)um(L-2)um(L-N)].Minimize the squared residuals(51)minJ=eTe=[y-Φθ]T[y-Φθ].Obtain the optimal estimate(52)θ^=[ΦTΦ]-1ΦTy.For the model predictive controller design, theFIR 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:
(53)gj=∑i=1jhj.Coefficients matrix ofFSR is
(54)Glu=[s11ls12l⋯s1mls21ls22l⋯s2ml⋮⋮⋱⋮sp1lsp2l⋯spml].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 areTsw=65, amp=0.1 the parameters of GBN applied to the second input are Tsw=65, amp=0.005, both the conversion probabilities are taken to be Psw=1/Tsw. Simulation time t=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 rateFo and coolant flow rate Fc, respectively, predicted value, actual value, and the relative error of two outputs biodiesel concentration CE and reactor temperature T were shown in Figures 4, 5, 6, and 7. Figures 4–7 show that relative error is small enough, and the model can describe CE and T change trends under Fo and Fc.Biodiesel concentration prediction result and relative error under reactor flow rateFo action.
(a)
(b)Reactor temperature prediction result and relative error under reactor flow rateFo action.
(a)
(b)Biodiesel concentration prediction result and relative error under reactor flow rateFc action.
(a)
(b)Reactor temperature prediction result and relative error under reactor flow rateFc action.
(a)
(b)Figures8 and 9 give the two output step response curves under two input Fo, Fc 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.Step response curve of biodiesel concentration and reactor temperature, respectively, underFo action.
(a)
(b)Step response curve of biodiesel concentration and reactor temperature, respectively, underFc action.
(a)
(b)


## 5.2. Dynamic Simulation

To validate the effectiveness and immunity in two-layer predictive control, the models obtained in Section5.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 Fo is constrained between 0 and 0.2 m3/s, the input Fc is constrained between 0 and 0.1 m3/s, and the output CE is constrained between 3.0536 kmol/m3 and 3.196 kmol/m3, the output T 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 domainN=200, prediction horizon P=200, control horizon M=20. The weight coefficient values of weight matrix Q and R 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 forCE with Fo control loop are kp=-6e-5, ki=-0.05, and kd=0; the parameters for T with Fc control loop are kp=-0.02, ki=-0.001, and kd=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.Biodiesel concentration and controller moves of two-layer predictive controller and PID controller.
(a)
(b)Reactor temperature and controller moves of two-layer predictive controller and PID controller.
(a)
(b)As Figures10 and 11 show, the two-layer predictive controller starts running at the time t=0. The results of steady state optimization are
(55)yss=[3.196,337.77],uss=[0.073,0.0062].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 timet=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 Figures10 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 (Tc0), feed temperature (T0), triglyceride initial concentration (CTG0), and stirrer rotational speed (N). After the system has attained the steady state, The nominal values of Tc0, T0 were increased 3 K, respectively, and CTG0, N were increased 5%, respectively, at the time t=1000s. Figures 12 and 13 show the biodiesel concentration and reactor temperature profiles when these disturbance variables were introduced.Biodiesel concentration and controller moves of four individual disturbance variables effects.
(a)
(b)Reactor temperature and controller moves of four individual disturbance variables effects.
(a)
(b)Figures12 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 CTG0 has the highest effect, with an overshoot of less than 0.01 kmol/m3. For the reactor temperature loop, the feed temperature T0 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.



---

*Source: 587841-2013-10-10.xml*