# Optimal Control of Complex Systems Based on Improved Dual Heuristic Dynamic Programming Algorithm

Authors: Hui Li; Yongsui Wen; Wenjie Sun
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

When applied to solving the data modeling and optimal control problems of complex systems, the dual heuristic dynamic programming (DHP) technique, which is based on the BP neural network algorithm (BP-DHP), has difficulty in prediction accuracy, slow convergence speed, poor stability, and so forth. In this paper, a dual DHP technique based on Extreme Learning Machine (ELM) algorithm (ELM-DHP) was proposed. Through constructing three kinds of network structures, the paper gives the detailed realization process of the DHP technique in the ELM. The controller designed upon the ELM-DHP algorithm controlled a molecular distillation system with complex features, such as multivariability, strong coupling, and nonlinearity. Finally, the effectiveness of the algorithm is verified by the simulation that compares DHP and HDP algorithms based on ELM and BP neural network. The algorithm can also be applied to solve the data modeling and optimal control problems of similar complex systems.


---

## Body

## 1. Introduction

With the increase of the dynamic complexity of controlled objects and their widespread application, some dynamic models of controlled objects can be obtained by mathematical methods. But some of the controlled objects are too complex to establish an accurate mathematical model, such as a robot system and a chemical generation process. Even though the mathematical model of complex systems has been established, it will be usually a high-order nonlinear time-varying complex differential equation. So, it could not describe the system accurately, and it is too difficult to analyze and process the data, and hence the unknown model must be learned by the observational data [1]. This is because of the great fault tolerance, self-adaptation, self-organization, and learning and memory ability of neural networks, which provides a new method for the modeling of complex systems [2–5]. However, BP networks, RBF networks, and SVM have some defects such as slow convergence, which causes a big gap between the approximation ability and the actual demand of complex systems [6, 7]. In order to solve this problem, Huang GB proposed single hidden layer feedforward network (SLFN) training—Extreme learning Machine (ELM)—in 2006 [8]. The ELM gives the weights and thresholds of the weights randomly and then calculates the output weights by the regularization principle, which can still approach any continuous system [6, 7]. It has been proved that the SLFN hidden layer node parameter randomly accessed does not affect the convergence ability, and it also makes the speed of learning the ELM thousands of times faster than the traditional BP network and SVM.In recent years, dynamic programming has been used to solve the optimal control problem, the “curse of dimensionality” problem, so no optimal solution could be obtained [9]. In 1977, Werbos proposed a heuristic dynamic programming (HDP) and dual heuristic programming (DHP) concept and proposed a method of approximate dynamic programming (adaptive/approximate dynamic programming, ADP) to solve the “curse of dimensionality” problem [10]. Werbos defines “intelligence” as the brain’s ability to learn a utility function maximally in a complex, unknown, nonlinear environment [11]. ADP is the general scheme for learning approximate optimal action strategies. Therefore, ADP can be regarded as a key method which is able to design the intelligent system of a brain. According to the basic principle, realization structure, and current development of the ADP method, Lewis and others gave a summary and prospect of the research and pointed out that ADP is an effective data-driven method [12–15]. ADP can realize the optimization control of the nonlinear system by using neural networks based on online data and control information to approximate the performance index function of the optimal control law, without a mathematical model of a nonlinear control system [16]. To develop neural dynamic programming results and relax the system dynamic requirements, Zhong et al. proposed a new goal representation ADP online optimization control structure for nonlinear systems. But the network could not output the derivative function information of the cost function directly based on the implementation structure of HDP, and the control effect of the HDP structure needed improvement [17]. In fact, some studies show that DHP and GDHP can be controlled better than HDP in the structure of ADP method to some extent [18, 19].In general, the study on the optimal control of nonlinear systems of ADP based on the traditional neural network has made great progress. But the ADP still has problems of slow response and poor stability. In this paper, ELM algorithm gives random input weights and thresholds to improve the response speed of the DHP algorithm and the stability of the DHP algorithm is improved by calculating the output weights by the regularization principle. In order to verify the validity of the algorithm, the ELM-DHP was designed to control the molecular distillation system with multivariability, nonlinearity, strong coupling, and large delay.


## 2. Algorithm Principle

### 2.1. DHP Algorithm Principle

The discrete-time nonlinear dynamic system is described as follows:(1)xk+1=Fxk,uk,k,k=0,1,….In formula (1), xk∈Rn represents the state vector of the system, uk∈Rn represents control variables, and F represents the system function.The performance index function (also called the cost function) corresponding to the system is(2)Jxk,k=∑i=k∞γi-kηxi,ui,i,where η is the utility function, γ is the discount factor (0<γ≤1), J is the cost function of state xk, and J depends on the initial time k and the initial state xk. For DHP, the purpose of dynamic programming is to select a control sequence ui,i=k,k+1,…,l, which minimizes the function ∂Jxk/∂xk.The DHP structure is shown in Figure1, which contains three neural networks: model network, critic network, and action network. The neural network has a powerful function of universal approximation, so the model network can be used to model the unknown nonlinear or complex nonlinear system and make the DHP method widely used. The input of the critic network is a state variable. The output of the critic network is approximation performance index functionJ on the statex derivative, which is also known as the costate. The action network, also known as “Actor,” represents the mapping between system state variables and control variables [20–22].Figure 1
The structural diagram of the DHP algorithm.∂
J
x
k
/
∂
x
k is based on the iteration of the derivative for performance index function and utility function to state.(3)∂Jxk∂xk=∂ηxk,uxk∂xk+γ∂Jxk+1∂xk.In (3), uxk is a feedback control variable, and costates ∂Jxk/∂xk and ∂Jxk+1/∂xk are the outputs of the critic network. If the weight of the critic network is set to ω, the right type of formula (1) is set to(4)exk,ω=∂ηxk,uxk∂xk+γ∂Jxk+1∂xk.At the same time, the left type of formula (1) can be written as ∂Jxk,ω/∂xk. By adjusting the weights ω of the critic network, the least-mean-square-error function is as follows:(5)ω∗=argminω⁡∂Jxk,ω∂xk-exk,ω2.According to the principle of optimality, the optimal control should satisfy the first-order differential necessary condition.(6)∂J∗xk∂uk=∂ηxk,uk∂uk+∂J∗xk+1∂uk=∂ηxk,uk∂uk+∂J∗xk+1∂xk+1∂Fxk,uk∂uk.So, the optimal quantity is obtained:(7)u∗=argminu⁡∂Jxk∂uk-∂ηxk,uk∂uk-∂J∗xk+1∂xk+1∂Fxk,uk∂uk.In formula (7), ∂J∗xk+1/∂xk+1 is the optimal costate, satisfying formula (5).From (1) to (7), we can conclude that the optimal control quantity of the DHP method can be obtained directly by the costate. Compared with the HDP method which obtained the optimal control by the relationship between the weights (ω) of the critic network and the input–output, the method of DHP has more computational efforts, but better control effect [9].


