# Optimal Tracking Control for a Discrete Time Nonlinear Nuclear Power System

Authors: Zhenhua Luan; Mengxuan Wang; Yuzhen Zhang; Qinglai Wei; Tianmin Zhou; Zhiwu Guo; Jun Ling
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

Recently, increasing attention has been paid to nuclear power control with the appeals of clean energy and demands of power regulation to integrate into the power grid. However, a nuclear power system is a discrete-time (DT) nonlinear and complicated system, where the parameters entangle with intrinsic states. Furthermore, the need for huge computational ability due to the high-level order property in the nuclear reactor model causes many difficulties in the power control of nuclear industries. In this study, a new scheme of optimal tracking control for DT nonlinear nuclear power systems is provided to accomplish the power control of a 2500-MW pressurized water reactor (PWR) nuclear power plant. The proposed approach based on the value iteration method is a novel algorithm in the human intelligence community, which has a basic actor-critic structure with neural networks (NNs). The new approach has some modifications, where the cost function is redefined by leveraging the higher-order polynomial to substitute neural networks in the entire actor critic architecture. Simulation results of the 2500-MW PWR nuclear power plant are given to demonstrate the effectiveness of the developed method.


---

## Body

## 1. Introduction

Considering the issues of environmental deterioration, e.g., air pollution due to excessive fossil fuel consumption, it is significant that humans develop clean energy technology to ease this situation. Nuclear energy is almost the most rapidly developing clean power to provide power to the power grid. However, currently adopted control strategies have problems such as the intrinsic nonlinearity of nuclear reactor systems and varying parameters following the power level. In fact, over decades of development for nuclear power industry technology and control policy, many outstanding researchers have made excellent progress in this field.Since the last century, one mature control strategy called PID control policy has deeply affected the power control in the nuclear industry [1, 2]. With the advancement in control technology, some model predictive control (MPC) and multimodel adaptive control theories focused on local linearization to approximate the nonlinearity of nuclear power systems and have also been applied in this area [3–5]. The control algorithm extensively applied in nuclear power control is fuzzy control or its combinations with other control theories to address different demands. Wu leveraged the parallel distribution compensation (PDC) method-based T-S fuzzy control to restrict the nonlinearity of a nuclear power system [6], and Eliasi designed an appropriate controller for the UTSG water level in nuclear power plants using a fuzzy control policy and MPC algorithm [7]. Many researchers have used different methods to tackle various problems. For example, Gang employed a radial basis function neural network (RBFN) to guarantee the correctness in identifying the nuclear steam generator process dynamics [8]. Wang applied the adaptive control method and guaranteed a cost control method in nuclear power control problems [9].The aforementioned approaches mostly relate to the linearization procedure, which largely omits the numerical error toward a nonlinear model. The intrinsic nonlinearity of the high-order nuclear model is ignored. To better satisfy the demands of tracking control problems, it is necessary to propose a new control strategy in this area. With the development in the intelligent control community, researchers are pursuing reinforcement learning (RL) algorithms to solve nonlinear problems in practice. Adaptive dynamic programming (ADP), which was proposed by Werbos, plays an important role in reinforcement learning-based control policy [10–13], and it is well known as a self-learning optimal control policy. The well-studied iteration method is the policy iteration (PI) algorithm and value iteration (VI) algorithm. Among the iteration methods, the value iteration algorithm is one type of the most crucial iterative ADP algorithms. It has been studied in many types of research [14–16]. To find the optimal control policy of discrete-time affine nonlinear systems, Al-Tamimi and Lewis used heuristic dynamic programming (HDP) to fulfill the design purpose [17]. Wei. Q. [18] proposed a new value iteration method, which mainly focuses on optimal control for DT nonlinear systems. This study also provided detailed proof of the iterative control policy and illustrated that the value function was monotonically nondecreasing, which implies that it will converge to the extremum.To satisfy the demands of industrial systems, optimal tracking control ADP methods have been deeply investigated. There are also ADP techniques [19–21] to obtain solutions of optimal tracking problems with various system dynamics, such as partially unknown system models or completely unknown system models. Related optimal tracking control techniques have been applied in many industrial plants in recent years [22–26].In this study, a value-iteration optimal tracking control method is developed for DT nonlinear systems. The main contributions of this study are summarized as follows:(1)
Compared with the traditional control methods dealing with DT nonlinear models of the 2500 MW pressurized water reactor (PWR) systems [1, 2, 4, 5], a self-learning optimal tracking controller is designed to satisfy complex nonlinear behaviors of the 2500 MW PWR nuclear system(2)
The developed value-iteration method guarantees the control law converges to a near-optimal control solution and the admissibility of iterative control laws is analyzedIn this study, our major work is to design an optimal tracking controller for a 2500 MW PWR nuclear power plant by combining the properties of the value iteration and actor critic algorithm. The 2500-MW PWR nuclear power plant is introduced in Section2, and the discrete definition is given. In Section 3, the details of this proposed algorithm are thoroughly described. The implementation of the proposed method and simulation works are provided in Section 4. Finally, the conclusions are drawn in Section 5.


