For the ADP method, the main thing is the value function approximation. In general, there are three methods to describe the value function approximation (Salas and Powell, 2013; Li and Jayaweera, 2015): lookup table, parametric approximation and nonparametric approximation. For example, in (Das and Ni, 2018), authors research about the battery storage systems operation in islanded microgrid considering battery lifetime characteristics, and the approximate value function is formulated based on lookup table idea. In (Li and Jayaweera, 2015), the authors use Q-learning method to define the approximate value function. In (Keerthisinghe et al., 2018), the piecewise linear function is used to build the approximate value function. In (Zeng et al., 2018), deep recurrent neural network learning is adopted to describe the approximate value function. The reference papers showed that ADP has better performance and lower computational burden.

Using the ADP method to optimal control the operation of microgrids has also attracted lots of attention.

The above section reviews the related work about scheduling algorithms for microgrid. In addition, when microgrid interconnects with the electricity/heat/gas utility grids, the operation of the electricity/heat/gas utility grids should also be considered.

For the coupled multi-energy networks operation, centralized optimization algorithm is often used to solve the optimal power flow. For example, authors in (Qin et al., 2020) study the operation of integrated energy systems consisting of electricity and natural gas utility networks, a multi-objective optimization method is used to solve the coordinated operation of the coupling network. In (Sun et al., 2020), authors study the day-ahead scheduling of gas-electric integrated energy system considering the bi-directional energy flow. The goal is to minimize the operation cost, and a second-order cone programming method is utilized to solve the problem. In (Fang et al., 2018), authors study the operation of the integrated gas and electrical power system considering the different response times of the gas and power systems. The problem is transformed into a single-stage linear programming. In (Chen et al., 2017), authors study the optimal operation of electricity-gas integrated energy system. The goal is to minimize the operation costs for both electrical and natural gas systems while satisfying steady-state operational constraints.

To model the electricity/heat/gas utility grid networks. The steady-state operational equations are often built as the constraints, and added to the previous optimization problem. For example, in (Liu et al., 2020), authors present a sequential reliability assessment method considering multi-energy flow and thermal inertia. Hydraulic circulation and heat exchange equations are used to model the thermal network. Conventional power flow equations are adopted to describe distribution network model. In (Martínez Ceseña et al., 2020), the electricity network model is represented as conventional power flow equations, as well as thermal and voltage limits. The gas network is represented as steady-state equations. The conventional steady-state equations and a thermal module are utilized to model the heat network. In (Yang et al., 2020), authors present a planning strategy for a district energy sector considering the coupling of power, gas, and heat systems. An optimal multi-energy flow model is developed, and the objective is to minimize operational costs. Distflow equations are used to describe the power distribution system, steady gas flow equations are adopted to model the gas distribution system, steady-state model is deployed to describe the distribution heat system. In (Martínez Ceseña and Mancarella, 2019), authors present a robust optimization framework for smart districts with multi-energy devices and electricity/heat/gas energy networks. The electricity network is modelled with typical power flow equations. The heat network is described based on nodal balance and cumulative head losses equations. The gas network is represented based on nodal balance, pressure drops, and head losses equations.

Based on the above reviews, optimization method is often used to calculate the power flow of the electricity/heat/gas energy networks. The electricity network is modelled based on typical power flow equations. The heat network is modelled based on nodal balance and heat losses equations. The gas network is represented based on nodal balance, pressure drops equations.

The above review shows that the operation problem of multi-energy supply microgrid and the operation problem of coupled electricity/heat/gas energy networks have drawn a lot of attention. However, using the ADP algorithm to solve the dispatching problem of the hydrogen-based multi-energy supply microgrids considering electricity/heat/gas energy networks has not drawn a lot of attention. The complexity of the whole model increases the difficulty of the control, especially the large numbers of constraints. Motivated by the aforementioned references, we present an ADP-based computationally efficient algorithm for the real-time operation of multi-energy supply grid-connected microgrids. A similar study is our previous work (Li et al., 2018a), in which only MPC algorithm is used, no other algorithms are compared.

Compared to previous works, the contribution of this paper can be concluded as follows:

• First, we build an ADP-based one-step decision model for the optimal operation of multi-energy supply grid-connected microgrids. In the one-step decision model, we consider large numbers of logical and physical constraints, and formulate the problem as a mixed-integer programming model;

The remainder of this paper is organized as follows. Section 2 describes the microgrid scheduling problem. Section 3 describes the electricity/heat/gas utility grids model. Section 6 presents the simulation results. Finally, Section 7 concludes the paper.

In fact, to operate the electricity/heat/gas integrated microgrids system, three aspects should be considered: 1) scheduling of the grid-connected microgrid; 2) utility grids operation; 3) the operation of the whole system.

In order to make the problem more readable, we use the simple model to describe the problem, and the detailed model is attached in Supplementary Material. The scheduling problem can be described as follows: m i n x t , x t + 1 , … , x t + T     ∑ τ = t t + T f ( x τ ) s . t .     A x i ≤ b ;   B x i = c   ( continuous variables )   l b ≤ x i ≤ u b   C x j ≤ d ;   D x j = e   ( integer / logical variables )   x i ∈ Z ;   x j ∈ { 0,1 , integer } (1) where x i are the continuous variables, x j are the integer/logical variables; A , B , C , D , b , c , d , e are the constraints matrix; f ( . ) is the operation cost function; T is the time horizon.

By solving the above mixed integer programming problem, we can obtain the scheduling results. However, due to the uncertainty of the load demand and the output of renewable energy, some parameters in constraints are not deterministic parameters. The above problem is then transferred to the following problem: m i n x t , x t + 1 , … , x t + T     ∑ τ = t t + T f ( x τ ) s . t .     A x i ≤ b ;   B x i = c ˜   ( continuous variables )   l b ≤ x i ≤ u b   C x j ≤ d ;   D x j = e   ( integer / logical variables )   x i ∈ Z ;   x j ∈ { 0,1 , integer } (2) where c ˜ are the uncertainty parameters. For example, in power balance constraints, generated power must equal to load demand, but the predicted load demand is uncertain.

The common method to solve the above uncertainty problem is stochastic optimization. The above problem can be transferred as follows: m i n x t 1 , x t + 1 1 , … , x t + T 1 x t 2 , x t + 1 2 , … , x t + T 2 … x t N s , x t + 1 N s , … , x t + T N s     ∑ s = 1 N s p s ⋅ ∑ τ = t t + T f ( x τ s ) s . t .     A x i s ≤ b ;   B x i s = c s ˜   ( continuous variables )   l b ≤ x i s ≤ u b   C x j s ≤ d ;   D x j s = e   ( integer / logical variables )   x i s ∈ Z ;   x j s ∈ { 0,1 , integer }   s = 1,2 , … , N s (3)

In the above stochastic problem, we use a scenario-based method to transfer the uncertainty parameters c ˜ to typical scenarios N s , and the probability of each scenario is p s . Lastly, to solve the above problem, we can obtain the scheduling results in each scenario.

Assume that the variables that exchanged energy with utility grids are x ex ∈ x i . Then the expected exchanged energy is: x e x ∗ = ∑ s = 1 N s p s ⋅ x e x s ,   s = 1,2 , … , N s (4)

Based on the day-ahead scheduling results, we can then implement real-time dispatching. Due to the real-time short-term prediction uncertainty, the real-time exchanged energy with the utility grid may not equal to the day-ahead scheduling results. In order to reduce this error, it is necessary for the real-time exchanged energy to follow the day-ahead scheduling results as close as possible. The sliding window model predictive control method is then adopted to deploy the real-time dispatching, the detailed model is attached in Supplementary Material. The real-time dispatching problem can be described as follows: m i n x t , x t + 1 , … , x t + t n     ∑ τ = t t + t n g ( x τ , x e x ∗ ) s . t .     A x i ≤ b ;   B x i = c   ( continuous variables )   l b ≤ x i ≤ u b   C x j ≤ d ;   D x j = e   ( integer / logical variables )   x i ∈ Z ;   x j ∈ { 0,1 , integer } (5) where g ( . ) is the real-time operation cost function; x e x ∗ are the day-ahead scheduling results; t n is the time horizon.