### 2.2. ELM Algorithm Principle

For a standard SLFN withL hidden layer neurons learning N arbitrary distinct samples xj,tj, xj=[xj1,xj2,…,xjn]T∈Rn, tj=[tj1,tj2,…,tjm]T∈Rm and activation function g· are mathematically modeled as [23](8)∑i=1Lβigwi·xj+bi=oj,j=1,2,…N,bi,βi∈R,ai∈R,where oj=[oj1,oj2,…,ojm]T∈Rm is the model output of the network, wi=wi1,wi2,…,winT is the input weight matrix between the input layer neuron and the ith hidden layer neuron, βi is the output weight matrix between the ith hidden layer neuron and the output layer neurons, bi is the threshold of the ith neuron in the hidden layer, and wi·xj is the inner product of wi and xj.The learning objective of the SLFN is to minimize the output error. Error can be expressed as(9)∑j=1Ntj-oj=0.The presence ofwi, xj, and bi makes(10)∑i=1Lβigwi·xj+bi=tj,j=1,2,…,N.So, (10) can be written as Hβ=T, where(11)Hw1,…,wL,b1,…,bN,x1,…,xN=gw1·x1+b1⋯gwL·x1+bL⋯⋯⋯gw1·xN+b1⋯gwL·xN+bLN×L,β=β1T⋮βLTL×m,T=t1T⋮tNTN×m.H is a hidden layer output matrix of ELM. So, the training of ELM is equivalent to the least-squares solution β^ of linear system Hβ=T.(12)Hw^i,b^iβ^-T=minw,b,β⁡Hw^i,b^i·β-T.In (12), i=1,2,…,L, (12) is equivalent to minimizing the loss function(13)E=∑j=1N∑i=1Lβi·gwi·xj+bi-tj2.Huang et al. [23] proved that the minimum value of the least-squares solution of the linear system satisfies the following.(
1) Minimum Training Error. The special solution β^=H-1T is one of the least-squares solutions of a general linear system Hβ=T,H-1 which is a generalized inverse matrix of H.(
2) Smallest Norm of Weights and Best Generalization Capability. Further, the special solution β^=H-1T has the smallest norm among all of the least-squares solutions of Hβ=T:β^=H-1T≤β, ∀β∈{β:Hβ-T≤Hz-T,∀z∈RN×N}. The generalization ability of SLFN with minimum weight is independent of the number of parameters [24]. The smaller the weight, the stronger the generalization ability of SLFN.(
3) Special Solution. The least-squares solution of Hβ=T is unique.


### 2.3. Proof the Stability of ELM-DHP

The stability of the ELM-DHP algorithm is proved (i.e., the output error of the system is 0). The discrete nonlinear system is controlled by the ELM-DHP algorithm, and the three networks of the ELM-DHP algorithm are all based on the fixed ELM implementation. Therefore, it just needs to be proved that ELM can approximate the discrete nonlinear system by 0 error.The ELM learning algorithm is chosen as a SLFN withL hidden layer neurons. The N arbitrary distinct samples xj,tj, where L≤N, xj=[xj1,xj2,…,xjn]T∈Rn, and tj=[tj1,tj2,…,tjm]T∈Rm of the nonlinear discrete system and nonlinear activation function g·, are mathematically modeled as formula (8).ELM learns a large number of samples generally, and the number of neurons in the hidden layer is far less than the number of samples,L≪N. So, we only need to prove that the learning error of ELM was 0 when L≤N. Huang et al. [7, 23, 25] proved in detail that the SLFN with L neurons can approximate any arbitrary sample xj,tj at any small error; that is,(14)e=∑j=1Ntj-oj=Hβ-T=0.The work above proves that the learning error of ELM is 0 (i.e., the stability of the ELM-DHP algorithm).


## 2.1. DHP Algorithm Principle