## 2. Nonlinear 2500-MW PWR Nuclear Power Plant

The famous nuclear system model is based on Mann’s model without xenon poisoning, which consists of a core full lump and two coolant lumps. The discrete version and its transformation are also given in this section.


### 2.1. Nonlinear 2500-MW PWR Nuclear Power Plant

This fifth-order nonlinear PWR model includes the point kinetics equations, six delayed neutron groups, two equations for the lumped coolant outlet temperature and average fuel temperature, and the reactive equation of the control rod [27, 28]. Multiple sets of delayed neutron point reactor dynamic equations are described as follows:(1)dntdt=ρ−βΛnt+∑i=16λicit,dcitdt=βiΛnt−λicit.To reduce the computational work caused by six delayed neutron point reactor dynamic equations, the simple method is to use single delayed neutron point kinetics equations to approximate multiple sets of delayed neutron point reactor dynamic equations [29]. Thus, the entire PWR model is summarized as follows:(2)dnrdt=ρ−βΛnr+βΛcr,dcrdt=λnr−λcr,dTfdt=ffP0aμfnr−OμfTf+O2μfTl+O2μfTe,dTldt=1−ffP0aμcnr+OμcTf−2M+O2μcTl+2M−O2μcT,edρrdt=GrZr,ρ=ρr+αfTf−Tf0+αcTl−Tl02+αcTe−Te02,where nr is the neutron density relative to the initial equilibrium density, %, cr is the delay neutron density relative to its initial equilibrium density, %, Tf is the average fuel temperature, °C, Tl is the coolant temperature at the core outlet, °C, ρr is reactivity contributed by the control rod movement, and Zr is the speed of the control rod. The remaining specifications are illustrated in Nomenclature [30]. Te always approximates Te0 in the PWR model. The state nr can be described as a percentage factor of the full power level, since the reactor power is expressed as Pt=P0anr.In addition, five parameters vary withnr, which causes severe instability of the nuclear reactor power model and increases the control complexity. The remaining parameters and specific relation are shown in Tables 1 and 2. With the lifting (lowering) load of the PWR model, the varying parameters will lead to a sharp difference in the solutions of the dynamic model. Thus, the model will be uncontrollable, and the solutions of this dynamic model may become divergent. Linearizing nuclear systems to realize various control goals has been a common method in traditional control policies in the past few years. However, there should be a new approach in nonlinear systems to solve these problems.Table 1
Nomenclature.
NomenclatureGrReactivity worth of the control rod bankP0aNominal core power (MW)ΛMean neutron lifetime (s)βThe total fraction of effective theith group of delayed neutronsOThe heat transfer coefficient between fuel and coolant (MW·s/°C)MHeat capacity of mass flux of coolant (MW·s/ °C)αcThe reactivity coefficient of coolant temperatureαfThe reactivity coefficient of fuel temperature.μcHeat capacity of coolant (MW·s/°C)μfHeat capacity of fuel (MW·s/°C)λEquivalent single-group delayed neutron precursor nuclear decay constantρCore reactivityTeAverage inlet temperature of coolant (°C)Te0Initial equilibrium inlet temperature of coolant (°C)Tl0Initial equilibrium outlet temperature of coolant (°C).Tf0Initial equilibrium fuel temperature (°C).Table 2
Parameters of the PWR nuclear reactor.
ParametersValueParametersValueβ0.0065μf26.3λ0.15s−1ff0.98Λ0.00002(s)Te0290°CGr0.0145Tf0673.8°CP02500(MW)Tl0302.2°Cμc160nr/9+54.022O5nr/3+4.933αfnr−4.24×10−5M28.0nr+74.0αc−4.0nr−17.3×10−5


### 2.2. System Discretization and Transformation

