Optimal energy management of microgrid based on multi-parameter dynamic programming

Abstract

With the wide application of microgrid system, fluctuation and randomness are the characteristics of distributed generation output. The traditional energy management system can’t meet the requirements to ensure the security and stability of the grid. The microgrid energy management is of great significance to the stable operation of power grid. In order to obtain higher economic benefits and pay less environmental costs, reasonable scheduling of various distributed power sources is able to achieve this goal. In this article, microgrid energy management including distributed generation is studied. The objective function includes the economic objective and the environmental objective. The model of energy management is considered as a multi-objectives and multi-parametric optimization problem. The multi-parameter dynamic programming is used to optimize the energy management of microgrid. Finally, the efficiency of the proposed method is examined by the simulation studies.

Introduction

Microgrid is an effective form of distributed energy grid-connected, which becomes an important part of smart grid. ¹ The use of clean energy is encouraged by countries that further promote the development of microgrids. The microgrid mainly includes distributed generators (DGs), diesel generators, energy storage systems (ESSs), as well as AC and DC loads. The operation of the microgrid has two modes of operation: islanded and grid-connected mode. In the islanded mode, the output power of distributed power in the microgrid must meet the requirements of local loads. The microgrid is connected to the main network via a point of common coupling (PCC). In the grid-connected operation model, the microgrid is able to provide energy to the main network or absorb power from power network. ² The application of multiple distributed power sources, bidirectional power and information flow, energy storage and other equipment makes the microgrid more flexible and intelligent. ³ The two-way flow of energy and information is the core feature of the microgrid, which makes the complex power control and management of the microgrid possible. ⁴ Under the condition of satisfying the power balance, it is necessary to study the energy optimization management strategy of the microgrid for the economical operation and the coordinated control of the distributed power source, improving the efficiency of energy utilization and the optimal distribution of power.

The operation goal of the microgrid is to obtain minimum operating costs, as well as, high quality of power supply. However, A lot of distributed power will change the balance of supply and demand in the microgrid, especially in the islanded mode. Distributed generations also increase the voltage of the system and reduce the quality of voltage. In general, the energy management system (EMS) can satisfy the objective function as well as other operational constraints. In the previous research, many algorithms are applied to the energy management of microgrids. The impact of the service life of lead-acid batteries, and the independent operation cost of the microgrid and the maximum utilization of the battery are considered as the objective functions. ⁵ An optimization calculation program is designed for the energy management strategy. In Wang et al., ⁶ an online EMS that takes into account distributed hybrid fuel cells was introduced. The system structure consists of three layers: the first layer captures the possible modes of operation. The second layer is based on the energy distribution between the fuzzy controller of battery and the fuel cell. The third layer adjusts each subsystem. Mathematical model of EMS based on the reliability of battery capacity and the increase in load demand were presented to get maximize profits.^7,8 A novel networked EMS to distributed power generation scheduling schemes was introduced in Doostizadeh et al. ⁹ The energy exchange problem of microgrid is optimized to participate in the traditional electricity market to gain minimal global energy production costs. In Huang et al., ¹⁰ a method for coordinated control of ESSs and distributed power sources has been employed. Users can adjust electricity consumption according to market prices to meet economic requirements. At the same time, many algorithms are applied to the energy management of microgrids. The model and cost function of different distributed power sources are given in Chen at al., ¹¹ and the optimal solution is obtained by genetic algorithm and game theory. In Sortomme and EI-Sharkawi, ¹² a new effective algorithm is addressed to complete the optimal energy system of microgrids. In summary, the heuristic algorithm can only obtain the near-local optimal solution, and not the global optimal solution. Some of these algorithms depend on the initial value setting with low efficiency.

The purpose of this article is to propose a global optimal solution to solve the optimal energy management problem under different operation modes of microgrid. Considering the characteristics of different distributed generation and the randomness of load, an optimized scheduling scheme is proposed. The optimal energy management of microgrids is a multi-stage decision process problem. The multi-stage decision process is expressed as the dynamic process of successive stages in a discrete-time system to find the optimal way of each process to maximize the benefits or minimize the cost of each stage. Therefore, this article proposes a multi-parameter dynamic programming method to solve the optimal energy management in the grid, taking into the dynamic performance of the battery consideration. This algorithm can solve the problem of falling into local optimal solution. The simulation results show the feasibility of the optimization method proposed in this article.

In summary, this article consists of the following sections. Section ‘The optimization model of microgrid energy management’ describes objective function and limitation factors of microgrid energy management. This section defines the model of the microgrid energy system under different operating modes. Section ‘Multi-parameter dynamic programming’ presents characteristics of optimization algorithms and their applications in EMSs. Section ‘Simulation results’ evaluates the multi-parameter dynamic programming in the simulation experiments. Section ‘Conclusion’ concludes this article and summarizes the effectiveness and economics of the optimization algorithm.

The optimization model of microgrid energy management

In this article, the objective function and limitation factors of the energy management are suitable for two possible modes of operation in the microgrid, islanded and grid-connected mode. In the islanded mode, the objective function of microgrid is that minimizes energy output cost (fuel costs) and operation and maintenance costs of the microgrid. When is connected to the grid, the microgrid is able to supply or absorb power from the distribution network. When absorbing power, the objective function is to minimize energy consumption, manufacturing cost and the operation maintenance cost. The electricity purchase price is the smallest; when supplying power, the objective function is to obtain maximum profits, that is, energy sales revenue minus the cost of energy and operation maintenance cost. The EMS in microgrid is to work out the scheduling plan of each unit under the condition of satisfying the objective function and constraints. The scheduling plan can guarantee the economic and efficient operation of the unit. In this article, the dynamic performance of the battery is considered. Hence, a new dynamic programming algorithm is used to solve the problem.

Objective function

Some of the necessary information about 24-h energy management needs to be obtained ahead of time. Relevant information includes the following:

1. Predictive value of hourly load;

2. Hourly wind and photovoltaic power forecasting;

3. Grid price forecasting;

4. The cost and parameters of distributed power generation;

5. Power limits;

6. Initial charging conditions of the energy storage device.

The islanded operation mode

In the islanded mode, the internal DGs supply power to the load, which minimizes costs of production, the operation and maintenance cost. Its objective function as equation (1)

$$\min \sum _{\mathrm{t}=1}^{\mathrm{n}}\{\sum _{\mathrm{i}=1}^{\mathrm{m}}{\mathrm{F}}_{\mathrm{i}}({\mathrm{P}}_{\mathrm{i}}(\mathrm{t})){\mathrm{\tau }}_{\mathrm{i}}(\mathrm{t})+{\mathrm{S}}_{\mathrm{i}}(\mathrm{t})\}$$ (1)

where n is the number of time steps scheduled for the day, m is the number of schedulable DGs. ${\mathrm{P}}_{\mathrm{i}}(\mathrm{t})$ is the output power of the ith schedulable DG of t, ${\mathrm{F}}_{\mathrm{i}}({\mathrm{P}}_{\mathrm{i}}(\mathrm{t}))$ is corresponding function of fuel cost.

According to different structure of consumption energy, the fuel cost function of power supply can be summed up in two basic forms. For microturbine or fuel cell, its corresponding fuel cost function ${\mathrm{F}}_{\mathrm{i}}({\mathrm{P}}_{\mathrm{i}}(\mathrm{t}))$ as follows

$${\mathrm{F}}_{\mathrm{i}}({\mathrm{P}}_{\mathrm{i}}(\mathrm{t}))={\mathrm{b}}_{\mathrm{i}}{\mathrm{P}}_{\mathrm{i}}(\mathrm{t})+{\mathrm{c}}_{\mathrm{i}}$$ (2)

For diesel engines, the corresponding fuel cost function as follows