The discrete-time nonlinear dynamic system is described as follows:(1)xk+1=Fxk,uk,k,k=0,1,….In formula (1), xk∈Rn represents the state vector of the system, uk∈Rn represents control variables, and F represents the system function.The performance index function (also called the cost function) corresponding to the system is(2)Jxk,k=∑i=k∞γi-kηxi,ui,i,where η is the utility function, γ is the discount factor (0<γ≤1), J is the cost function of state xk, and J depends on the initial time k and the initial state xk. For DHP, the purpose of dynamic programming is to select a control sequence ui,i=k,k+1,…,l, which minimizes the function ∂Jxk/∂xk.The DHP structure is shown in Figure1, which contains three neural networks: model network, critic network, and action network. The neural network has a powerful function of universal approximation, so the model network can be used to model the unknown nonlinear or complex nonlinear system and make the DHP method widely used. The input of the critic network is a state variable. The output of the critic network is approximation performance index functionJ on the statex derivative, which is also known as the costate. The action network, also known as “Actor,” represents the mapping between system state variables and control variables [20–22].Figure 1
The structural diagram of the DHP algorithm.∂
J
x
k
/
∂
x
k is based on the iteration of the derivative for performance index function and utility function to state.(3)∂Jxk∂xk=∂ηxk,uxk∂xk+γ∂Jxk+1∂xk.In (3), uxk is a feedback control variable, and costates ∂Jxk/∂xk and ∂Jxk+1/∂xk are the outputs of the critic network. If the weight of the critic network is set to ω, the right type of formula (1) is set to(4)exk,ω=∂ηxk,uxk∂xk+γ∂Jxk+1∂xk.At the same time, the left type of formula (1) can be written as ∂Jxk,ω/∂xk. By adjusting the weights ω of the critic network, the least-mean-square-error function is as follows:(5)ω∗=argminω⁡∂Jxk,ω∂xk-exk,ω2.According to the principle of optimality, the optimal control should satisfy the first-order differential necessary condition.(6)∂J∗xk∂uk=∂ηxk,uk∂uk+∂J∗xk+1∂uk=∂ηxk,uk∂uk+∂J∗xk+1∂xk+1∂Fxk,uk∂uk.So, the optimal quantity is obtained:(7)u∗=argminu⁡∂Jxk∂uk-∂ηxk,uk∂uk-∂J∗xk+1∂xk+1∂Fxk,uk∂uk.In formula (7), ∂J∗xk+1/∂xk+1 is the optimal costate, satisfying formula (5).From (1) to (7), we can conclude that the optimal control quantity of the DHP method can be obtained directly by the costate. Compared with the HDP method which obtained the optimal control by the relationship between the weights (ω) of the critic network and the input–output, the method of DHP has more computational efforts, but better control effect [9].


## 2.2. ELM Algorithm Principle

For a standard SLFN withL hidden layer neurons learning N arbitrary distinct samples xj,tj, xj=[xj1,xj2,…,xjn]T∈Rn, tj=[tj1,tj2,…,tjm]T∈Rm and activation function g· are mathematically modeled as [23](8)∑i=1Lβigwi·xj+bi=oj,j=1,2,…N,bi,βi∈R,ai∈R,where oj=[oj1,oj2,…,ojm]T∈Rm is the model output of the network, wi=wi1,wi2,…,winT is the input weight matrix between the input layer neuron and the ith hidden layer neuron, βi is the output weight matrix between the ith hidden layer neuron and the output layer neurons, bi is the threshold of the ith neuron in the hidden layer, and wi·xj is the inner product of wi and xj.The learning objective of the SLFN is to minimize the output error. Error can be expressed as(9)∑j=1Ntj-oj=0.The presence ofwi, xj, and bi makes(10)∑i=1Lβigwi·xj+bi=tj,j=1,2,…,N.So, (10) can be written as Hβ=T, where(11)Hw1,…,wL,b1,…,bN,x1,…,xN=gw1·x1+b1⋯gwL·x1+bL⋯⋯⋯gw1·xN+b1⋯gwL·xN+bLN×L,β=β1T⋮βLTL×m,T=t1T⋮tNTN×m.H is a hidden layer output matrix of ELM. So, the training of ELM is equivalent to the least-squares solution β^ of linear system Hβ=T.(12)Hw^i,b^iβ^-T=minw,b,β⁡Hw^i,b^i·β-T.In (12), i=1,2,…,L, (12) is equivalent to minimizing the loss function(13)E=∑j=1N∑i=1Lβi·gwi·xj+bi-tj2.Huang et al. [23] proved that the minimum value of the least-squares solution of the linear system satisfies the following.(
1) Minimum Training Error. The special solution β^=H-1T is one of the least-squares solutions of a general linear system Hβ=T,H-1 which is a generalized inverse matrix of H.(
2) Smallest Norm of Weights and Best Generalization Capability. Further, the special solution β^=H-1T has the smallest norm among all of the least-squares solutions of Hβ=T:β^=H-1T≤β, ∀β∈{β:Hβ-T≤Hz-T,∀z∈RN×N}. The generalization ability of SLFN with minimum weight is independent of the number of parameters [24]. The smaller the weight, the stronger the generalization ability of SLFN.(
3) Special Solution. The least-squares solution of Hβ=T is unique.


## 2.3. Proof the Stability of ELM-DHP

The stability of the ELM-DHP algorithm is proved (i.e., the output error of the system is 0). The discrete nonlinear system is controlled by the ELM-DHP algorithm, and the three networks of the ELM-DHP algorithm are all based on the fixed ELM implementation. Therefore, it just needs to be proved that ELM can approximate the discrete nonlinear system by 0 error.The ELM learning algorithm is chosen as a SLFN withL hidden layer neurons. The N arbitrary distinct samples xj,tj, where L≤N, xj=[xj1,xj2,…,xjn]T∈Rn, and tj=[tj1,tj2,…,tjm]T∈Rm of the nonlinear discrete system and nonlinear activation function g·, are mathematically modeled as formula (8).ELM learns a large number of samples generally, and the number of neurons in the hidden layer is far less than the number of samples,L≪N. So, we only need to prove that the learning error of ELM was 0 when L≤N. Huang et al. [7, 23, 25] proved in detail that the SLFN with L neurons can approximate any arbitrary sample xj,tj at any small error; that is,(14)e=∑j=1Ntj-oj=Hβ-T=0.The work above proves that the learning error of ELM is 0 (i.e., the stability of the ELM-DHP algorithm).