The optimal tracking control problem can be considered minimizing the real dynamic trajectory with the desired trajectory. Depending on the model, we letx=x1,x2,x3,x4,x5T; thus, the control-oriented nuclear power system can be defined as(3)x˙=fx+gxu,where x1=nr,x2=cr,x3=Tf,x4=Tl,andx5=ρr; thus, fx and gx can be derived as(4)fx=x5+αfx3−Tf0−βΛ+αcTl−Tl02Λx1+βΛx2λx1−λx2ffP0aμfx1−Oμfx3+O2μfx4+O2μfTe1−ffP0aμcx1+Oμcx3−2M+O2μcx4+2M−O2μcTe0(5)gx=0000Gr.As the above descriptions show, there is only one controllable signalZr, i.e., the speed of the control rod. We let control variable u=Zr.According to the definition,xk.=xk+1−xk/ΔT, and we have the discretization version of the PWR power model:(6)xk+1=fxk+gxkuk.Basically, we define the optimal tracking control problem as obtaining an optimal control strategy so that the system can track the reference state, and the desired trajectoryxd can be expressed as(7)xd,k+1=fxd,k+gxd,kud,k.The tracking error of the state is defined asεk=xk−xd,k. Thus, the dynamic error system can be described as follows:(8)εk+1=fxk+gxkuk−xd,k.Additionally, we must confirm an initial steady-state control policyud, and the error of the controller can be defined as uε,k=uk−ud,k, so the dynamic tracking error system can also be written as(9)εk+1=fεεk+gεεkuε,k,where fεεk=fεk+xd,k−xd,k and gεεk=gεk+xd,k.We define the utility function as follows:(10)Uεk,uε,k=εkTQεk+uε,kTRuε,k.Thus, the tracking error cost function is written as(11)ℑεk,uε,k=∑k∞Uεk,uε,k.From the principle of optimality, the DT Hamiltonian function is derived as(12)Hεk,uε,k,∇ℑε,k=Uεk,uε,k+ℑε,k+1−ℑε,k.Thus, the HJB equation can be written as(13)minuεHεk,uε,k,∇ℑε,k∗=0.Then, the optimal tracking control law for the error system is derived:(14)uε,k∗=12R−1gεTεk∂ℑεεk+1∂εk+1.Finally, we obtain the standard optimal control law as(15)uk∗=ud,k+uε,k∗..For linear systems, the HJB equation is reduced to the Riccati equation. However, due to the nonlinearity of the nuclear power system, it is extremely intractable to solve the HJB equation (13) for the nonlinear system. Thus, the value iteration (VI) method based on the actor critic NN algorithm will be adopted to find an approximate optimal solution of the HJB equation (13), which implies that nothing is required about the knowledge of the model drifts or the command generator.


## 2.1. Nonlinear 2500-MW PWR Nuclear Power Plant

This fifth-order nonlinear PWR model includes the point kinetics equations, six delayed neutron groups, two equations for the lumped coolant outlet temperature and average fuel temperature, and the reactive equation of the control rod [27, 28]. Multiple sets of delayed neutron point reactor dynamic equations are described as follows:(1)dntdt=ρ−βΛnt+∑i=16λicit,dcitdt=βiΛnt−λicit.To reduce the computational work caused by six delayed neutron point reactor dynamic equations, the simple method is to use single delayed neutron point kinetics equations to approximate multiple sets of delayed neutron point reactor dynamic equations [29]. Thus, the entire PWR model is summarized as follows:(2)dnrdt=ρ−βΛnr+βΛcr,dcrdt=λnr−λcr,dTfdt=ffP0aμfnr−OμfTf+O2μfTl+O2μfTe,dTldt=1−ffP0aμcnr+OμcTf−2M+O2μcTl+2M−O2μcT,edρrdt=GrZr,ρ=ρr+αfTf−Tf0+αcTl−Tl02+αcTe−Te02,where nr is the neutron density relative to the initial equilibrium density, %, cr is the delay neutron density relative to its initial equilibrium density, %, Tf is the average fuel temperature, °C, Tl is the coolant temperature at the core outlet, °C, ρr is reactivity contributed by the control rod movement, and Zr is the speed of the control rod. The remaining specifications are illustrated in Nomenclature [30]. Te always approximates Te0 in the PWR model. The state nr can be described as a percentage factor of the full power level, since the reactor power is expressed as Pt=P0anr.In addition, five parameters vary withnr, which causes severe instability of the nuclear reactor power model and increases the control complexity. The remaining parameters and specific relation are shown in Tables 1 and 2. With the lifting (lowering) load of the PWR model, the varying parameters will lead to a sharp difference in the solutions of the dynamic model. Thus, the model will be uncontrollable, and the solutions of this dynamic model may become divergent. Linearizing nuclear systems to realize various control goals has been a common method in traditional control policies in the past few years. However, there should be a new approach in nonlinear systems to solve these problems.Table 1
Nomenclature.
NomenclatureGrReactivity worth of the control rod bankP0aNominal core power (MW)ΛMean neutron lifetime (s)βThe total fraction of effective theith group of delayed neutronsOThe heat transfer coefficient between fuel and coolant (MW·s/°C)MHeat capacity of mass flux of coolant (MW·s/ °C)αcThe reactivity coefficient of coolant temperatureαfThe reactivity coefficient of fuel temperature.μcHeat capacity of coolant (MW·s/°C)μfHeat capacity of fuel (MW·s/°C)λEquivalent single-group delayed neutron precursor nuclear decay constantρCore reactivityTeAverage inlet temperature of coolant (°C)Te0Initial equilibrium inlet temperature of coolant (°C)Tl0Initial equilibrium outlet temperature of coolant (°C).Tf0Initial equilibrium fuel temperature (°C).Table 2
Parameters of the PWR nuclear reactor.
ParametersValueParametersValueβ0.0065μf26.3λ0.15s−1ff0.98Λ0.00002(s)Te0290°CGr0.0145Tf0673.8°CP02500(MW)Tl0302.2°Cμc160nr/9+54.022O5nr/3+4.933αfnr−4.24×10−5M28.0nr+74.0αc−4.0nr−17.3×10−5