In the real-time sliding window dispatching, in the first time step t, the MPC problem is solved, then only the current time decisions (current time is t) are deployed, and the future decisions (future times are t + 1 , … , t + t n ) are abandoned. After that, the time slides to the next step t + 1 , and the MPC problem is solved again, then only the new current time decisions (new current time is t + 1 ) are deployed, and the future decisions (future times are t + 2 , … , t + t n + 1 ) are abandoned. This process is repeated until the last time is reached, the process can be seen in Figure 3A.

In the above section, the sliding window MPC method is adopted to deploy real-time dispatching. However, the solving time of the MPC method is long, because we need to solve the multiple windows optimization problem. In this section, the one-step decision model is developed to solve the real-time dispatching problem. With the one-step decision model, the solving time can be reduced. On the other hand, the ADP idea is also adopted, namely, integrating the future impacts into the current decision model, to make the current decision results more reasonable and effective.

In fact, the above MPC problem can be transferred into a series of smaller problems based on dynamic programming idea, which can be represented as follows: m i n x t     [ g ( x t , x e x t ) + m i n x t + 1   [ g ( x t + 1 , x e x t + 1 ) + m i n x t + 2   [ g ( x t + 2 , x e x t + 2 )   + ⋯ ∑ τ = t + t n − m t + t n g ( x τ , x e x ∗ ) ] ] ] s . t .     A x i ≤ b ;   B x i = c   ( continuous variables )   l b ≤ x i ≤ u b   C x j ≤ d ;   D x j = e   ( integer / logical variables )   x i ∈ Z ;   x j ∈ { 0,1 , integer } (6) We use V F t + 1 to describe the total future cost from t + 1 to t + t n , namely, V F t + 1 = m i n x t + 1 , … , x t + t n   ∑ τ = t + 1 t + t n g ( x τ , x e x ∗ ) (7) Then the above problem can be represented as: m i n x t     [ g ( x t , x e x t ) + V F t + 1 ] s . t .     A x i ≤ b ;   B x i = c   ( continuous variables )   l b ≤ x i ≤ u b   C x j ≤ d ;   D x j = e   ( integer / logical variables )   x i ∈ Z ;   x j ∈ { 0,1 , integer } (8) Because the future cost V F t + 1 is dependent on the current decisions x t and post-decision states S t + 1 , then the general one-step ADP decision model can be described as follows, and the detailed model is attached in Supplementary Material: m i n x t     g ( x t , x e x t , V F ( S t + 1 ) ) s . t .     A x i ≤ b ;   B x i = c   ( continuous variables )   l b ≤ x i ≤ u b   C x j ≤ d ;   D x j = e   ( integer / logical variables )   x i ∈ Z ;   x j ∈ { 0,1 , integer }   x i , x j   ∈ x t ;   S t + 1 = S F ( S t , x t ) ; (9) where V F is the approximate value function, V F ( S t + 1 ) is the approximate future operation cost based on the state S t + 1 ; S F is the state transition function, which is used to describe how the current state S t is changed to the next time state S t + 1 .

By solving the above one-step decision model (namely, the decision variables are only at the current time), one can obtain the optimal dispatching results. However, it can be seen that the main thing in the above one-step decision model is the approximate value function V F . If we can find a good approximate value function V F to describe the relationship between the state S t + 1 and the future operation cost C f u t u r e , then we can obtain good and effective decision results.

For the electricity utility grid operation, it is a classical optimal power flow (OPF) problem. The OPF problem can be seen as follows: m i n P g , Q g , V ∑ i = 1 n g { f P i ( P g i ) + f Q i ( Q g i ) } s . t . ( 12 ) , ( 13 ) (11) where the P g i , Q g i are the real and reactive power of the i t h generator. f P i , f Q i are the individual polynomial cost function of the i t h generator.