## 3. Implementation of the ELM-DHP Algorithm

The ELM-DHP algorithm includes three networks: model network, critic network, and action network. The hidden layer of the three networks is a sigmoidal bipolar function and the output layer is a purelin linear function. The realization process of the ELM-DHP algorithm is studied by using the discrete-time nonlinear dynamic programming ofn-dimensional state vector and m-dimensional control vector as the research object.


### 3.1. Network Model

The model network adoptsm+n-lm-n structure. The m+n inputs are the n components of the state vector x(k) in the k moments and the m components of the predicted output u(k) of the action network to state x(k) in the system of k-1 moments. The n output is the n components of the prediction vector x^k+1 to the state vector xk+1 in the system of k+1 moments. The model network has lm hidden layer neurons. The structure of the model network is shown in Figure 2.Figure 2
The structure of the model network.The model network is trained offline, and the calculation process is as follows.The input layer to the hidden layer weight matrixWm1 and the hidden layer threshold matrix B=b1,b2…,blm are randomly generated. Define the input vector Mk and the expected output vector x^(k) of the model network in k moments:(15)Mk=u1k,u2k,…,umk,x1k,x2k,…,xnk,x^k=x^1k,x^2k,…,x^nk.Calculate the output matrixmh2k of the hidden layer in the model network(16)mh1jk=∑i=1n+mMik·Wm1ijk+bjk,j=1,2,…,lm,mh2jk=1-e-mh1jk1+e-mh1jk,j=1,2,…,lm,where mh1jk is the input of the jth node in the model network hidden layer, mh2jk is the output of the jth node in the model network hidden layer, and mh2=[mh21,mh22,…,mh2lm]∈Rm.Calculate the weightsWm2k from the hidden layer to the output layer:(17)x^k+1=mh2k×Wm2k.According to the idea of the ELM, the error is minimized as(18)Emk+1=∑i=1nxik+1-x^ik+12=0.In equality (18), xik+1 is the expected output ith output layer neurons of the model network.W
m
2
k is equivalent to solving the least-squares solution W^m2k of the linear system mh2k×Wm2k=x^k+1:(19)mh2k×W^m2-x^k+1=minWm2mh2k×Wm2-x^k+1.The special solutionW^m2k of the weight matrix of the hidden layer and output layer in the model network is as follows:(20)W^m2k=mh2-1k×x^k+1,where mh2-1k is a generalized inverse matrix of mh2k in k moments.


### 3.2. Critic Network

The critic network is composed ofn-lc-n. The n inputs are the n components of the state vector x^k, and the output is the estimation of the state λk=∂Jk/∂xk, Jk=γJk+1+ηk. lc is the number of hidden layer neurons in the critic network. In the critic network, the weight matrix from the input layer to the hidden layer, the weight matrix from the hidden layer to the output layer, and the hidden layer threshold matrix of k time are, respectively, defined as Wc1, Wc2, Ck=c1k,…,clck. Figure 3 shows the structure of the critic network.Figure 3
The structure of the critic network.The critic network uses the least-squares method of ELM, whose forward calculation process is(21)ch1jk=∑i=1nx^ikWc1ijk+cjk,j=1,2,…,lc,ch2jk=1-e-ch1jk1+e-ch1jk,j=1,2,…,lc,λk=ch2k×Wc2k,where ch1j is the input of the jth node in the critic network hidden layer, ch2j is the output of the jth node in the critic network hidden layer, ch2k=ch21,ch22,…,ch2lc, and λk is the output of the critic network output layer. The inputs x^ik of the critic network come from the output of the model network and the outputs of the critic network are costate function ∂Jk/∂xk in the DHP. λak is expressed to the expected output of the critic network, which can be written as(22)λak=∂Jk∂xk.The training error of DHP critic network is minimized based on the idea of ELM.(23)Ec=∑kEck=12∑k∑j=1nEj2k=0,Eck=12∑j=1nEj2k=12Ek×ETk=0,where Eck is the error of the critic network in k moments and Ec is the error of all the time points in the critic network.According to the DHP structure and the definition of the expected outputsλak of the critic network, we can obtain(24)Ek=λk-∂+ηk∂xk-γ∂+Jk+1∂xk.In formula (24), ∂+ηk/∂xk and ∂+Jk+1/∂xk represent the notion that ηk and Jk+1 take the derivative of composite function xk.Based on (23) and (24), we can acquire (25)λk=∂+ηk∂xk+γ∂+Jk+1∂xk.According toCh2(k)×Wc2=λ(k), the weight Wc2 from the hidden layer to the output layer is equal to the least-squares solution W^c2 of the linear system Ch2(k)×Wc2=λ(k), and hence we get W^c2:(26)Ch2k×W^c2-λk=minWc2Ch2k×Wc2-λk,W^c2=Ch2-1k×λk.In formula (26), Ch2-1(k) is a generalized inverse matrix of Ch2(k). Based on the DHP structure and the chain rule [9], we can obtain(27)∂+Jk+1∂xk=12·λk+1×Wm2T×Wm1xk⊗1-mh2k⊗mh2k;1-mh2k⊗mhh2k︸totalnT+14λk+1×Wm2Tk⊗1-mh2k⊗mh2k;1-mh2k⊗mh2k︸totaln×Wm1uTk×Wa2Tk×Wa1k⊗1-ah2k⊗ah2k;1-ah2k⊗ah2k︸totalnT.In (27), Wm1u=Wm11:m,: represents front m lines of the weight matrix Wm1, and Wm1x=Wm1m+1:m+n,: represents from m+1th to m+nth line of the weight matrix Wm1.(28)∂+ηk∂xk=2xk+uk×Wa2Tk×Wa1k⊗1-ah2k⊗ah2k;1-ah2k⊗ah2k︸totalnT.