## 2.2. System Discretization and Transformation

The optimal tracking control problem can be considered minimizing the real dynamic trajectory with the desired trajectory. Depending on the model, we letx=x1,x2,x3,x4,x5T; thus, the control-oriented nuclear power system can be defined as(3)x˙=fx+gxu,where x1=nr,x2=cr,x3=Tf,x4=Tl,andx5=ρr; thus, fx and gx can be derived as(4)fx=x5+αfx3−Tf0−βΛ+αcTl−Tl02Λx1+βΛx2λx1−λx2ffP0aμfx1−Oμfx3+O2μfx4+O2μfTe1−ffP0aμcx1+Oμcx3−2M+O2μcx4+2M−O2μcTe0(5)gx=0000Gr.As the above descriptions show, there is only one controllable signalZr, i.e., the speed of the control rod. We let control variable u=Zr.According to the definition,xk.=xk+1−xk/ΔT, and we have the discretization version of the PWR power model:(6)xk+1=fxk+gxkuk.Basically, we define the optimal tracking control problem as obtaining an optimal control strategy so that the system can track the reference state, and the desired trajectoryxd can be expressed as(7)xd,k+1=fxd,k+gxd,kud,k.The tracking error of the state is defined asεk=xk−xd,k. Thus, the dynamic error system can be described as follows:(8)εk+1=fxk+gxkuk−xd,k.Additionally, we must confirm an initial steady-state control policyud, and the error of the controller can be defined as uε,k=uk−ud,k, so the dynamic tracking error system can also be written as(9)εk+1=fεεk+gεεkuε,k,where fεεk=fεk+xd,k−xd,k and gεεk=gεk+xd,k.We define the utility function as follows:(10)Uεk,uε,k=εkTQεk+uε,kTRuε,k.Thus, the tracking error cost function is written as(11)ℑεk,uε,k=∑k∞Uεk,uε,k.From the principle of optimality, the DT Hamiltonian function is derived as(12)Hεk,uε,k,∇ℑε,k=Uεk,uε,k+ℑε,k+1−ℑε,k.Thus, the HJB equation can be written as(13)minuεHεk,uε,k,∇ℑε,k∗=0.Then, the optimal tracking control law for the error system is derived:(14)uε,k∗=12R−1gεTεk∂ℑεεk+1∂εk+1.Finally, we obtain the standard optimal control law as(15)uk∗=ud,k+uε,k∗..For linear systems, the HJB equation is reduced to the Riccati equation. However, due to the nonlinearity of the nuclear power system, it is extremely intractable to solve the HJB equation (13) for the nonlinear system. Thus, the value iteration (VI) method based on the actor critic NN algorithm will be adopted to find an approximate optimal solution of the HJB equation (13), which implies that nothing is required about the knowledge of the model drifts or the command generator.


## 3. Algorithm Analysis

The detailed convergence properties of this proposed value iteration algorithm are illustrated in this section, and the details of actor critic NN will also be discussed.


### 3.1. Analysis of the Value Iteration Algorithm for Tracking Control of Nonlinear Systems