Power balance constraints can be shown as the following: P i g = P i l o a d + ∑ j = 1 n b u s V i V j ( G i j l i n e c o s θ i j + B i j l i n e s i n θ i j ) Q i g = Q i l o a d + ∑ j = 1 n b u s V i V j ( G i j l i n e s i n θ i j − B i j l i n e c o s θ i j ) (12) where P i l o a d , Q i l o a d are the real and reactive load demand at bus i. G i j l i n e , B i j l i n e are the parameters of the power lines from bus i to bus j. V m i , m i n ≤ V m i ≤ V m i , m a x ; i = 1,2 , … , n b u s P g i , m i n ≤ P g i ≤ P g i , m a x ; i = 1,2 , … , n g Q g i , m i n ≤ Q g i ≤ Q g i , m a x ; i = 1,2 , … , n g (13) where V m i , V m i , m i n , V m i , m a x are the voltage magnitude, minimum voltage magnitude and maximum voltage magnitude at bus i. P g i , m i n , P g i , m a x , Q g i , m i n , Q g i , m a x are the minimum and maximum real and reactive power of i generator.

For the heating utility grid, it is a heating power flow problem. During the heating transportation, heat transportation loss should be considered. The heating transportation loss can be described as follows (Pirouti, 2013; Shabanpour-Haghighi and Seifi, 2015). Q h e a t l o s s = c p ⋅ m ˙ ( T s 1 − T s 2 ) (14) where c p is the specific heat capacity (KJ/kgK), m ˙ is the mass flow rate (kg/s), and T s 1 , T s 2 are the temperature at node s 1 and node s 2 .

The temperature drop through the heating flow system can be described as: T s 2 = ( T s 1 − T g ) ⋅ e − l U c p ⋅ m ˙ + T g (15) where l is the pipe length, U is the heat transition coefficient (W/mK), and T g is the ground temperature.

Based on (Eqs. 14 Eqs. 15), it can be seen that the heating loss during the transportation is a nonlinear equation. In order to reduce the complexity, in this paper, we choose a linear model to describe the heating transportation loss. We assume that the heating loss is a linear function of the transportation distance, which can be shown as the following: Q h e a t l o s s = k h e a t l o s s ⋅ l (16) where k h e a t l o s s is the coefficient of the heating loss.

Then, the heating power flow of the heating utility grid can be presented. For each heating pipeline, two state variables (binary variables, 0 or 1): U L i n e h e a t o u t , U L i n e h e a t i n are defined. Then the heating power flow in each pipeline can be described as the following constraints: 0 ≤ L i n e h e a t o u t ( i , t ) ≤ U L i n e h e a t o u t ( i , t ) ⋅ L i n e h e a t m a x ( i ) 0 ≤ L i n e h e a t i n ( i , t ) ≤ U L i n e h e a t i n ( i , t ) ⋅ L i n e h e a t m a x ( i ) U L i n e h e a t o u t ( i , t ) + U L i n e h e a t i n ( i , t ) ≤ 1 (17)

An example is presented here to explain the logical illustrated in Eq. 18. In Eq. 18, there are three nodes h 1 , h 2 , and h 3 , the connections are h 1 ↔ h 2 , and h 2 ↔ h 3 . The heating power flow at node h 2 can be described as in Eq. 19. h 1 L i n e h e a t i n ( 1 , t ) L i n e h e a t o u t ( 1 , t ) − − − ← → h 2 L i n e h e a t i n ( 2 , t ) L i n e h e a t o u t ( 2 , t ) − − − ← → h 3 L i n e h e a t i n ( 3 , t ) L i n e h e a t o u t ( 3 , t ) − − − ← → (18) L i n e h e a t o u t ( 1 , t ) ⋅ ( 1 − Q h e a t l o s s , 1 ) − L i n e h e a t i n ( 1 , t ) / ( 1 − Q h e a t l o s s , 1 ) = L i n e h e a t o u t ( 2 , t ) − L i n e h e a t i n ( 2 , t ) (19)

For the gas utility grid, it is a gas power flow problem. The gas flow can be described as follows (De Wolf and Smeers, 2000): s i g n ( f i j ) ⋅ f i j 2 = C i j 2 ( p i 2 − p j 2 ) (20) where f i j is the gas flow between nodes i and j, p i and p j are the pressure at nodes i and j, and C i j is a constant which depends on the length, the diameter and the absolute rugosity of the pipe and on the gas composition.