### 3.3. Action Network

The action network uses the structure ofn-ld-m. n inputs are the n components of the state vector xk of the system at k moments. m outputs are the m components of the control vector uk corresponding to the input state vector xk. ld represents the number of neurons in the action network hidden layer. Wa1 and Wa2 are, respectively, the weight matrix from the input layer to the hidden layer and the weight matrix from the hidden layer to the output layer in the action network. dk=d1k,…,dldk is the hidden layer threshold matrix of the action network. Figure 4 is the structure of the action network.Figure 4
The structure of the action network.The calculation process of the action network is as follows:(29)ah1jk=∑i=1nxik·Wa1ijk+djk,j=1,2,…,ld,ah2jk=1-e-ah1jk1+e-ah1jk,j=1,2,…,ld,uk=ah2k·Wa2k,where ah1jk is the input of the jth node and ah2jk is the output of the jth node in the action network hidden layer and ah2k=ah21,ah22,…,ah2ld. According to the idea of weight adjustment of ELM, the weight matrix Wa2 from the hidden layer to the output layer is obtained:(30)Wa2=ah2-1k×uk.In (30), ah2-1k is a generalized inverse matrix of ah2k and uk is the expected output of the action network. The weights of the network will be corrected if uk can be got. The inverse sigmoidal function is defined as g(·). The calculation process of uk is as follows: (31)A=∂+Uk∂xk+γ∂+Jk+1∂xk×Wc2-1,(32)B=gA×Wc1-1×Wm2-1,(33)ukxk=gB×Wm1-1.In (33), uk is the first m rows of matrix g(B)×Wm1-1.Defineux=g(B)∗Wm1-1; we have(34)uk=ux1:m,:.So,Wa2 can be got:(35)Wa2=ah2-1k×uk.


### 3.4. Training Strategy

In this paper, the model network of the DHP algorithm is trained by an offline method at first to obtain the weight matrix of the model network. Then, the action network and the critic network are trained simultaneously. Training strategies are as follows:(1) First, the model is trained by an offline method and the weight matrix of the model network is obtained.(2) Takingxk into the action network, uk can be obtained.(3) Takinguk and xk into the model network, x^k+1 will be obtained.(4) Takingx^k+1 into the critic network, λ^xk+1 can be obtained.(5) Calculate the expected output value of the critic network∂+ηk/∂xk.(6) Next, calculate the value of∂+ηk/∂xk.(7) Next, calculate and update the weights of the critic network.(8) Last, makek=k+1 and go back to the second step until λk=γ∗λ^k+1+∂+ηk/∂xk.


## 3.1. Network Model

The model network adoptsm+n-lm-n structure. The m+n inputs are the n components of the state vector x(k) in the k moments and the m components of the predicted output u(k) of the action network to state x(k) in the system of k-1 moments. The n output is the n components of the prediction vector x^k+1 to the state vector xk+1 in the system of k+1 moments. The model network has lm hidden layer neurons. The structure of the model network is shown in Figure 2.Figure 2
The structure of the model network.The model network is trained offline, and the calculation process is as follows.The input layer to the hidden layer weight matrixWm1 and the hidden layer threshold matrix B=b1,b2…,blm are randomly generated. Define the input vector Mk and the expected output vector x^(k) of the model network in k moments:(15)Mk=u1k,u2k,…,umk,x1k,x2k,…,xnk,x^k=x^1k,x^2k,…,x^nk.Calculate the output matrixmh2k of the hidden layer in the model network(16)mh1jk=∑i=1n+mMik·Wm1ijk+bjk,j=1,2,…,lm,mh2jk=1-e-mh1jk1+e-mh1jk,j=1,2,…,lm,where mh1jk is the input of the jth node in the model network hidden layer, mh2jk is the output of the jth node in the model network hidden layer, and mh2=[mh21,mh22,…,mh2lm]∈Rm.Calculate the weightsWm2k from the hidden layer to the output layer:(17)x^k+1=mh2k×Wm2k.According to the idea of the ELM, the error is minimized as(18)Emk+1=∑i=1nxik+1-x^ik+12=0.In equality (18), xik+1 is the expected output ith output layer neurons of the model network.W
m
2
k is equivalent to solving the least-squares solution W^m2k of the linear system mh2k×Wm2k=x^k+1:(19)mh2k×W^m2-x^k+1=minWm2mh2k×Wm2-x^k+1.The special solutionW^m2k of the weight matrix of the hidden layer and output layer in the model network is as follows:(20)W^m2k=mh2-1k×x^k+1,where mh2-1k is a generalized inverse matrix of mh2k in k moments.


## 3.2. Critic Network