Considering the nonaffine nonlinear system, for an infinite time optimal tracking problem, the goal is to obtain an optimal controller such that the statexk tracks the specified reference trajectory xd,k.Remark 1.
For many nonlinear systems, there is a feedback controlud,k that makes (9) work. For example, with regard to DT nonlinear systems (7) with invertible gxd,k, the desired control ud,k can be derived as(16)ud,k=g−1xd,kxd,k+1−fxd,k.
From equation (11) of Section 2.2, the quadratic cost function of tracking errors εk is defined as(17)ℑε0,uε,0¯=∑k=0∞εkQεk+uε,kTRuε,k,where uε,0¯=uε,0,uε,1,… and Q· and R· are positive definite functions.
To obtain an optimum tracking control law that tracks the reference statexd,k and minimizes the tracking error cost function (18), we can redefine the optimal tracking error cost function as follows:(18)ℑ∗εk=infuε,k¯ℑεk,uε,k¯.
Based on Bellman’s principle of optimality,ℑ∗εk satisfies Bellman equation:(19)ℑ∗εk=minεkUεk,uε,k+ℑ∗εk+1.
Then, the optimal tracking control law is obtained by(20)u∗εk=argminuε,kUεk,uε,k+ℑ∗εk.
Given the above formulation, we can derive the tracking error performance index function as follows:(21)ℑ∗εk=Uεk,uε,k∗+ℑ∗εk+1.
Letφεk be a positive definite function for εk∈ℝ5, and the initial tracking value function is(22)ℑ0εk=φεk.
The optimal control lawv0εk can be obtained by(23)v0εk=argminuεk,kUεk,uε,k+ℑ0εk+1,where ℑ0εk+1=φεk. For i=1,2,3…, in this iterative value function algorithm, the value function is updated through(24)ℑiεk=minuε,kUεk,uε,k+ℑi−1εk+1,and the control policy is improved by(25)πiεk=argminuε,kUεk,uε,k+ℑi−1εk+1.Theorem 1.
For the tracking error cost functionℑiεk and control law πiεk obtained by (22)–(25), we have α, β, γ, and η satisfying 0<η≤γ<∞ and 0≤α≤β<1, respectively. If ∀εk, we have(26)ηUεk,πk≤ℑ∗εk+1≤γUεk,πk,(27)αℑ∗εk≤ℑ0εk≤βℑ∗εk,are satisfied uniformly; thus, the iterative value function ℑiεk satisfies(28)1+α−11+γ−1iℑ∗εk≤ℑiεk≤1+β−11+η−1iℑ∗εk.Theorem 2.
For cost functionℑiεk and control law πiεk obtained by (22)–(25), we have α, β, γ, and η satisfying 1≤α≤β<∞. If ∀εk, indexes (26) and (27) hold uniformly; then, the value function Jiεk satisfies equation (28)
The proof is thoroughly provided in [17].Corollary 1.
Fori=0,1,…, πiεk and Viεk are obtained by (22)–(25). Let α, β, γ, and η be constants that satisfy 0<η≤γ<∞ and 0≤α≤β<∞, respectively. If ∀εk, inequalities (26) and (27) hold uniformly. Then, the iterative value function ℑiεk converges to the optimal cost function ℑ∗εk, i.e.,(29)limi⟶∞ℑiεk=ℑ∗εk.
Based on the aforementioned analysis, we can conclude that the iterative tracking error performance index function will converge to the optimum asi⟶∞, which is independent of the initial value function φεk.
According to Lyapunov stability principle,ℑiεk is a Lyapunov function. Since the utility function Uεk,πk is a positive definite function and ℑi0=0, it should be noticed that ℑiεk is a positive definite function as well. We let(30)ℑiεk+1−ℑiεk≤ℑiεk+1−ℑi+1εk,=−Uεk,πk≤0,where the error tracking control law πk is admissible.


### 3.2. Actor Critic NN Implementation of the Value Iteration Algorithm for DT Nonlinear Systems

The actor critic NN has been employed in various fields to approximate the cost function and optimal controller. For example, the optimal tracking applied on partially unknown DT nonlinear systems [31] and rigorous proof for this method are provided. The actor critic structure also has good performance in the tracking control problem for continuous-time nonlinear systems [32, 33].With regard to optimal tracking control algorithms of DT nonlinear systems, it is quite difficult to directly obtain solutions by solving (13). Thus, the actor critic network structure with the flowchart of the nuclear power system is given in Figure 1, which describes the inner procedures of the method.Figure 1
Actor-critic structure of the optimal tracking error system.While the tracking error system is fed with a specific initial state and desired trajectory, the error will be calculated by a utility function. Simultaneously, the critic network will be trained in a way that minimizes the utility function. Under the training process, the actor network will behave like an optimal controller. To avoid overfitting, a specified threshold is given at the beginning of the training procedure so that it can be stopped in time.Inspired by Abu-Khalaf [34], an optimal control algorithm with a high-order polynomial was proposed to substitute the neural unit. This technique is introduced in this actor critic NN structure to obtain a better approximate effect.To solve (24) and (25), we let the tracking error cost function V^i be approximated by a critic NN:(31)V^iε=∑j=1LwcjΘjε=WViΘε,where we have the approximate activation function Θε≜Θ1εk,Θ2εk,…,ΘLεkT, the weight vector WVi=wc1,wc2,…,wcL, and L is the number of neural units in the hidden structure of the critic NN.Then, the iterative formulation can be obtained as follows:(32)WVi+1Θεi+1=Uεk,uε,k+WViΘiε.For each sampleεk related to xk, we formulate the definition as follows:(33)WVi+1=ZTZ−1ZTζi,where Z=Θεi and ζi=Uεk,uε,k+WViΘiε. Then, the weights of critic network are obtained. We let uiε be approximated by an actor NN:(34)u^iε=∑j=1Mwajδjε=Wuiδε,where we have the approximate activation function δε≜δ1εk,δ2εk,…,δMεkT and weight vector Wui=wa1,wa2,…,waM. M is the number of neural units in actor NN.According to (24) and (30), we will tune the weights of the critic NN at each iteration of this VI algorithm. Our goal is to minimize the residual error between each V^iε to obtain a new target function as follows:(35)dεk,εk+1,Wui,WVi=εkTQεk+u^iTεkRu^iεk+V^iεk+1,=εkTQεk+u^iTεkRu^iεk,+V^iεk+1.Similarly, the actor NN is applied to evaluate and approximate the optimal tracking control policy. We will tune the weights of action NN to solve (20) at each iteration of this VI algorithm. According to u^iεk,Wui, from (33), we can rewrite (14) as(36)Wui=argminu^iεkTQεk+u^iTεk,δRu^iεk,δ+V^iεk+1,where εk+1=fεεk+gεεku^iεk,δ and δ are updated by the same method as the weights of the critic NN. For getting the approximation of the weights of actor critic NNs, the least square (LS) method is utilized to solve the weights of NNs.