During the gas transportation, the pressure will drop, which is modeled as in Eq. 21. p 1 → f d e p − − − − f l o s s → f i n p 2 (21)

Based on Eqs. 20 Eqs. 21, we can obtain f d e p 2 = C 12 2 ( p 1 2 − p 2 2 ) .

Then, the gas pressure drop can be described as: f i n 2 = C 12 2 ( p 1 2 − p 2 2 − H l o s s 2 ) = C 12 2 ( p 1 2 − p 2 2 ) − C 12 2 H l o s s 2 = f d e p 2 − C 12 2 ⋅ H l o s s 2 (22)

Assume that the loss C 12 2 ⋅ H l o s s 2 can be represented as C 12 2 ⋅ H l o s s 2 ≈ f d e p 2 ⋅ f l o s s , where f l o s s is a coefficient parameter to describe the pressure drop. Next, we can obtain f i n 2 ≈ f d e p 2 − f d e p 2 ⋅ f l o s s , namely, f i n ≈ f d e p ( 1 − f l o s s ) .

In (Martinez-Mares and Fuerte-Esquivel, 2012), it shows that the pressure drop H l o s s is a complex function related to the nonlinear effect of the pipeline distance L g a s p i p e and the weather conditions. Coefficient parameter f l o s s is also a nonlinear function. In order to reduce the complexity, here a linear model is adopted to describe the pressure drop. Assume that coefficient parameter f l o s s is a linear function of the gas pipeline distance, which can be shown as f l o s s = k g a s l o s s ⋅ L g a s p i p e (23) where k g a s l o s s is the coefficient of the gas loss.

Then the gas power flow in the gas utility grid can be presented. For each pipeline, two state variables (binary variables, 0 or 1) U L i n e g a s o u t , U L i n e g a s i n are defined. Then the gas flow constraints are: 0 ≤ L i n e g a s o u t ( i , t ) ≤ U L i n e g a s o u t ( i , t ) ⋅ L i n e g a s m a x ( i ) 0 ≤ L i n e g a s i n ( i , t ) ≤ U L i n e g a s i n ( i , t ) ⋅ L i n e g a s m a x ( i ) U L i n e g a s o u t ( i , t ) + U L i n e g a s i n ( i , t ) ≤ 1 (24)

Here we also use an example to explain the gas flow, which is shown in Eq. 25. There are three nodes g 1 , g 2 , and g 3 . The connections are g 1 ↔ g 2 , and g 2 ↔ g 3 . The gas flow at node g 2 can be described as: g 1 L i n e g a s i n ( 1 , t ) L i n e g a s o u t ( 1 , t ) − − − ← → g 2 L i n e g a s i n ( 2 , t ) L i n e g a s o u t ( 2 , t ) − − − ← → g 3 L i n e g a s i n ( 3 , t ) L i n e g a s o u t ( 3 , t ) − − − ← → (25) L i n e g a s o u t ( 1 , t ) ⋅ ( 1 − f g a s l o s s , 1 ) − L i n e g a s i n ( 1 , t ) / ( 1 − f g a s l o s s , 1 ) = L i n e g a s o u t ( 2 , t ) − L i n e g a s i n ( 2 , t ) (26)

The gas flow in a gas pipeline is restricted by the pressure of the beginning and end nodes. This constraint can be described as: − C i j 2 ( p j , m a x 2 − p i , m i n 2 ) ≤ f i j ≤ C i j 2 ( p i , m a x 2 − p j , m i n 2 ) (27) where p i , m i n , p i , m a x , p j , m i n , p j , m a x are the minimum and maximum pressure at node i and j.

Different cases are presented to research about the performance of each algorithm. Cases A D P l i n e a r b and A D P l i n e a r c are used to study the linear regression AVF. Cases A D P n o n l i n e a r A and A D P n o n l i n e a r B are used to study the nonlinear regression AVF. Cases A D P o n l i n e 30 , A D P o n l i n e n e g 1 , and A D P o n l i n e n e g 3 are compared to study the online process. In order to study the influence of optimization window numbers, cases M P C 12 , M P C 6 , and M P C 1 are set. Cases A D P N N n e g 1 , A D P N N n e g 5 , and A D P N N n e g 10 are presented to study the neural network regression AVF. All cases are compared and concluded as follows:

1. A D P l i n e a r b : the AVF is constructed based on linear regression, and the coefficient is C b = 10 − 2 ;

The simulation results of the real-time SOC of MG4 can be seen in Figure 4. Here SOC means the percentage of hydrogen in tanks. It can be seen that with different algorithms, the real-time dispatching results are significantly different. This is because in different algorithms, the future operation value functions are different, leading to different scheduling results.

We compare these different algorithms in the following: •   A D P I n d :   m i n x t   u c o s t e l ⋅ | E g r i d e l , T − Z g r i d e l , t | + u c o s t h e a t ⋅ | E g r i d h e a t , T − Z g r i d h e a t , t |   + V F I n d ( S ( t + 1 ) , L p r e ) ⋅ C I n d   + α ⋅ L S g a s t ˜ + β ⋅ L S e l t ˜ + γ ⋅ L S h e a t t ˜ ; (29) where I n d = { l i n e a r , n o n l i n e a r , N N , o n l i n e } represents different types of ADP algorithms. E g r i d e l , T is the day-ahead exchanged electricity power at time T, Z g r i d e l , t is the real-time exchanged electricity power at time t, E g r i d h e a t , T is the day-ahead exchanged heat power at time T, Z g r i d h e a t , t is the real-time exchanged heat power at time t. | E g r i d e l , T − Z g r i d e l , t | is used to describe the real-time electricity power deviation from the day-ahead results, and the unit is MW. u c o s t e l , u c o s t h e a t are the unit cost of electricity and heat power deviation from the day-ahead results, and the unit is €/MW. L S k t ˜ , k = ( g a s , e l , h e a t ) are the load shedding of the gas, electric, and heat load demands, the unit is M W . α, β, γ are penalty values of demands load shedding, the unit is €/MW. •   V F l i n e a r ( S ( t + 1 ) , L p r e ) = a 0 + a 1 ⋅ S t + 1 + a 2 ⋅ L p r e ; (30) where C l i n e a r ∈ { C b , C c } are coefficients, which is used to adjust the proportion of linear based AVF. •     V F n o n l i n e a r ( S ( t + 1 ) , L p r e ) = b 0 + b 1 ⋅ S t + 1 + b 2 ⋅ L p r e   + b 3 ⋅ S t + 1 ⋅ L p r e + b 4 ⋅ S t + 1 2 + b 5 ⋅ L p r e 2 ; (31) where C n o n l i n e a r ∈ { C A , C B } are coefficients, which is used to adjust the proportion of nonlinear based AVF. •     V F N N ( S ( t + 1 ) , L p r e ) = N N ( S t + 1 , L p r e P V , L p r e e l , L p r e h e a t , L p r e g a s ) ; (32) where C N N ∈ { C N N A , C N N B , C N N C } are coefficients, which is used to adjust the proportion of neural network based AVF. •     V F o n l i n e ( S ( t + 1 ) , L p r e ) = a 0 ^ + a 1 ^ ⋅ S t + 1 + a 2 ^ ⋅ L p r e ; (33) where a 0 ^ , a 1 ^ , a 2 ^ are changed in each iteration. C o n l i n e ∈ { C o n l i n e A , C o n l i n e B } are coefficients, which is used to adjust the proportion of AVF. •   M P C :   m i n x t , … , x t + s w   ∑ j = 0 s w u c o s t e l ⋅ | E g r i d e l , T − Z g r i d e l , t + j |   + ∑ j = 0 s w u c o s t h e a t ⋅ | E g r i d h e a t , T − Z g r i d h e a t , t + j |   + α ⋅ ∑ j = 0 s w L S g a s t + j ˜ + β ⋅ ∑ j = 0 s w L S e l t + j ˜ + γ ⋅ ∑ j = 0 s w L S h e a t t + j ˜ •     G l o b a l : s w = 288 ; M P C 12 : s w = 12 ;   M P C 6 : s w = 6 ; M P C 1 : s w = 1 ; (34) where s w is the optimization window number in MPC algorithm.