The critic network is composed ofn-lc-n. The n inputs are the n components of the state vector x^k, and the output is the estimation of the state λk=∂Jk/∂xk, Jk=γJk+1+ηk. lc is the number of hidden layer neurons in the critic network. In the critic network, the weight matrix from the input layer to the hidden layer, the weight matrix from the hidden layer to the output layer, and the hidden layer threshold matrix of k time are, respectively, defined as Wc1, Wc2, Ck=c1k,…,clck. Figure 3 shows the structure of the critic network.Figure 3
The structure of the critic network.The critic network uses the least-squares method of ELM, whose forward calculation process is(21)ch1jk=∑i=1nx^ikWc1ijk+cjk,j=1,2,…,lc,ch2jk=1-e-ch1jk1+e-ch1jk,j=1,2,…,lc,λk=ch2k×Wc2k,where ch1j is the input of the jth node in the critic network hidden layer, ch2j is the output of the jth node in the critic network hidden layer, ch2k=ch21,ch22,…,ch2lc, and λk is the output of the critic network output layer. The inputs x^ik of the critic network come from the output of the model network and the outputs of the critic network are costate function ∂Jk/∂xk in the DHP. λak is expressed to the expected output of the critic network, which can be written as(22)λak=∂Jk∂xk.The training error of DHP critic network is minimized based on the idea of ELM.(23)Ec=∑kEck=12∑k∑j=1nEj2k=0,Eck=12∑j=1nEj2k=12Ek×ETk=0,where Eck is the error of the critic network in k moments and Ec is the error of all the time points in the critic network.According to the DHP structure and the definition of the expected outputsλak of the critic network, we can obtain(24)Ek=λk-∂+ηk∂xk-γ∂+Jk+1∂xk.In formula (24), ∂+ηk/∂xk and ∂+Jk+1/∂xk represent the notion that ηk and Jk+1 take the derivative of composite function xk.Based on (23) and (24), we can acquire (25)λk=∂+ηk∂xk+γ∂+Jk+1∂xk.According toCh2(k)×Wc2=λ(k), the weight Wc2 from the hidden layer to the output layer is equal to the least-squares solution W^c2 of the linear system Ch2(k)×Wc2=λ(k), and hence we get W^c2:(26)Ch2k×W^c2-λk=minWc2Ch2k×Wc2-λk,W^c2=Ch2-1k×λk.In formula (26), Ch2-1(k) is a generalized inverse matrix of Ch2(k). Based on the DHP structure and the chain rule [9], we can obtain(27)∂+Jk+1∂xk=12·λk+1×Wm2T×Wm1xk⊗1-mh2k⊗mh2k;1-mh2k⊗mhh2k︸totalnT+14λk+1×Wm2Tk⊗1-mh2k⊗mh2k;1-mh2k⊗mh2k︸totaln×Wm1uTk×Wa2Tk×Wa1k⊗1-ah2k⊗ah2k;1-ah2k⊗ah2k︸totalnT.In (27), Wm1u=Wm11:m,: represents front m lines of the weight matrix Wm1, and Wm1x=Wm1m+1:m+n,: represents from m+1th to m+nth line of the weight matrix Wm1.(28)∂+ηk∂xk=2xk+uk×Wa2Tk×Wa1k⊗1-ah2k⊗ah2k;1-ah2k⊗ah2k︸totalnT.


## 3.3. Action Network

The action network uses the structure ofn-ld-m. n inputs are the n components of the state vector xk of the system at k moments. m outputs are the m components of the control vector uk corresponding to the input state vector xk. ld represents the number of neurons in the action network hidden layer. Wa1 and Wa2 are, respectively, the weight matrix from the input layer to the hidden layer and the weight matrix from the hidden layer to the output layer in the action network. dk=d1k,…,dldk is the hidden layer threshold matrix of the action network. Figure 4 is the structure of the action network.Figure 4
The structure of the action network.The calculation process of the action network is as follows:(29)ah1jk=∑i=1nxik·Wa1ijk+djk,j=1,2,…,ld,ah2jk=1-e-ah1jk1+e-ah1jk,j=1,2,…,ld,uk=ah2k·Wa2k,where ah1jk is the input of the jth node and ah2jk is the output of the jth node in the action network hidden layer and ah2k=ah21,ah22,…,ah2ld. According to the idea of weight adjustment of ELM, the weight matrix Wa2 from the hidden layer to the output layer is obtained:(30)Wa2=ah2-1k×uk.In (30), ah2-1k is a generalized inverse matrix of ah2k and uk is the expected output of the action network. The weights of the network will be corrected if uk can be got. The inverse sigmoidal function is defined as g(·). The calculation process of uk is as follows: (31)A=∂+Uk∂xk+γ∂+Jk+1∂xk×Wc2-1,(32)B=gA×Wc1-1×Wm2-1,(33)ukxk=gB×Wm1-1.In (33), uk is the first m rows of matrix g(B)×Wm1-1.Defineux=g(B)∗Wm1-1; we have(34)uk=ux1:m,:.So,Wa2 can be got:(35)Wa2=ah2-1k×uk.


## 3.4. Training Strategy

In this paper, the model network of the DHP algorithm is trained by an offline method at first to obtain the weight matrix of the model network. Then, the action network and the critic network are trained simultaneously. Training strategies are as follows:(1) First, the model is trained by an offline method and the weight matrix of the model network is obtained.(2) Takingxk into the action network, uk can be obtained.(3) Takinguk and xk into the model network, x^k+1 will be obtained.(4) Takingx^k+1 into the critic network, λ^xk+1 can be obtained.(5) Calculate the expected output value of the critic network∂+ηk/∂xk.(6) Next, calculate the value of∂+ηk/∂xk.(7) Next, calculate and update the weights of the critic network.(8) Last, makek=k+1 and go back to the second step until λk=γ∗λ^k+1+∂+ηk/∂xk.


## 4. Simulation Analysis

### 4.1. Simulation Example Analysis