$${\mathrm{F}}_{\mathrm{i}}({\mathrm{P}}_{\mathrm{i}}(\mathrm{t}))={\mathrm{a}}_{\mathrm{i}}{\mathrm{P}}_{\mathrm{i}}(\mathrm{t}{)}^2+{\mathrm{b}}_{\mathrm{i}}{\mathrm{P}}_{\mathrm{i}}(\mathrm{t})+{\mathrm{c}}_{\mathrm{i}}$$ (3)

where ${\mathrm{a}}_{\mathrm{i}}$ , ${\mathrm{b}}_{\mathrm{i}}$ and ${\mathrm{c}}_{\mathrm{i}}$ are parameters of the fuel cost function

$${\mathrm{\tau }}_{\mathrm{i}}(\mathrm{t})=\{\begin{array}{c}1\mathrm{The}\mathrm{ith}\mathrm{power}\mathrm{supply}\mathrm{is}\mathrm{in}\mathrm{operation}\hfill \\ \multicolumn{1}{c}{0\mathrm{The}\mathrm{ith}\mathrm{power}\mathrm{supply}\mathrm{is}\mathrm{out}\mathrm{of}\mathrm{operation}}\end{array}$$ (4)

${\mathrm{S}}_{\mathrm{i}}(\mathrm{t})$ is the start-up cost function of the schedulable DG, the expression is

$${\mathrm{S}}_{\mathrm{i}}(\mathrm{t})=\{\begin{array}{cc}\hfill & \mathrm{s}{\mathrm{c}}_{\mathrm{i}}{\mathrm{\tau }}_{\mathrm{i}}(\mathrm{t})-{\mathrm{\tau }}_{\mathrm{i}}(\mathrm{t}-1)=1\hfill \\ \multicolumn{1}{c}{0}& \mathrm{otherwise}\hfill \end{array}$$ (5)

where $\mathrm{s}{\mathrm{c}}_{\mathrm{i}}$ is the starting cost of each unit.

As a result, the operation of power grid need to satisfy power balance mentioned below

$$\sum _{\mathrm{i}=1}^{\mathrm{m}}{\mathrm{P}}_{\mathrm{i}}(\mathrm{t})={\mathrm{P}}_{\mathrm{load}}(\mathrm{t})-{\mathrm{P}}_{\mathrm{wind}}(\mathrm{t})-{\mathrm{P}}_{\mathrm{pv}}(\mathrm{t})-{\mathrm{P}}_{\mathrm{ess}}(\mathrm{t})$$ (6)

where ${\mathrm{P}}_{\mathrm{load}}(\mathrm{t})$ is the predetermination of load power at time t.

The grid-connected operation mode

There are two types of operating scenarios for the microgrid in networked operation mode. Microgrid can both transmit power to the grid and absorb power. For the sake of clarity, the following are explained according to two different operating scenarios. Under the two operation modes, the state of power grid and the flow direction of energy are different. Therefore, two different objective functions will be obtained

Scenario 1: purchase electricity

In this condition, the objective function is given by

$$\min \sum _{\mathrm{t}=1}^{\mathrm{n}}\{\sum _{\mathrm{i}=1}^{\mathrm{m}}{\mathrm{F}}_{\mathrm{i}}({\mathrm{P}}_{\mathrm{i}}(\mathrm{t})){\mathrm{\tau }}_{\mathrm{i}}(\mathrm{t})+{\mathrm{S}}_{\mathrm{i}}(\mathrm{t})+{\mathrm{C}}_{\mathrm{grid}}(\mathrm{t}){\mathrm{P}}_{\mathrm{grid}}(\mathrm{t})\}$$ (7)

where ${\mathrm{C}}_{\mathrm{grid}}(\mathrm{t})$ is real-time price of the grid. ${\mathrm{P}}_{\mathrm{grid}}(\mathrm{t})$ is the power exchanged with the grid at time t, ${\mathrm{P}}_{\mathrm{grid}}(\mathrm{t})>0$ .

As a result, the operation of power grid need to satisfy power balance mentioned below

$$\sum _{\mathrm{i}=1}^{\mathrm{m}}{\mathrm{P}}_{\mathrm{i}}(\mathrm{t})={\mathrm{P}}_{\mathrm{load}}(\mathrm{t})-{\mathrm{P}}_{\mathrm{wind}}(\mathrm{t})-{\mathrm{P}}_{\mathrm{pv}}(\mathrm{t})-{\mathrm{P}}_{\mathrm{grid}}(\mathrm{t})-{\mathrm{P}}_{\mathrm{ess}}(\mathrm{t})$$ (8)

Scenario 2: sell electricity to the main grid

In this condition, the objective function is given by

$$\max \sum _{\mathrm{t}=1}^{\mathrm{n}}\{{\mathrm{S}}_{\mathrm{i}}(\mathrm{t})-\sum _{\mathrm{i}=1}^{\mathrm{m}}{\mathrm{F}}_{\mathrm{i}}({\mathrm{P}}_{\mathrm{i}}(\mathrm{t})){\mathrm{\tau }}_{\mathrm{i}}(\mathrm{t})-{\mathrm{C}}_{\mathrm{grid}}(\mathrm{t}){\mathrm{P}}_{\mathrm{grid}}(\mathrm{t})\}$$ (9)

where ${\mathrm{C}}_{\mathrm{grid}}(\mathrm{t})$ is real-time price of the grid. ${\mathrm{P}}_{\mathrm{grid}}(\mathrm{t})$ is the power exchanged with the grid at time t, ${\mathrm{P}}_{\mathrm{grid}}(\mathrm{t})<0$ .

As a result, the operation of power grid need to satisfy power balance mentioned below

$$\sum _{\mathrm{i}=1}^{\mathrm{m}}{\mathrm{P}}_{\mathrm{i}}(\mathrm{t})={\mathrm{P}}_{\mathrm{load}}(\mathrm{t})-{\mathrm{P}}_{\mathrm{wind}}(\mathrm{t})-{\mathrm{P}}_{\mathrm{pv}}(\mathrm{t})-{\mathrm{P}}_{\mathrm{grid}}(\mathrm{t})-{\mathrm{P}}_{\mathrm{ess}}(\mathrm{t})$$ (10)

Constraints

The various constraints defined for optimal energy management are formulated below

• 1. Power constraints for dispatchable DGs

$${\mathrm{P}}_{\mathrm{imin}}\le {\mathrm{P}}_{\mathrm{i}}(\mathrm{t})\le {\mathrm{P}}_{\mathrm{imax}}$$ (11)

• 2. Grid power exchange constraint

$${\mathrm{P}}_{\mathrm{gridmin}}\le {\mathrm{P}}_{\mathrm{grid}}(\mathrm{t})\le {\mathrm{P}}_{\mathrm{gridmax}}$$ (12)

• 3. ESS charging or discharging power constraint

$${\mathrm{P}}_{\mathrm{essmin}}\le {\mathrm{P}}_{\mathrm{ess}}(\mathrm{t})\le {\mathrm{P}}_{\mathrm{essmax}}$$ (13)

when the ESS is discharged, ${\mathrm{P}}_{\mathrm{ess}}(\mathrm{t})>0$ ; when the ESS is charged, ${\mathrm{P}}_{\mathrm{ess}}(\mathrm{t})<0$ .

• 4. Dynamic performance constraints of ESS

$$\mathrm{SOC}(\mathrm{t}+1)=\mathrm{SOC}(\mathrm{t})-\frac{\mathrm{\eta }{\mathrm{P}}_{\mathrm{ess}}(\mathrm{t})}{{\mathrm{W}}_{\mathrm{ess}}}$$ (14)

$$\mathrm{SO}{\mathrm{C}}_{\min }\le \mathrm{SOC}(\mathrm{t}+1)\le \mathrm{SO}{\mathrm{C}}_{\max }$$ (15)

where $\mathrm{SOC}(\mathrm{t})$ is the charge state of the ESS at time t. $\mathrm{\eta }$ is the charging or discharging efficiency of the energy storage unit. ${\mathrm{W}}_{\mathrm{ess}}$ is the capacity of the ESS. The decision variables that need to be determined are the charging/discharging power ${\mathrm{P}}_{\mathrm{ess}}(\mathrm{t})$ and their state of charges $\mathrm{SOC}(\mathrm{t})$ .

Multi-parameter dynamic programming

Battery is the most widely used energy storage decoration in microgrid. In the actual EMS, the charging level of the battery depends on the charging level of the adjacent time. Therefore, the battery charging level is related in the adjacent time. The optimal model of microgrid is considered as a dynamic programming problem. In the previous research, heuristic algorithm is used to solve the optimal solution problem. ¹³ However, the heuristic algorithm cannot guarantee the optimal solution. Therefore, multi-parameter dynamic programming is adopted in this article.

Considering different optimization goals is another important design in process. Not only multiple goals can be considered in technology selection and unit size determine part of the design process, but also in practice microgrid goals such as minimizing costs and emissions. ¹⁴ Dynamic programming algorithm can deal with optimization problem by numerical analysis method. Set ${\mathrm{x}}_{\mathrm{i}}$ as state variables which contain output powers of dispatchable DGs at time t. ${\mathrm{u}}_{\mathrm{i}}$ is the decision vector for each plan generation adjustment. The function ${\mathrm{g}}_{\mathrm{i}}$ shows the relationship between neighbouring states. The function ${\mathrm{\Gamma }}_{\mathrm{i}}$ represents the cost function at time t. In order to get minimum absorption cost, the following formula needs to be met

$$min\mathrm{\Gamma }=\sum _{\mathrm{i}=1}^{\mathrm{N}-1}{\mathrm{\Gamma }}_{\mathrm{i}}({\mathrm{x}}_{\mathrm{i}},{\mathrm{u}}_{\mathrm{i}})+{\mathrm{\Gamma }}_{\mathrm{N}}({\mathrm{x}}_{\mathrm{N}})$$ (16)

Subject to the following constraints

$$\{\begin{array}{c}{\mathrm{x}}_{\mathrm{i}+1}={\mathrm{g}}_{\mathrm{i}}({\mathrm{x}}_{\mathrm{i}},{\mathrm{u}}_{\mathrm{i}})\hfill \\ \multicolumn{1}{c}{{\mathrm{c}}_{\mathrm{iE}}({\mathrm{x}}_{\mathrm{i}},{\mathrm{u}}_{\mathrm{i}})=0,}\\ \multicolumn{1}{c}{{\mathrm{c}}_{\mathrm{iI}}({\mathrm{x}}_{\mathrm{i}},{\mathrm{u}}_{\mathrm{i}})\le 0,{\mathrm{c}}_{\mathrm{NI}}({\mathrm{x}}_{\mathrm{N}})\le 0}\end{array}{\mathrm{c}}_{\mathrm{NE}}({\mathrm{x}}_{\mathrm{N}})=0$$ (17)

where i = 1, 2, …, N − 1.

In discrete time coordinates, the optimal control variables ${\mathrm{u}}_{\mathrm{i}}(\mathrm{i}=1,2,\dots ,\mathrm{N}-1)$ as well as the initial state variable ${\mathrm{x}}_1$ are sought. ${\mathrm{g}}_{\mathrm{i}},{\mathrm{\Gamma }}_{\mathrm{i}}$ , ${\mathrm{c}}_{\mathrm{iE}}$ and ${\mathrm{c}}_{\mathrm{iI}}$ are second-order derivable functions. The objective function of optimal energy management in the islanded mode can be converted into

$${\mathrm{x}}_{\mathrm{i}}={[{\mathrm{P}}_1(\mathrm{i}){\mathrm{P}}_2(\mathrm{i})\cdots {\mathrm{P}}_{\mathrm{m}}(\mathrm{i}){\mathrm{P}}_{\mathrm{ess}}(\mathrm{i})\mathrm{SOC}(\mathrm{i})]}^{\mathrm{T}}$$ (18)

$${\mathrm{u}}_{\mathrm{i}}={[\mathrm{\Delta }{\mathrm{P}}_1(\mathrm{i})\mathrm{\Delta }{\mathrm{P}}_2(\mathrm{i})\cdots \mathrm{\Delta }{\mathrm{P}}_{\mathrm{m}}(\mathrm{i})\mathrm{\Delta }{\mathrm{P}}_{\mathrm{ess}}(\mathrm{i})]}^{\mathrm{T}}$$ (19)

$${\mathrm{\Gamma }}_{\mathrm{i}}({\mathrm{x}}_{\mathrm{i}},{\mathrm{u}}_{\mathrm{i}})={\mathrm{x}}_{\mathrm{i}}^{\mathrm{T}}{\mathrm{y}}_{\mathrm{i}}+{\mathrm{x}}_{\mathrm{i}}^{\mathrm{T}}{\mathrm{Q}}_{\mathrm{i}}{\mathrm{x}}_{\mathrm{i}}+{\mathrm{v}}_{\mathrm{i}}$$ (20)

$${\mathrm{\Gamma }}_{\mathrm{N}}({\mathrm{x}}_{\mathrm{i}})={\mathrm{x}}_{\mathrm{N}}^{\mathrm{T}}{\mathrm{y}}_{\mathrm{N}}+{\mathrm{x}}_{\mathrm{N}}^{\mathrm{T}}{\mathrm{Q}}_{\mathrm{N}}{\mathrm{x}}_{\mathrm{N}}+{\mathrm{v}}_{\mathrm{N}}$$ (21)

$${\mathrm{g}}_{\mathrm{i}}({\mathrm{x}}_{\mathrm{i}},{\mathrm{u}}_{\mathrm{i}})=\left[\begin{array}{cccccc}1\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{\cdots }& \cdots \hfill & \cdots \hfill & \cdots \hfill & \cdots \hfill & \cdots \hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \end{array}\right]{\mathrm{x}}_{\mathrm{i}}+\left[\begin{array}{ccccc}1\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & \cdots \hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 1\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill \end{array}\right]{\mathrm{u}}_{\mathrm{i}}$$ (22)

$$\begin{array}{cc}\multicolumn{1}{c}{{\mathrm{c}}_{\mathrm{iE}}({\mathrm{x}}_{\mathrm{i}},{\mathrm{u}}_{\mathrm{i}})}& =\left[\begin{array}{cccccc}1\hfill & 1\hfill & \cdots \hfill & \cdots \hfill & 1\hfill & 0\hfill \end{array}\right]{\mathrm{x}}_{\mathrm{i}}\hfill \\ \multicolumn{1}{c}{}& +[-{\mathrm{P}}_{\mathrm{load}}(\mathrm{i})+{\mathrm{P}}_{\mathrm{wind}}(\mathrm{i})+{\mathrm{P}}_{\mathrm{pv}}(\mathrm{i})]\hfill \end{array}$$ (23)

$$\begin{array}{cc}\multicolumn{1}{c}{{\mathrm{c}}_{\mathrm{iI}}({\mathrm{x}}_{\mathrm{i}},{\mathrm{u}}_{\mathrm{i}})}& =\left[\begin{array}{cccccccc}1\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{-1}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{\cdots }& \cdots \hfill & \cdots \hfill & \cdots \hfill & \cdots \hfill & \cdots \hfill & \cdots \hfill & \cdots \hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill & -1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & \frac{-{\mathrm{\eta }}_{\mathrm{i}}}{{\mathrm{W}}_{\mathrm{ess}}}\hfill & 1\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & \frac{{\mathrm{\eta }}_{\mathrm{i}}}{{\mathrm{W}}_{\mathrm{ess}}}\hfill & 1\hfill \end{array}\right]{\mathrm{x}}_{\mathrm{i}}+\left[\begin{array}{c}-{\mathrm{P}}_{1\max }\hfill \\ \multicolumn{1}{c}{{\mathrm{P}}_{1\min }}\\ \multicolumn{1}{c}{-{\mathrm{P}}_{2\max }}\\ \multicolumn{1}{c}{\cdots }\\ \multicolumn{1}{c}{-{\mathrm{P}}_{\mathrm{essmax}}}\\ \multicolumn{1}{c}{{\mathrm{P}}_{\mathrm{essmin}}}\\ \multicolumn{1}{c}{-\mathrm{SO}{\mathrm{C}}_{\max }}\\ \multicolumn{1}{c}{\mathrm{SO}{\mathrm{C}}_{\max }}\end{array}\right]\hfill \end{array}$$ (24)

$$\begin{array}{cc}\multicolumn{1}{c}{{\mathrm{c}}_{\mathrm{NI}}({\mathrm{x}}_{\mathrm{N}})=}& \left[\begin{array}{cccccccc}1\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{-1}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{\cdots }& \cdots \hfill & \cdots \hfill & \cdots \hfill & \cdots \hfill & \cdots \hfill & \cdots \hfill & \cdots \hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill & -1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & \frac{-{\mathrm{\eta }}_{\mathrm{N}}}{{\mathrm{W}}_{\mathrm{ess}}}\hfill & 1\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & \frac{{\mathrm{\eta }}_{\mathrm{N}}}{{\mathrm{W}}_{\mathrm{ess}}}\hfill & 1\hfill \end{array}\right]{\mathrm{x}}_{\mathrm{i}}\hfill \\ \multicolumn{1}{c}{}& +\left[\begin{array}{c}-{\mathrm{P}}_{1\max }\hfill \\ \multicolumn{1}{c}{{\mathrm{P}}_{1\min }}\\ \multicolumn{1}{c}{-{\mathrm{P}}_{2\max }}\\ \multicolumn{1}{c}{\cdots }\\ \multicolumn{1}{c}{-{\mathrm{P}}_{\mathrm{essmax}}}\\ \multicolumn{1}{c}{{\mathrm{P}}_{\mathrm{essmin}}}\\ \multicolumn{1}{c}{-\mathrm{SO}{\mathrm{C}}_{\max }}\\ \multicolumn{1}{c}{\mathrm{SO}{\mathrm{C}}_{\max }}\end{array}\right]\hfill \end{array}$$ (25)

where ${\mathrm{\eta }}_{\mathrm{N}}$ is η the ESS charging or discharging efficiency, ${\mathrm{W}}_{\mathrm{ess}}$ denotes the capacity of ESS. Both parameters can be obtained from the factory parameters of ESS.

In the network operation mode, it is essential to increase the number of state quantity and decision amount of power exchanged between the microgrid and the distribution network

$${\mathrm{x}}_{\mathrm{i}}={[{\mathrm{P}}_1(\mathrm{i}){\mathrm{P}}_2(\mathrm{i})\cdots {\mathrm{P}}_{\mathrm{m}}(\mathrm{i}){\mathrm{P}}_{\mathrm{ess}}(\mathrm{i}){\mathrm{P}}_{\mathrm{grid}}(\mathrm{i})\mathrm{SOC}(\mathrm{i})]}^{\mathrm{T}}$$ (26)

$${\mathrm{u}}_{\mathrm{i}}={[\mathrm{\Delta }{\mathrm{P}}_1(\mathrm{i})\mathrm{\Delta }{\mathrm{P}}_2(\mathrm{i})\cdots \mathrm{\Delta }{\mathrm{P}}_{\mathrm{m}}(\mathrm{i})\mathrm{\Delta }{\mathrm{P}}_{\mathrm{ess}}(\mathrm{i})\mathrm{\Delta }{\mathrm{P}}_{\mathrm{grid}}(\mathrm{i})]}^{\mathrm{T}}$$ (27)

The multi-stage decision process describes the dynamic process of successive phases in a discrete-time model, thereby getting the optimal path of each phase to maximize the benefits or the cost of the entire phase. In a multi-stage decision-making process, for a particular phase, it contains multiple system states to determine the current phase cost. A set of decision variables acts on the current state and produces a new state variable for the next phase. This process is repeated until the last stage. The final overall goal is to get the optimal path to the smallest sum of costs in all phase.

A classical multi-decision process including a model of discrete time and a convex phase additive cost function is expressed as

$${\mathrm{x}}_{\mathrm{k}+1}={\mathrm{f}}_{\mathrm{k}}({\mathrm{x}}_{\mathrm{k}},{\mathrm{u}}_{\mathrm{K}}),{\mathrm{x}}_{\mathrm{k}}\in \mathrm{X},{\mathrm{u}}_{\mathrm{k}}\in {\mathrm{U}}_{\mathrm{k}},\mathrm{k}\in \{0,1,\dots ,\mathrm{N}-1\}$$ (28)

$$\mathrm{J}(\mathrm{U})={\mathrm{g}}_{\mathrm{N}}({\mathrm{x}}_{\mathrm{N}})+\sum _{\mathrm{k}=0}^{\mathrm{N}-1}{\mathrm{g}}_{\mathrm{k}}({\mathrm{x}}_{\mathrm{k}},{\mathrm{u}}_{\mathrm{k}})$$ (29)

where ${\mathrm{x}}_{\mathrm{k}}$ is the state variable of k, $\mathrm{X}\in {\mathbb{R}}^{\mathrm{n}}$ ; ${\mathrm{u}}_{\mathrm{k}}$ is the optimal decision variable of k, $\mathrm{U}\stackrel{\mathrm{\Delta }}{=}\{{\mathrm{u}}_0,{\mathrm{u}}_1,\dots ,{\mathrm{u}}_{\mathrm{N}-1}\}$ , ${\mathrm{U}}_{\mathrm{k}}\subseteq {\mathbb{R}}^{\mathrm{m}}$ ; ${\mathrm{f}}_{\mathrm{k}}$ describes the dynamic characteristics; ${\mathrm{g}}_{\mathrm{k}}$ is the value function of k.

Using the dynamic programming algorithm, the optimization system is decomposed into a series of low-dimensional systems based on the optimization principle. ¹⁵ Hence, the complexity of the system can be greatly reduced. The objective function for an ordinary function multi-stage optimization problem can be formulated as

$${\mathrm{V}}_{\mathrm{k}}({\mathrm{x}}_{\mathrm{k}})=\underset{{\mathrm{u}}_{\mathrm{k}},\dots ,{\mathrm{u}}_{\mathrm{N}}}{min}[{\mathrm{g}}_{\mathrm{N}}({\mathrm{x}}_{\mathrm{N}})+\sum _{\mathrm{i}=\mathrm{k}}^{\mathrm{N}-1}{\mathrm{g}}_{\mathrm{i}}({\mathrm{u}}_{\mathrm{i}},{\mathrm{x}}_{\mathrm{i}})]$$ (30)

where ${\mathrm{u}}_{\mathrm{i}}={\mathrm{\mu }}_{\mathrm{i}}({\mathrm{x}}_{\mathrm{i}})\in {\mathcal{U}}_{\mathrm{i}}$ . Applying the optimization principle, the new recursive expressions can be given below

$${\mathrm{V}}_{\mathrm{k}}({\mathrm{x}}_{\mathrm{k}})=\underset{{\mathrm{u}}_{\mathrm{k}}\in {\mathcal{U}}_{\mathrm{k}}}{min}[{\mathrm{g}}_{\mathrm{k}}({\mathrm{u}}_{\mathrm{k}},{\mathrm{x}}_{\mathrm{k}})+{\mathrm{V}}_{\mathrm{k}+1}({\mathrm{x}}_{\mathrm{k}+1})]$$ (31)

It can be seen from the above formula that the existing cost function is total future cost functions, and the completed cost function and the set of cost functions homologous to the current decision. This method can solve the global optimal variable of the initial state. The obvious advantage is that based on the premise that all future phases are optimized to the existing phase, the decision of each phase is only the corresponding phase. Although dynamic programming is a mature method, there are some restrictions. This method can’t be used in many cases. The dynamic programming algorithm is applied to some problems with specific constraints, ¹³ ${\mathrm{u}}_0={\mathrm{K}}_0{\mathrm{x}}_0$ , ${\mathrm{u}}_1={\mathrm{K}}_1{\mathrm{x}}_1$ , …, ${\mathrm{u}}_{\mathrm{N}-1}={\mathrm{K}}_{\mathrm{N}-1}{\mathrm{x}}_{\mathrm{N}-1}((-\mathrm{b}\pm \sqrt{{\mathrm{b}}^2-4\mathrm{ac}})/2\mathrm{a})$ , ${\mathrm{u}}_{\mathrm{k}}\in {\mathcal{U}}_{\mathrm{k}}$ , ${\mathrm{K}}_{\mathrm{i}}$ is a real matrix. Under the constraints of inequality, non-linear decision rules are generated. In order to ensure that the solution obtained is the global optimal solution, non-convexity is introduced. The objective function in every decision process can be defined as $\{{\mathrm{x}}_{\mathrm{k}},{\mathrm{u}}_{\mathrm{k}}\}$ . In order to get a concrete solution, it is necessary to obtain global optimization parameters.

This article proposes a new algorithm to solve constrained dynamic problems. This method can effectively avoid solving no global optimal solution. This method compromises the basic definitions of dynamic programming algorithms with multi-parameter programming, which can handle uncertain data in the model. The calculation process is shown in Figure 1.

Figure 1.

Multi-stage process.

In Figure 2, k is the index of the time, ${\mathrm{x}}_{\mathrm{k}}$ is the state of the system at time k. Based on a sequence of stage-wise optimal decisions, the system transforms from its original state into a final state.

Figure 2.

Time horizon.

The last decision-making stage is shown in Figure 1, the optimization problem of the nth stage as given below

$$\underset{{\mathrm{u}}_{\mathrm{N}-1}}{min}{\mathrm{J}}_{\mathrm{N}}({\mathrm{u}}_{\mathrm{N}-1},{\mathrm{x}}_{\mathrm{N}-1})={\mathrm{g}}_{\mathrm{N}}·({\mathrm{x}}_{\mathrm{N}})+{\mathrm{g}}_{\mathrm{N}-1}({\mathrm{x}}_{\mathrm{N}-1},{\mathrm{u}}_{\mathrm{N}-1})$$ (32)

$${\mathrm{u}}_{\mathrm{N}-1}^{\min }\le {\mathrm{u}}_{\mathrm{N}-1}\le {\mathrm{u}}_{\mathrm{N}-1}^{\max }$$ (33)

$${\mathrm{x}}_{\mathrm{N}}^{\min }\le {\mathrm{x}}_{\mathrm{N}}\le {\mathrm{x}}_{\mathrm{N}}^{\max }$$ (34)

$${\mathrm{x}}_{\mathrm{N}}={\mathrm{f}}_{\mathrm{N}-1}({\mathrm{x}}_{\mathrm{N}-1},{\mathrm{u}}_{\mathrm{N}-1})$$ (35)

where ${\mathrm{x}}_{\mathrm{N}}$ , ${\mathrm{x}}_{\mathrm{N}-1}\in {\mathbb{R}}^{\mathrm{n}}$ , ${\mathrm{u}}_{\mathrm{N}-1}\in {\mathbb{R}}^{\mathrm{m}}$ . According to Kuhn–Tucker (KKT) conditions

$$\mathrm{\Delta }\mathcal{L}({\mathrm{u}}_{\mathrm{N}-1},\mathrm{\lambda },{\mathrm{x}}_{\mathrm{N}-1})=0$$ (36)

$${\mathrm{\lambda }}_{\mathrm{i}}{\mathrm{\psi }}_{\mathrm{i}}({\mathrm{u}}_{\mathrm{N}-1},{\mathrm{x}}_{\mathrm{N}-1})=0,\forall \mathrm{i}=1,\dots ,\mathrm{p}$$ (37)

$${\omega }_{\mathrm{j}}({\mathrm{u}}_{\mathrm{N}-1},{\mathrm{x}}_{\mathrm{N}-1})=0,\forall \mathrm{j}=1,\dots ,\mathrm{q}$$ (38)

$$\begin{array}{cc}\multicolumn{1}{c}{\mathcal{L}={\mathrm{J}}_{\mathrm{N}}({\mathrm{u}}_{\mathrm{N}-1},{\mathrm{x}}_{\mathrm{N}-1})}& +\sum _{\mathrm{i}=1}^{\mathrm{p}}{\mathrm{\lambda }}_{\mathrm{i}}{\mathrm{\psi }}_{\mathrm{i}}({\mathrm{u}}_{\mathrm{N}-1},{\mathrm{x}}_{\mathrm{N}-1})\hfill \\ \multicolumn{1}{c}{}& +\sum _{\mathrm{j}=1}^{\mathrm{q}}{\mathrm{\mu }}_{\mathrm{j}}{\omega }_{\mathrm{j}}({\mathrm{u}}_{\mathrm{N}-1})\hfill \end{array}$$ (39)

where λ and μ are the Lagrange multiplier phasors of the inequality and the equality. Due to ${\mathrm{x}}_{\mathrm{N}-1}$ is a variable, the quantity of the above formula ${\mathrm{u}}_{\mathrm{N}-1}^{*}$ is not a point, but an optimization function ${\mathrm{u}}_{\mathrm{N}-1}^{*}={\mathrm{\mu }}_{\mathrm{N}-1}^{*}({\mathrm{x}}_{\mathrm{N}-1})$ . This optimal function is based on Theorem 1 below, with $\mathrm{x}={\mathrm{x}}_{\mathrm{N}-1}$ and $\mathrm{u}={\mathrm{u}}_{\mathrm{N}-1}$ .

• Theorem 1: Reliability definition

Set ${\mathrm{x}}_0$ as a parameter vector function, $({\mathrm{u}}_0,{\mathrm{\lambda }}_0,{\mathrm{\mu }}_0)$ be the KKT optimization condition, where ${\mathrm{\lambda }}_0$ is non-negative and ${\mathrm{u}}_0$ is the feasible point. The assumptions are as follows:

1. Strict Complementary Slackness (SCS) is established;

2. Constraint gradients are linearly independent;

3. Second-order sufficiency condition (SOSC) holds.

On the frontier of ${\mathrm{x}}_0$ , the existence of the unique continuous differential function $\mathrm{Z}(\mathrm{x})=[\mathrm{u}(\mathrm{x}),\mathrm{\lambda }(\mathrm{x}),\mathrm{\mu }(\mathrm{x})]$ , satisfying the above formula and $\mathrm{Z}({\mathrm{x}}_0)=[\mathrm{u}({\mathrm{x}}_0),\mathrm{\lambda }({\mathrm{x}}_0),\mathrm{\mu }({\mathrm{x}}_0)]$

$$\left(\begin{array}{c}\frac{\mathrm{du}({\mathrm{x}}_0)}{\mathrm{dx}}\hfill \\ \multicolumn{1}{c}{\frac{\mathrm{d}\mathrm{\lambda }({\mathrm{x}}_0)}{\mathrm{dx}}}\\ \multicolumn{1}{c}{\frac{\mathrm{d}\mathrm{\mu }({\mathrm{x}}_0)}{\mathrm{dx}}}\end{array}\right)=-({\mathrm{M}}_0{)}^{-1}{\mathrm{N}}_0$$ (40)

where ${\mathrm{M}}_0$ and ${\mathrm{N}}_0$ are Jacobian matrices of the system of equations

$${\mathrm{M}}_0=\left(\begin{array}{ccccccc}{\nabla }^2\mathcal{L}\hfill & \nabla {\mathrm{\psi }}_1\hfill & \cdots \hfill & \nabla {\mathrm{\psi }}_{\mathrm{p}}\hfill & \nabla {\omega }_1\hfill & \cdots \hfill & \nabla {\omega }_{\mathrm{q}}\hfill \\ \multicolumn{1}{c}{-{\mathrm{\lambda }}_1{\nabla }^{\mathrm{T}}{\mathrm{\psi }}_1}& -{\mathrm{\psi }}_1\hfill & \hfill & \hfill & \hfill & \hfill & \hfill \\ \multicolumn{1}{c}{⋮}& \hfill & \ddots \hfill & \hfill & \hfill & \hfill & \hfill \\ \multicolumn{1}{c}{-{\mathrm{\lambda }}_{\mathrm{p}}{\nabla }^{\mathrm{T}}{\mathrm{\psi }}_{\mathrm{p}}}& \hfill & \hfill & -{\mathrm{\psi }}_{\mathrm{p}}\hfill & \hfill & \hfill & \hfill \\ \multicolumn{1}{c}{{\nabla }^{\mathrm{T}}{\omega }_1}& \hfill & \hfill & \hfill & \hfill & \hfill & \hfill \\ \multicolumn{1}{c}{⋮}& \hfill & \hfill & \hfill & \hfill & \hfill & \hfill \\ \multicolumn{1}{c}{{\nabla }^{\mathrm{T}}{\omega }_{\mathrm{q}}}& \hfill & \hfill & \hfill & \hfill & \hfill & \hfill \end{array}\right)$$ (41)

$${\mathrm{N}}_0={({\nabla }_{\mathrm{xu}}^2\mathcal{L},-{\mathrm{\lambda }}_1{\nabla }_{\mathrm{x}}^{\mathrm{T}}{\mathrm{\psi }}_1,\dots ,-{\mathrm{\lambda }}_{\mathrm{p}}{\nabla }_{\mathrm{x}}^{\mathrm{T}}{\mathrm{\psi }}_{\mathrm{p}},{\nabla }_{\mathrm{x}}^{\mathrm{T}}{\omega }_1,\dots ,{\nabla }_{\mathrm{x}}^{\mathrm{T}}{\omega }_{\mathrm{q}})}^{\mathrm{T}}$$ (42)

• Corollary 1: In the field of ${\mathrm{x}}_0$ , the first-order function estimates of $\mathrm{Z}(\mathrm{x})=[\mathrm{u}(\mathrm{x}),\mathrm{\lambda }(\mathrm{x}),\mathrm{\mu }(\mathrm{x})]$ are as follows

$$\left[\begin{array}{c}\mathrm{u}(\mathrm{x})\hfill \\ \multicolumn{1}{c}{\mathrm{\lambda }(\mathrm{x})}\\ \multicolumn{1}{c}{\mathrm{\mu }(\mathrm{x})}\end{array}\right]=\left[\begin{array}{c}{\mathrm{u}}_0\hfill \\ \multicolumn{1}{c}{{\mathrm{\lambda }}_0}\\ \multicolumn{1}{c}{{\mathrm{\mu }}_0}\end{array}\right]+({\mathrm{M}}_0{)}^{-1}{\mathrm{N}}_0\mathrm{x}+\mathrm{o}(\parallel \mathrm{x}\parallel )$$ (43)

where ${(\mathrm{u}}_0,{\mathrm{\lambda }}_0,{\mathrm{\mu }}_0)=[{\mathrm{u}(\mathrm{x}}_0),{\lambda (\mathrm{x}}_0),{\mathrm{\mu }(\mathrm{x}}_0)]$ , ${\mathrm{M}}_0=\mathrm{M}({\mathrm{x}}_0)$ , ${\mathrm{N}}_0=\mathrm{N}({\mathrm{x}}_0)$ and $\varphi (\mathrm{x})=\mathrm{o}(\parallel \mathrm{x}\parallel )$ , when x → x₀, $((\varphi (\mathrm{x}))/(\parallel \mathrm{x}\parallel ))\to 0$ . Every exact linear optimal solution is limited by feasible solution and optimal condition. $\stackrel{⌣}{\psi }$ represents a stable constraint, and $\tilde{\mathrm{\lambda }}$ represents an unstable constraint

$$\{\begin{array}{c}\stackrel{⌣}{\psi }(\mathrm{u}({\mathrm{x}}_{\mathrm{N}-1}),{\mathrm{x}}_{\mathrm{N}-1})\le 0\\ \tilde{\lambda }({\mathrm{x}}_{\mathrm{N}-1})\ge 0\end{array}$$ (44)

Therefore, an accurate linear conditional function can be obtained

$${\mathrm{u}}_{\mathrm{N}-1}={\mathrm{\mu }}_{\mathrm{N}-1}^{*}({\mathrm{x}}_{\mathrm{N}-1})$$ (45)

$$\{\begin{array}{c}\begin{array}{cc}{\mathrm{u}}_{\mathrm{N}-1}={\mathrm{K}}_{\mathrm{N}-1}^1{\mathrm{x}}_{\mathrm{N}-1}+{\mathrm{C}}_{\mathrm{N}-1}^1\hfill & {\mathrm{x}}_{\mathrm{N}-1}\in {\mathrm{CR}}_{\mathrm{N}-1}^1\hfill \end{array}\hfill \\ \multicolumn{1}{c}{\begin{array}{cc}{\mathrm{u}}_{\mathrm{N}-1}={\mathrm{K}}_{\mathrm{N}-1}^2{\mathrm{x}}_{\mathrm{N}-1}+{\mathrm{C}}_{\mathrm{N}-1}^2\hfill & {\mathrm{x}}_{\mathrm{N}-1}\in {\mathrm{CR}}_{\mathrm{N}-1}^2\hfill \end{array}}\\ \multicolumn{1}{c}{}& \hfill & ⋮\hfill \\ \multicolumn{1}{c}{\begin{array}{cc}{\mathrm{u}}_{\mathrm{N}-1}={\mathrm{K}}_{\mathrm{N}-1}^{{\mathrm{L}}_{\mathrm{N}-1}}{\mathrm{x}}_{\mathrm{N}-1}+{\mathrm{C}}_{\mathrm{N}-1}^{{\mathrm{L}}_{\mathrm{N}-1}}\hfill & {\mathrm{x}}_{\mathrm{N}-1}\in {\mathrm{CR}}_{\mathrm{N}-1}^{{\mathrm{L}}_{\mathrm{N}-1}}\hfill \end{array}}\end{array}$$ (46)

where ${\mathrm{K}}_{\mathrm{N}-1}^{\mathrm{i}}$ and ${\mathrm{C}}_{\mathrm{N}-1}^{\mathrm{i}}$ are real matrices, ${\mathrm{CR}}_{\mathrm{N}-1}^{\mathrm{i}}\subset {\mathbb{R}}^{\mathrm{n}}$ . Considering the N − 1 phase, the corresponding optimization problem turns into

$$\begin{array}{cc}\multicolumn{1}{c}{\underset{{\mathrm{u}}_{\mathrm{N}-2}}{min}{\mathrm{J}}_{\mathrm{N}-1}({\mathrm{u}}_{\mathrm{N}-2},{\mathrm{x}}_{\mathrm{N}-2})}& ={\mathrm{g}}_{\mathrm{N}}({\mathrm{x}}_{\mathrm{N}})+{\mathrm{g}}_{\mathrm{N}-1}({\mathrm{x}}_{\mathrm{N}-1},{\mathrm{u}}_{\mathrm{N}-1})\hfill \\ \multicolumn{1}{c}{}& +{\mathrm{g}}_{\mathrm{N}-2}({\mathrm{x}}_{\mathrm{N}-2},{\mathrm{u}}_{\mathrm{N}-2})\hfill \end{array}$$ (47)

$$\{\begin{array}{c}{\mathrm{u}}_{\mathrm{N}-2}^{\min }\le {\mathrm{u}}_{\mathrm{N}-2}\le {\mathrm{u}}_{\mathrm{N}-2}^{\max }\hfill \\ \multicolumn{1}{c}{{\mathrm{x}}_{\mathrm{N}-1}^{\min }\le {\mathrm{x}}_{\mathrm{N}-1}\le {\mathrm{x}}_{\mathrm{N}-1}^{\max }}\end{array}$$ (48)

where ${\mathrm{x}}_{\mathrm{N}}={\mathrm{f}}_{\mathrm{N}-1}({\mathrm{x}}_{\mathrm{N}-1},{\mathrm{u}}_{\mathrm{N}-1})$ , ${\mathrm{x}}_{\mathrm{N}-1}={\mathrm{f}}_{\mathrm{N}-2}({\mathrm{x}}_{\mathrm{N}-2},{\mathrm{u}}_{\mathrm{N}-2})$ , ${\mathrm{x}}_{\mathrm{N}}$ , ${\mathrm{x}}_{\mathrm{N}-1}$ , ${\mathrm{x}}_{\mathrm{N}-2}\in {\mathbb{R}}^{\mathrm{n}}$ , ${\mathrm{u}}_{\mathrm{N}-1}$ , ${\mathrm{u}}_{\mathrm{N}-2}\in {\mathbb{R}}^{\mathrm{m}}$ .

In the traditional dynamic programming process, it is able to get complete model information ${\mathrm{x}}_{\mathrm{k}+1}=\mathrm{A}{\mathrm{x}}_{\mathrm{k}}+\mathrm{B}{\mathrm{u}}_{\mathrm{k}}$ . $\{{\mathrm{x}}_{\mathrm{N}-2},{\mathrm{u}}_{\mathrm{N}-2}\}$ represents an existing cost function. A and B are system constant matrices. Due to ${\mathrm{x}}_{\mathrm{N}-1}$ is a non-linear function, ${\mathrm{u}}_{\mathrm{N}-1}={\mathrm{\mu }}_{\mathrm{N}-1}^{*}({\mathrm{f}}_{\mathrm{N}-2}({\mathrm{x}}_{\mathrm{N}-2},{\mathrm{u}}_{\mathrm{N}-2}))$ , so ${\mathrm{x}}_{\mathrm{N}-2}$ and ${\mathrm{u}}_{\mathrm{N}-2}$ are also non-linear functions. The global optimization problem based on the above formula cannot be solved. ${\mathrm{u}}_{\mathrm{N}-1}$ , ${\mathrm{u}}_{\mathrm{N}-2}$ and ${\mathrm{x}}_{\mathrm{N}-2}$ are convex functions.

• Lemma 1: If a dynamic system is represented by a complex function and taking cost function as objective function that minimizes the stages, the recursive formula for the dynamic programming of the k-stage value function is given as

$${\mathrm{V}}_{\mathrm{k}}({\mathrm{x}}_{\mathrm{k}})=\underset{{\mathrm{u}}_{\mathrm{k}}\in {\mathcal{U}}_{\mathrm{k}}}{min}[{\mathrm{g}}_{\mathrm{k}}({\mathrm{u}}_{\mathrm{k}},{\mathrm{x}}_{\mathrm{k}})+{\mathrm{V}}_{\mathrm{k}+1}({\mathrm{x}}_{\mathrm{k}+1})]$$ (49)

The variables $({\mathrm{x}}_{\mathrm{k}},{\mathrm{u}}_{\mathrm{k}}^{*}({\mathrm{x}}_{\mathrm{k}}),{\mathrm{u}}_{\mathrm{k}+1}^{*}({\mathrm{x}}_{\mathrm{k}+1}),\dots ,{\mathrm{u}}_{\mathrm{N}-1}^{*}({\mathrm{x}}_{\mathrm{N}-1}))$ meet constraints in formulas (50) and (51)

$$\begin{array}{c}{\mathrm{V}}_{\mathrm{k}}({\mathrm{x}}_{\mathrm{k}},{\mathrm{u}}_{\mathrm{k}}^{*}\left({\mathrm{x}}_{\mathrm{k}}\right),{\mathrm{u}}_{\mathrm{k}+1}^{*}\left({\mathrm{x}}_{\mathrm{k}+1}\right),\dots ,{\mathrm{u}}_{\mathrm{N}-1}^{*}\left({\mathrm{x}}_{\mathrm{N}-1}\right))\hfill \\ \le {\mathrm{V}}_{\mathrm{k}}({\mathrm{x}}_{\mathrm{k}},{\mathrm{u}}_{\mathrm{k}}\left({\mathrm{x}}_{\mathrm{k}}\right),{\mathrm{u}}_{\mathrm{k}+1}^{*}\left({\mathrm{x}}_{\mathrm{k}+1}\right),\dots ,{\mathrm{u}}_{\mathrm{N}-1}^{*}\left({\mathrm{x}}_{\mathrm{N}-1}\right))\hfill \end{array}$$ (50)

$$\begin{array}{c}{\mathrm{V}}_{\mathrm{k}+1}({\mathrm{x}}_{\mathrm{k}},{\mathrm{u}}_{\mathrm{k}}^{*}\left({\mathrm{x}}_{\mathrm{k}}\right),{\mathrm{u}}_{\mathrm{k}+1}^{*}\left({\mathrm{x}}_{\mathrm{k}+1}\right),\dots ,{\mathrm{u}}_{\mathrm{N}-1}^{*}\left({\mathrm{x}}_{\mathrm{N}-1}\right))\hfill \\ \le {\mathrm{V}}_{\mathrm{k}+1}({\mathrm{x}}_{\mathrm{k}},{\mathrm{u}}_{\mathrm{k}}^{*}\left({\mathrm{x}}_{\mathrm{k}+1}\right),{\mathrm{u}}_{\mathrm{k}+1}\left({\mathrm{x}}_{\mathrm{k}}\right),\dots ,{\mathrm{u}}_{\mathrm{N}-1}^{*}\left({\mathrm{x}}_{\mathrm{N}-1}\right))\hfill \end{array}$$ (51)

where ${\mathrm{u}}_{\mathrm{i}}^{*}$ is the optimal solution at this stage, $\mathrm{i}\in \{\mathrm{k},\mathrm{k}+1,\dots ,\mathrm{N}-1\}$ . Since ${\mathrm{V}}_{\mathrm{k}}$ and ${\mathrm{V}}_{\mathrm{k}+1}$ are convex functions, their optimal solutions is the only one. Therefore, the optimal value of the k-phase in the function can be acquired by contrasting the optimal value of the kth phase and the (k+1)th phase

$${\mathrm{u}}_{\mathrm{N}-2}={\mathrm{\mu }}_{\mathrm{N}-2}^{*}({\mathrm{x}}_{\mathrm{N}-2},{\mathrm{x}}_{\mathrm{N}-1})$$ (52)

$$\{\begin{array}{c}\begin{array}{cc}{\mathrm{u}}_{\mathrm{N}-2}={\mathrm{K}}_{\mathrm{N}-2}^1{\mathrm{x}}_{\mathrm{N}-2}+{\mathrm{H}}_{\mathrm{N}-2}^1{\mathrm{u}}_{\mathrm{N}-1}+{\mathrm{C}}_{\mathrm{N}-2}^1\hfill & \hfill \\ \multicolumn{1}{c}{{\mathrm{x}}_{\mathrm{N}-2},{\mathrm{u}}_{\mathrm{N}-1}\in {\mathrm{CR}}_{\mathrm{N}-2}^1}\end{array}\hfill \\ \multicolumn{1}{c}{\begin{array}{cc}{\mathrm{u}}_{\mathrm{N}-2}={\mathrm{K}}_{\mathrm{N}-2}^2{\mathrm{x}}_{\mathrm{N}-2}+{\mathrm{H}}_{\mathrm{N}-2}^2{\mathrm{u}}_{\mathrm{N}-1}+{\mathrm{C}}_{\mathrm{N}-2}^2\hfill & \hfill \\ \multicolumn{1}{c}{{\mathrm{x}}_{\mathrm{N}-2},{\mathrm{u}}_{\mathrm{N}-1}\in {\mathrm{CR}}_{\mathrm{N}-2}^2}\end{array}}\\ \multicolumn{1}{c}{}& \hfill & ⋮\hfill \\ \multicolumn{1}{c}{\begin{array}{cc}{\mathrm{u}}_{\mathrm{N}-2}={\mathrm{K}}_{\mathrm{N}-2}^{{\mathrm{L}}_{\mathrm{N}-2}}{\mathrm{x}}_{\mathrm{N}-2}+{\mathrm{H}}_{\mathrm{N}-2}^{{\mathrm{L}}_{\mathrm{N}-2}}{\mathrm{u}}_{\mathrm{N}-1}+{\mathrm{C}}_{\mathrm{N}-2}^{{\mathrm{L}}_{\mathrm{N}-2}}\hfill & \hfill \\ \multicolumn{1}{c}{{\mathrm{x}}_{\mathrm{N}-2},{\mathrm{u}}_{\mathrm{N}-1}\in {\mathrm{CR}}_{\mathrm{N}-2}^{{\mathrm{L}}_{\mathrm{N}-2}}}\end{array}}\end{array}$$ (53)

where ${\mathrm{K}}_{\mathrm{N}-2}^{\mathrm{j}}{\mathrm{H}}_{\mathrm{N}-2}^{\mathrm{j}}$ and ${\mathrm{C}}_{\mathrm{N}-2}^2$ are real matrices. According to the specific optimization strategy, at this stage, this method cannot solve the global optimization problem. This article examines all possible combinations using the optimization strategy. In summary, the algorithm flow is as follows:

• Step 1. Initialization j = 1, in the Nth stage of solving the problem, this process is regarded as the multi-parameter optimization problem of the existing state space.

• Step 2. Set j ← j + 1, solve the (N – j) + 1 stage of the problem, which is considered as a multi-parameter optimization problem, the parameter ${\mathrm{x}}_{\mathrm{N}-\mathrm{j}}$ is the existing state space, and the future optimization control variables are ${\mathrm{u}}_{\mathrm{N}-\mathrm{j}+1},\dots ,{\mathrm{u}}_{\mathrm{N}-1}$ .

• Step 3. Calculate the optimal control amount at the sampling time j, compare ${\mathrm{u}}_{\mathrm{N}-\mathrm{j}+1}={\mathrm{\mu }}_{\mathrm{N}-\mathrm{j}+1}({\mathrm{u}}_{\mathrm{N}-\mathrm{j}+2},\dots ,{\mathrm{u}}_{\mathrm{N}-1},{\mathrm{x}}_{\mathrm{N}-\mathrm{j}+1})$ with ${\mathrm{u}}_{\mathrm{N}-\mathrm{j}}={\mathrm{f}}_{\mathrm{N}-\mathrm{j}}({\mathrm{u}}_{\mathrm{N}-\mathrm{j}+1},\dots ,{\mathrm{u}}_{\mathrm{N}-1},{\mathrm{x}}_{\mathrm{N}-\mathrm{j}})$ , and calculate ${\mathrm{u}}_{\mathrm{N}-\mathrm{j}}={\mathrm{\mu }}_{\mathrm{N}-\mathrm{j}}({\mathrm{x}}_{\mathrm{N}-\mathrm{j}})$ .

• Step 4. If j = N, the calculation stops. Otherwise, return to Step 1.

Simulation results

To illustrate the effectiveness of the proposed algorithm in energy management of microgrid, the following microgrid system is used for simulation calculation. As shown in the structure of microgrid below (Figure 3), ¹⁶ the rated voltage of microgrid is 400 V. Microgrid is connected to the 20-kV grid via a public connection point (PCC). The power supply in the microgrid is composed of micro gas turbines, fuel cells, wind turbines, photovoltaic power generation and loads.

Figure 3.

The structure of microgrid.

The parameters of each micro power supply are given in Table 1. The parameters of cost function are given in Table 2.

Table 1.

The parameters of each unit.

Unit type	Maximum power (kW)	Minimum power (kW)
Micro gas turbine	30	6
Fuel cell	50	6
Wind turbines	15	0
PV	13	0
Energy storage unit	30	−30

Table 2.

The parameters of cost function.

Unit type	b (Ect/kW h)	c (Ect/kW h)	Start-up cost
Micro gas turbine	4.37	85.06	9
Fuel cell	2.88	255.18	16

Gas emission cost parameters of distributed power (see Table 3):

Table 3.

The parameters of gas emission cost.

Gas type	External cost ($/1b)	Fuel cell emissions ($/1b)	Micro gas turbine ($/1b)
CO2	0.014	1.078	1.596
SO2	0.99	0.006	0.008
NOX	4.2	0.03	0.44

This article determines the time-of-use electricity price with reference to Wu et al. ¹⁴ (see Table 4):

Table 4.

The time-of-use electricity price.

Time	The purchase price of power (c$/kW h)	The selling price of power (c$/kW h)
(11 pm–7 am)	5.7323	5.9492
(7 am–10 am); (3 pm–6 pm); (9 pm–11 pm)	9.7385	5.9492
(10 am–3 pm); (6 pm–9 pm)	13.852	5.9492

The islanded mode

In the islanded mode, the predicted solar and wind power generation capacity are shown in Figures 4 and 5 below. The sum of the maximum output capacity of all power sources can meet the maximum load demand.

Figure 4.

Forecasted wind generation.

Figure 5.

Forecasted PV generation.

In isolated grid operation mode, the load forecast is shown in the following Figure 6. The sum of the maximum output capacities of all power sources can meet the maximum load demand.

Figure 6.

Load demand.

In the results in Figure 7, it can be seen that the charge and discharge of the energy storage units have a consistent trend with load change. When the load is small, the load is powered by distributed power and the battery pack is charged. When the load is large, the energy storage unit discharges, which will reduce fuel costs to meet economic operation requirements.

Figure 7.

Simulation results in islanded mode.

Purchasing electricity from the main network

In this operating mode, when the microgrid is connected to the network, the output of all power sources in the microgrid cannot meet requirements of the load. The load curve is shown in the following Figure 8.

Figure 8.

Load demand.

The result in Figure 9 calculated by the algorithm shows that at a time when electricity prices are low, the energy storage unit is charged. During the time when the electricity price is high, the energy storage unit inside the microgrid will be discharged, so the unit can reduce overall system operating costs.

Figure 9.

Simulation results in purchasing electricity mode.

Selling electricity to the main network

In this mode of operation, the power in the microgrid sends electrical energy to the main grid (see Figure 10).

Figure 10.

Load demand.

In Figure 11, it can be seen that the energy storage unit is charged when the electricity price is low. When the electricity price is high, the energy storage unit is discharged and delivering energy to the main network. The coordination of energy storage units and distributed power sources can improve energy efficiency to reduce operating costs.

Figure 11.

Simulation results in selling electricity mode.

In Table 5, genetic algorithm is applied to the optimal energy management model proposed in this article. Under the operation mode of island and network, the calculation time of the optimization algorithm proposed in this article is compared with that of genetic algorithm. In comparison, the efficiency of the multi-parameter dynamic programming can be verified.

Table 5.

Computation time for multi-parameter dynamic programming and genetic algorithm.

Optimization algorithm	Total computation time (s)
	Isolated mode	Network mode
Genetic algorithm	2.46	3.52
Multi-parameter dynamic programming	8.61	11.23

Conclusion

The purpose of this article is to study the optimal energy management of the microgrid in order to obtain the power generation plan of each unit in every hour of the next day. Under different operation modes of microgrid, battery performance and coordination with other distributed energy resource (DER) are considered in this article. Therefore, the multi-parameter dynamic programming algorithm proposed considers the distributed power and the fluctuations of load in the microgrid. The simulation results show that the method can guarantee the solution of the optimization problem. According to the energy management strategy, the fuel cost is minimized, and the gas emissions are reduced with high energy efficiency. The proposed method converges quickly and obtains the global optimal solution.