In Figure 4, we set the case G l o b a l as the basic case, because in case G l o b a l , the scheduling results are “global optimization”; however, in the other cases, the results are “local optimization”. Compare cases A D P l i n e a r and case G l o b a l , the SOC curves are very different, especially, in cases A D P l i n e a r , the SOC value reaches at the maximum point. And it can also be seen that in cases A D P l i n e a r , different coefficient values C b , C c lead to different dispatching results. Compare cases A D P n o n l i n e a r and case G l o b a l , the SOC curves have a similar tendency, but the values are different. Compare cases M P C and case G l o b a l , it can be seen that with different optimization window numbers, the SOC curves are very different. For example, in case M P C 1 , it has a similar tendency SOC; however, in case M P C 12 , the SOC reaches at the minimum point.

We compare different A D P l i n e a r cases in Figure 5A, namely, we choose different coefficients A D P 1 : C a = 10 0 ; A D P n e g 2 : C b = 10 − 2 ; A D P p o s 2 : C c = 10 2 ; A D P n e g 4 : C d = 10 − 4 ; A D P p o s 4 : C e = 10 4 (35) The linear regression of value function is shown in Figure 5B.

In fact, case A D P 1 and case A D P n e g 2 have very similar SOC curve, and they overlap together. It can be seen that the scheduling results based on linear approximate value function A D P l i n e a r deviate from the “global optimization” curve. This means that the linear approximate value function can not describe the future operation cost well. One important reason is that the dataset which is used to regress the linear value function is not completely, the other reason is that the linear function can not regress the value function well, and at last, leading to inaccuracy approximate value function.

Then, we adopt the nonlinear function to regress the dataset. And we compare different A D P n o n l i n e a r cases in Figure 6A, namely, we choose different coefficients A D P n o n l i n e a r A : C A = 1.2 * 10 − 8 ; A D P n o n l i n e a r B : C B = 10 − 8 ; A D P n o n l i n e a r C : C C = 10 − 9 ; A D P n o n l i n e a r D : C D = 10 − 10 ; A D P n o n l i n e a r E : C e = 10 − 13 . The nonlinear regression of value function is shown in Figure 6B.

It can be seen that based on the nonlinear approximate value function, the scheduling results have similar tendency to the global results. And with different coefficients C A , C B , C C , C D , C E , the scheduling results are similar to each other. However, the SOC curve values are still far away from the Global optimization results.

After that we adopt the neural network to regress the dataset. And we compare different A D P N N cases in Figure 7A, namely, we choose different coefficients A D P N N A : C N N A = 10 − 1 ; A D P N N B : C N N B = 10 − 5 ; A D P N N C : C N N C = 10 − 10 .

It can be seen that based on the neural network approximate value function, if we choose the approximate coefficients, the scheduling results are very close to the global optimization results, which means that the neural network can regress the value function well.

After that we develop an online simulation process, namely, at each time, the one-step decision model is iteratively simulated 30 times. The simulated operation cost of MG4 is shown in Figure 7B.

At each time, the one-step optimization model is solved for 30 times, and in each iteration, the parameters of the approximate value function is updated. Based on Figure 7B, it can be seen that at each time step, after 30 times iteration, the operation costs keep constantly, which means that the iteration process is convergence.

In this section, we analyze the operation cost of MGs based on different algorithms. Operation costs are the results of the problem Eq. 29 and problem Eq. 34. We use a 2-norm error to describe the difference between real-time operation cost of different algorithms and global optimization. The 2-norm error can be represented as: e r r o r r c = ∑ t = 1 T | R C m e t h o d t − R C g l o b a l t | 2 (37) where e r r o r r c is the 2-norm error of real-time operation cost, R C m e t h o d t is the real-time operation cost under method m e t h o d = { A D P l i n e a r , A D P n o n l i n e a r , A D P N N , A D P o n l i n e , M P C } at time t, R C g l o b a l t is the real-time operation cost under global optimization at time t.