The molecular distillation technology was also called short films. When enough energy is obtained, the average free path that escapes from the surface of a liquid of light molecules differs from that of heavy molecules, which achieve the nonequilibrium liquid–liquid separation process under high vacuum conditions [26]. The molecular distillation technology has advantages of low temperature distillation, short heating time, and high separation efficiency, and it is conducive to separate the material, that is, high boiling point, heat sensitivity, and high viscosity material separation. This technology is widely used in food, medicine, oil processing, and petrochemical industry [27–29]. Molecular distillation equipment can be divided into four types: stationary, falling film, scraped film, and centrifugal type [30]. At present, wiped film molecular distillation is the most widely used technology in scientific research and industrial production. The evaporation effect of the molecular distillation system is not only related to the size and shape of the evaporator and space, the distance to the surface evaporation condensation, the manufacturing process, and other types of equipment, but also connected with the pressure within the parameters of the feed flow rate, temperature of the evaporator, scraping, and other devices running the motor speed film process parameters [31]. In order to enhance the purification effect of molecular distillation, Wang et al. found that the head wave has an effect on the separation efficiency of molecular distillation by the study of the head wave [32]. Micov et al. studied the separation factors of the wiped film molecular distillation process and established a one-dimensional mathematical model [30]. Cvengros and Tkac established a mathematical equation which can be used to calculate the one-dimensional analysis mathematical equation of micro unit movement velocity in distillation equipment through the DSMC method and summarized the effects of evaporation temperature, distance, and vacuum degree and other related factors on the separation results [33]. Wu studied the simulation of the temperature, pressure, and reflux ratio on yield and purity by using the central response surface method combined with thin film evaporation and rectification coupling technology [34]. Although much research has been made, there are still many problems in molecular distillation system with multivariability, nonlinearity, strong coupling, and large delay. Therefore, the effectiveness of the ELM-DHP algorithm was verified by controlling the scraping film molecular distillation system.The current state variables of the molecular distillation system are determined by the amount of state variables in the preceding section of the system and the control variables in the previous stage. So, distillation temperature, evaporation pressure, wiper motor speed, feeding speed, andSchisandra yield and purity of the front section were used as the input of the ELM-DHP controller, and the currentSchisandra yield and purity were used as the output of the ELM-DHP controller.


### 4.2. Simulation Comparison

The structures of the model network, critic network, and action network were set as 6-20-2, 2-14-2, and 2-5-4 through experiment, respectively. In the process of system identification, the weight values of the three networks between the input layer and the hidden layer are selected in the range[-0.1,01]. 600 groups of data are collected to study, and 150 groups of data were used as the test set. Firstly, we need to train the model network offline; the least-squares solutions were calculated as the weight matrix between the hidden layer and the output layer. Then, we complete the training of the model network and keep its weight unchanged. The 50 time steps of the model network are shown in Figures 5 and 6.Figure 5
BP network model fitting ability test.
(a)
(b)Figure 6
ELM network model fitting ability test.
(a)
(b)Figures5 and 6 show that the predicted values of the BP network and the ELM algorithm are in good agreement with the expected values. Figures 7 and 8 show that the maximum error of the BP network in the prediction of the state is 0.4, but the maximum error of ELM for the state prediction is about 0.06. Thus, it can be concluded that ELM has higher prediction accuracy and better generalization ability.Figure 7
BP network model prediction error.
(a)
(b)Figure 8
ELM network model prediction error.
(a)
(b)Parameter setting will affect the convergence speed of the algorithm to a certain extent. After the experiment, the discount factor was chosen asγ=0.9. Next, the weights of the critic network and the action network from the hidden layer to the output layer are calculated. Then, the training of the critic network and action network is set to 150 steps with 100 training epochs for each step.In addition, in order to compare with the HDP and DHP technology based on BP neural network, controllers designed by BP-HDP, BP-DHP, and ELM-DHP were proposed. Four controllers are used to control the wiped film molecular distillation system, respectively, and the 50 time steps of the simulation results are shown in Figures9–14. Figures 9–12 show that the control quantities of BP-HDP, BP-DHP, ELM-HDP, and ELM-DHP controllers achieve stable control in 45 steps, 35 steps, 18 steps, and 7 steps individually. Thus, it can be concluded that the HDP and DHP algorithms based on ELM can achieve faster response speed. There will be a larger fluctuation when the controlled variables of the HDP controller achieve stability. So, it can be concluded that the DHP algorithm has a higher stability. The results of Figures 13 and 14 are shown in Table 1. The purification effect increases with yield and purity and the best purification effect is 100%, but it is impossible to achieve. It can be seen in Table 1 that the optimal state quantities derived by ELM-HDP and ELM-DHP were 5% higher than BP-HDP and BP-DHP, and the optimal state of ELM-DHP is slightly higher than that of ELM-HDP. In the above analysis, the superiority and effectiveness of the ELM algorithm can be demonstrated clearly.Table 1
Optimal state of controller.
BP-HDP
BP-DHP
ELM-HDP
ELM-DHP
Yield (%)
91.95
92.39
96.58
97.35
Purity (%)
89.25
91.11
97.12
97.18Figure 9
The molecular distillation temperature of optimum control quantity.Figure 10
The variable feed rate of optimum control quantity.Figure 11
The wiped film motor speed of optimum control quantity.Figure 12
The molecular distillation pressure of optimum control quantity.Figure 13
The yield of optimum state.Figure 14
The purity of optimum state.


## 4.1. Simulation Example Analysis