## 3.1. Analysis of the Value Iteration Algorithm for Tracking Control of Nonlinear Systems

Considering the nonaffine nonlinear system, for an infinite time optimal tracking problem, the goal is to obtain an optimal controller such that the statexk tracks the specified reference trajectory xd,k.Remark 1.
For many nonlinear systems, there is a feedback controlud,k that makes (9) work. For example, with regard to DT nonlinear systems (7) with invertible gxd,k, the desired control ud,k can be derived as(16)ud,k=g−1xd,kxd,k+1−fxd,k.
From equation (11) of Section 2.2, the quadratic cost function of tracking errors εk is defined as(17)ℑε0,uε,0¯=∑k=0∞εkQεk+uε,kTRuε,k,where uε,0¯=uε,0,uε,1,… and Q· and R· are positive definite functions.
To obtain an optimum tracking control law that tracks the reference statexd,k and minimizes the tracking error cost function (18), we can redefine the optimal tracking error cost function as follows:(18)ℑ∗εk=infuε,k¯ℑεk,uε,k¯.
Based on Bellman’s principle of optimality,ℑ∗εk satisfies Bellman equation:(19)ℑ∗εk=minεkUεk,uε,k+ℑ∗εk+1.
Then, the optimal tracking control law is obtained by(20)u∗εk=argminuε,kUεk,uε,k+ℑ∗εk.
Given the above formulation, we can derive the tracking error performance index function as follows:(21)ℑ∗εk=Uεk,uε,k∗+ℑ∗εk+1.
Letφεk be a positive definite function for εk∈ℝ5, and the initial tracking value function is(22)ℑ0εk=φεk.
The optimal control lawv0εk can be obtained by(23)v0εk=argminuεk,kUεk,uε,k+ℑ0εk+1,where ℑ0εk+1=φεk. For i=1,2,3…, in this iterative value function algorithm, the value function is updated through(24)ℑiεk=minuε,kUεk,uε,k+ℑi−1εk+1,and the control policy is improved by(25)πiεk=argminuε,kUεk,uε,k+ℑi−1εk+1.Theorem 1.
For the tracking error cost functionℑiεk and control law πiεk obtained by (22)–(25), we have α, β, γ, and η satisfying 0<η≤γ<∞ and 0≤α≤β<1, respectively. If ∀εk, we have(26)ηUεk,πk≤ℑ∗εk+1≤γUεk,πk,(27)αℑ∗εk≤ℑ0εk≤βℑ∗εk,are satisfied uniformly; thus, the iterative value function ℑiεk satisfies(28)1+α−11+γ−1iℑ∗εk≤ℑiεk≤1+β−11+η−1iℑ∗εk.Theorem 2.
For cost functionℑiεk and control law πiεk obtained by (22)–(25), we have α, β, γ, and η satisfying 1≤α≤β<∞. If ∀εk, indexes (26) and (27) hold uniformly; then, the value function Jiεk satisfies equation (28)
The proof is thoroughly provided in [17].Corollary 1.
Fori=0,1,…, πiεk and Viεk are obtained by (22)–(25). Let α, β, γ, and η be constants that satisfy 0<η≤γ<∞ and 0≤α≤β<∞, respectively. If ∀εk, inequalities (26) and (27) hold uniformly. Then, the iterative value function ℑiεk converges to the optimal cost function ℑ∗εk, i.e.,(29)limi⟶∞ℑiεk=ℑ∗εk.
Based on the aforementioned analysis, we can conclude that the iterative tracking error performance index function will converge to the optimum asi⟶∞, which is independent of the initial value function φεk.
According to Lyapunov stability principle,ℑiεk is a Lyapunov function. Since the utility function Uεk,πk is a positive definite function and ℑi0=0, it should be noticed that ℑiεk is a positive definite function as well. We let(30)ℑiεk+1−ℑiεk≤ℑiεk+1−ℑi+1εk,=−Uεk,πk≤0,where the error tracking control law πk is admissible.


## 3.2. Actor Critic NN Implementation of the Value Iteration Algorithm for DT Nonlinear Systems