The molecular distillation technology was also called short films. When enough energy is obtained, the average free path that escapes from the surface of a liquid of light molecules differs from that of heavy molecules, which achieve the nonequilibrium liquid–liquid separation process under high vacuum conditions [26]. The molecular distillation technology has advantages of low temperature distillation, short heating time, and high separation efficiency, and it is conducive to separate the material, that is, high boiling point, heat sensitivity, and high viscosity material separation. This technology is widely used in food, medicine, oil processing, and petrochemical industry [27–29]. Molecular distillation equipment can be divided into four types: stationary, falling film, scraped film, and centrifugal type [30]. At present, wiped film molecular distillation is the most widely used technology in scientific research and industrial production. The evaporation effect of the molecular distillation system is not only related to the size and shape of the evaporator and space, the distance to the surface evaporation condensation, the manufacturing process, and other types of equipment, but also connected with the pressure within the parameters of the feed flow rate, temperature of the evaporator, scraping, and other devices running the motor speed film process parameters [31]. In order to enhance the purification effect of molecular distillation, Wang et al. found that the head wave has an effect on the separation efficiency of molecular distillation by the study of the head wave [32]. Micov et al. studied the separation factors of the wiped film molecular distillation process and established a one-dimensional mathematical model [30]. Cvengros and Tkac established a mathematical equation which can be used to calculate the one-dimensional analysis mathematical equation of micro unit movement velocity in distillation equipment through the DSMC method and summarized the effects of evaporation temperature, distance, and vacuum degree and other related factors on the separation results [33]. Wu studied the simulation of the temperature, pressure, and reflux ratio on yield and purity by using the central response surface method combined with thin film evaporation and rectification coupling technology [34]. Although much research has been made, there are still many problems in molecular distillation system with multivariability, nonlinearity, strong coupling, and large delay. Therefore, the effectiveness of the ELM-DHP algorithm was verified by controlling the scraping film molecular distillation system.The current state variables of the molecular distillation system are determined by the amount of state variables in the preceding section of the system and the control variables in the previous stage. So, distillation temperature, evaporation pressure, wiper motor speed, feeding speed, andSchisandra yield and purity of the front section were used as the input of the ELM-DHP controller, and the currentSchisandra yield and purity were used as the output of the ELM-DHP controller.


## 4.2. Simulation Comparison

The structures of the model network, critic network, and action network were set as 6-20-2, 2-14-2, and 2-5-4 through experiment, respectively. In the process of system identification, the weight values of the three networks between the input layer and the hidden layer are selected in the range[-0.1,01]. 600 groups of data are collected to study, and 150 groups of data were used as the test set. Firstly, we need to train the model network offline; the least-squares solutions were calculated as the weight matrix between the hidden layer and the output layer. Then, we complete the training of the model network and keep its weight unchanged. The 50 time steps of the model network are shown in Figures 5 and 6.Figure 5
BP network model fitting ability test.
(a)
(b)Figure 6
ELM network model fitting ability test.
(a)
(b)Figures5 and 6 show that the predicted values of the BP network and the ELM algorithm are in good agreement with the expected values. Figures 7 and 8 show that the maximum error of the BP network in the prediction of the state is 0.4, but the maximum error of ELM for the state prediction is about 0.06. Thus, it can be concluded that ELM has higher prediction accuracy and better generalization ability.Figure 7
BP network model prediction error.
(a)
(b)Figure 8
ELM network model prediction error.
(a)
(b)Parameter setting will affect the convergence speed of the algorithm to a certain extent. After the experiment, the discount factor was chosen asγ=0.9. Next, the weights of the critic network and the action network from the hidden layer to the output layer are calculated. Then, the training of the critic network and action network is set to 150 steps with 100 training epochs for each step.In addition, in order to compare with the HDP and DHP technology based on BP neural network, controllers designed by BP-HDP, BP-DHP, and ELM-DHP were proposed. Four controllers are used to control the wiped film molecular distillation system, respectively, and the 50 time steps of the simulation results are shown in Figures9–14. Figures 9–12 show that the control quantities of BP-HDP, BP-DHP, ELM-HDP, and ELM-DHP controllers achieve stable control in 45 steps, 35 steps, 18 steps, and 7 steps individually. Thus, it can be concluded that the HDP and DHP algorithms based on ELM can achieve faster response speed. There will be a larger fluctuation when the controlled variables of the HDP controller achieve stability. So, it can be concluded that the DHP algorithm has a higher stability. The results of Figures 13 and 14 are shown in Table 1. The purification effect increases with yield and purity and the best purification effect is 100%, but it is impossible to achieve. It can be seen in Table 1 that the optimal state quantities derived by ELM-HDP and ELM-DHP were 5% higher than BP-HDP and BP-DHP, and the optimal state of ELM-DHP is slightly higher than that of ELM-HDP. In the above analysis, the superiority and effectiveness of the ELM algorithm can be demonstrated clearly.Table 1
Optimal state of controller.
BP-HDP
BP-DHP
ELM-HDP
ELM-DHP
Yield (%)
91.95
92.39
96.58
97.35
Purity (%)
89.25
91.11
97.12
97.18Figure 9
The molecular distillation temperature of optimum control quantity.Figure 10
The variable feed rate of optimum control quantity.Figure 11
The wiped film motor speed of optimum control quantity.Figure 12
The molecular distillation pressure of optimum control quantity.Figure 13
The yield of optimum state.Figure 14
The purity of optimum state.


## 5. Summary

For those problems which the BP-DHP algorithm has, such as poor prediction accuracy, slow convergence speed, and poor stability, the ELM-DHP algorithm was studied in this paper to solve the data modeling and optimal control problem of the wiped film molecular distillation system with complex features such as multivariability, strong coupling, nonlinearity, and large time delay as an example. The ELM-DHP controller was designed to control the molecular distillation system and a simulation verification was carried out. When compared with the ELM-HDP, BP-HDP, and BP-DHP algorithms, the prediction accuracy of ELM is higher than that of the BP neural network, and the response speed and stability of the ELM- HDP and ELM-DHP algorithms are higher than those achieved by the BP network, which shows the superiority of ELM. Compared with other algorithms, the response speed of ELM-DHP is more than two times that of the other algorithms, and the optimal state achieved by ELM-DHP is closer to the ideal result. Thus, the ELM-DHP algorithm is better than BP-HDP, BP-DHP, and ELM-HDP algorithms. The ELM-DHP algorithm does not depend on the specific mechanism model and is only in accordance with the relevant experimental data, so the algorithm can also solve the optimal control problem of similar complex systems which have features such as multiple variables, strong coupling, nonlinearity, and large time delay.



---

*Source: 5476415-2017-10-09.xml*