The actor critic NN has been employed in various fields to approximate the cost function and optimal controller. For example, the optimal tracking applied on partially unknown DT nonlinear systems [31] and rigorous proof for this method are provided. The actor critic structure also has good performance in the tracking control problem for continuous-time nonlinear systems [32, 33].With regard to optimal tracking control algorithms of DT nonlinear systems, it is quite difficult to directly obtain solutions by solving (13). Thus, the actor critic network structure with the flowchart of the nuclear power system is given in Figure 1, which describes the inner procedures of the method.Figure 1
Actor-critic structure of the optimal tracking error system.While the tracking error system is fed with a specific initial state and desired trajectory, the error will be calculated by a utility function. Simultaneously, the critic network will be trained in a way that minimizes the utility function. Under the training process, the actor network will behave like an optimal controller. To avoid overfitting, a specified threshold is given at the beginning of the training procedure so that it can be stopped in time.Inspired by Abu-Khalaf [34], an optimal control algorithm with a high-order polynomial was proposed to substitute the neural unit. This technique is introduced in this actor critic NN structure to obtain a better approximate effect.To solve (24) and (25), we let the tracking error cost function V^i be approximated by a critic NN:(31)V^iε=∑j=1LwcjΘjε=WViΘε,where we have the approximate activation function Θε≜Θ1εk,Θ2εk,…,ΘLεkT, the weight vector WVi=wc1,wc2,…,wcL, and L is the number of neural units in the hidden structure of the critic NN.Then, the iterative formulation can be obtained as follows:(32)WVi+1Θεi+1=Uεk,uε,k+WViΘiε.For each sampleεk related to xk, we formulate the definition as follows:(33)WVi+1=ZTZ−1ZTζi,where Z=Θεi and ζi=Uεk,uε,k+WViΘiε. Then, the weights of critic network are obtained. We let uiε be approximated by an actor NN:(34)u^iε=∑j=1Mwajδjε=Wuiδε,where we have the approximate activation function δε≜δ1εk,δ2εk,…,δMεkT and weight vector Wui=wa1,wa2,…,waM. M is the number of neural units in actor NN.According to (24) and (30), we will tune the weights of the critic NN at each iteration of this VI algorithm. Our goal is to minimize the residual error between each V^iε to obtain a new target function as follows:(35)dεk,εk+1,Wui,WVi=εkTQεk+u^iTεkRu^iεk+V^iεk+1,=εkTQεk+u^iTεkRu^iεk,+V^iεk+1.Similarly, the actor NN is applied to evaluate and approximate the optimal tracking control policy. We will tune the weights of action NN to solve (20) at each iteration of this VI algorithm. According to u^iεk,Wui, from (33), we can rewrite (14) as(36)Wui=argminu^iεkTQεk+u^iTεk,δRu^iεk,δ+V^iεk+1,where εk+1=fεεk+gεεku^iεk,δ and δ are updated by the same method as the weights of the critic NN. For getting the approximation of the weights of actor critic NNs, the least square (LS) method is utilized to solve the weights of NNs.


## 4. Numerical Simulations

In this section, numerical results are given to demonstrate the validity of this value iteration optimal tracking method. Experimental simulations of the performance index and weights of actor critic NNs are provided. This developed method is an offline policy with an initial random control policy.


### 4.1. Actor Critic NN’s Implementation of the Value Iteration Algorithm

It is generally known that NNs can be leveraged to approximate any functions on prescribed compact sets. We choose an error compact set to train this actor critic NN to obtain an offline tracking policy. As predefined in Section3, the critic NN is approximated as ℑε=WcLϕe with 15 neurons L=1,2,…,15, and the weights are WcL=Wc1,Wc2,…,Wc15. The actor NN is chosen with 5 neurons, and the weights are WaL=Wa1,Wa2,Wa3,Wa4,Wa5.At the beginning of this algorithm execution, generally, the matrices areQ=10I5×5 and R=0.1I1×1, where I is the identity matrix. The tracking error compact sets are randomly set as the difference between the initial state and the desired trajectory. We chose the control period as 1s periteration, and the total iteration epoch was 50. Simultaneously, we chose the %60 initial working point withx60=0.6,0.6,530.8572,306.5198,−0.0047T as state inputs and the %80 power working point with x80=0.8,0.8,594.0060,310.7468,−0.0018T as the desired tracking trajectory. Under the proposed value iteration algorithm, the tracking error performance index value, as shown in Figure 2, is monotonically nondecreasing and converges to ℑε∗, which is identical to our analysis results above.Figure 2
Performance index function of the tracking error.Additionally, when the tracking error performance index function converges, the training process of the actor-critic NN will stop. As shown in Figure3 and 4, the weights of the actor-critic NN converge to a steady solution, which implies that a perfect approximator of the optimal tracking controller is obtained.Figure 3
Weights of the critic network.Figure 4
Weights of the actor network.


### 4.2. Application on the PWR Nuclear Power System

In this section, we apply the calculated control law on this DT nonlinear PWR power model, and the implementation results are shown in Figure5–10.Figure 5
Power-level load-varying curve of this PWR power plant.Figure 6
Delay neutron density relative curve of this PWR power plant.Figure 7
Average temperature of the reactor core of this PWR power plant.Figure 8
Average temperature of the coolant outlet of this PWR power plant.Generally, the optimal tracking control of PWR power plants focuses on power-level adjustment, but there is high interference among different states in this nuclear power model. Thus, it is necessary to track all states in the PWR power model, and the tracking objects should guarantee the stability of each state and safety of nuclear plants.Figure 9
Reactor coefficient of the control rod of this PWR power plant.Figure 10
Five tracking errors of this PWR power system.As shown in Figures5–9, we give a %20 step increase signal to this PWR power system. These 5 figures demonstrate that the 5 states catch the desired trajectory in less than 50s. We use a 5th-order DT nonlinear PWR power system in this study, which implies that Xenon poison is without consideration. Although the PWR power model is quite simplified, all states of this model are difficult to track.As shown in Figure5, the power level tracks the desired states without overshoot and oscillation. The average temperature of the reactor core and coolant outlet approximate the desired temperature; compared with the reference temperature, the maximum deviation is 0.031°C and 0.002°C, respectively. In addition, the reactor coefficient of the control rod tracks the desired curve in less than 40s. Based the aforementioned results, we also can see that our tracking method applied this nuclear system has no steady-state error and less regulation time. As shown in Figure 10, the tracking errors progressively converge to 0 with the proposal of this value iteration optimal tracking method.


## 4.1. Actor Critic NN’s Implementation of the Value Iteration Algorithm

It is generally known that NNs can be leveraged to approximate any functions on prescribed compact sets. We choose an error compact set to train this actor critic NN to obtain an offline tracking policy. As predefined in Section3, the critic NN is approximated as ℑε=WcLϕe with 15 neurons L=1,2,…,15, and the weights are WcL=Wc1,Wc2,…,Wc15. The actor NN is chosen with 5 neurons, and the weights are WaL=Wa1,Wa2,Wa3,Wa4,Wa5.At the beginning of this algorithm execution, generally, the matrices areQ=10I5×5 and R=0.1I1×1, where I is the identity matrix. The tracking error compact sets are randomly set as the difference between the initial state and the desired trajectory. We chose the control period as 1s periteration, and the total iteration epoch was 50. Simultaneously, we chose the %60 initial working point withx60=0.6,0.6,530.8572,306.5198,−0.0047T as state inputs and the %80 power working point with x80=0.8,0.8,594.0060,310.7468,−0.0018T as the desired tracking trajectory. Under the proposed value iteration algorithm, the tracking error performance index value, as shown in Figure 2, is monotonically nondecreasing and converges to ℑε∗, which is identical to our analysis results above.Figure 2
Performance index function of the tracking error.Additionally, when the tracking error performance index function converges, the training process of the actor-critic NN will stop. As shown in Figure3 and 4, the weights of the actor-critic NN converge to a steady solution, which implies that a perfect approximator of the optimal tracking controller is obtained.Figure 3
Weights of the critic network.Figure 4
Weights of the actor network.


## 4.2. Application on the PWR Nuclear Power System

In this section, we apply the calculated control law on this DT nonlinear PWR power model, and the implementation results are shown in Figure5–10.Figure 5
Power-level load-varying curve of this PWR power plant.Figure 6
Delay neutron density relative curve of this PWR power plant.Figure 7
Average temperature of the reactor core of this PWR power plant.Figure 8
Average temperature of the coolant outlet of this PWR power plant.Generally, the optimal tracking control of PWR power plants focuses on power-level adjustment, but there is high interference among different states in this nuclear power model. Thus, it is necessary to track all states in the PWR power model, and the tracking objects should guarantee the stability of each state and safety of nuclear plants.Figure 9
Reactor coefficient of the control rod of this PWR power plant.Figure 10
Five tracking errors of this PWR power system.As shown in Figures5–9, we give a %20 step increase signal to this PWR power system. These 5 figures demonstrate that the 5 states catch the desired trajectory in less than 50s. We use a 5th-order DT nonlinear PWR power system in this study, which implies that Xenon poison is without consideration. Although the PWR power model is quite simplified, all states of this model are difficult to track.As shown in Figure5, the power level tracks the desired states without overshoot and oscillation. The average temperature of the reactor core and coolant outlet approximate the desired temperature; compared with the reference temperature, the maximum deviation is 0.031°C and 0.002°C, respectively. In addition, the reactor coefficient of the control rod tracks the desired curve in less than 40s. Based the aforementioned results, we also can see that our tracking method applied this nuclear system has no steady-state error and less regulation time. As shown in Figure 10, the tracking errors progressively converge to 0 with the proposal of this value iteration optimal tracking method.


## 5. Conclusion

It is well known that optimal tracking power-level control for DT nonlinear nuclear power systems is crucial for both regular operation and safety problems, and manual control is inefficient. However, the intrinsic nonlinearity and parameters that vary with the states cause difficulties in power-level control, and there are tracking issues.In this study, a value iteration-based actor critic NN algorithm is designed to obtain an optimal tracking control policy for the DT nonlinear nuclear power plant. The proposed algorithm performs well in tracking states, as shown in the simulation results, and can also swiftly calculate the optimal control law. Thus, we formulate the tracking control problem as HJB equations to solve it.



---

*Source: 7953358-2022-06-22.xml*