Grouped consensus power allocation method for multiple energy storage units

Abstract

A grouped consensus power allocation method for multiple energy storage units is provided, including determining a grouped coordination control strategy of multiple battery energy storage units (BESUs) of a BESS (battery energy storage system), where the multiple BESUs are communicated through a communication topology and grouped into a charging group and a discharging group, performing power coordination control power coordination control between the charging group and the discharging group based on a distributed algorithm, determining whether a switching between the charging group and the discharging group is triggered during an operation process of the BESS, and switching when it is triggered. The method can reduce the number of charging/discharging state conversions of the BESUs, prolong the service life of energy storage, ensure that the power allocation results meet the power constraints of the BESUs, reduce low-power operation times of the BESUs, and improve the efficiency of energy storage operation.

Description

FIELD OF THE DISCLOSURE

The disclosure relates to the field of power coordination control of multiple energy storage units, and more particularly to a grouped consensus power allocation method for multiple energy storage units.

BACKGROUND OF THE DISCLOSURE

With the increasingly prominent energy and environmental problems, new energy power generation technologies such as wind energy and solar energy have been rapidly developed. However, affected by climate, environment and other factors, new energy power generation has strong randomness, volatility and intermittency, which restricts the large-scale development and utilization of new energy to a great extent. It is a hot spot to study how to alleviate the adverse impact of new energy on the safe and stable operation of power system, so as to improve its permeability in power system. Battery energy storage has flexible source/load characteristics, and has the characteristics of high energy density, strong regulation ability, flexible allocation and relatively low cost. Using large-capacity battery energy storage to assist new energy power generation is one of effective technical solutions to solve the above problems. Large-capacity energy storage system is often composed of multiple energy storage units, and states of respective energy storage units are different. If the state differences among the respective energy storage units are ignored in the design of the power allocation strategy between the respective energy storage units, the operating conditions of the respective energy storage units will deteriorate, which will greatly deteriorate the regulation ability and operation economy of the energy storage system. Therefore, considering the characteristics of the respective energy storage units, it is of great significance to study how to refine and allocate the power between the respective energy storage units and improve the operating conditions of the respective energy storage units, so as to give full play to the regulation potential of energy storage and improve the operation economy of energy storage.

At present, scholars at home and abroad have done a lot of research on the coordinated power allocation strategy among multiple energy storage units, and the coordinated power allocation strategy can be divided into centralized control and distributed control in terms of control mode. Among them, the centralized control is difficult to meet the control requirements of power grid energy storage system under the high permeability of renewable energy. On the one hand, under the centralized control mode, a control center needs to establish communication with each of the energy storage units, the calculation amount of the control center is large, the calculation efficiency of the control center is low, there are risks such as control center failure and communication failure, and the system reliability is relatively low. On the other hand, under the background of large-scale grid connection of new energy and the gradual complexity and variability of energy structure in the power grid, the flexibility and scalability of the energy storage system are becoming more and more important, and the centralized control mode is difficult to meet the development needs. Based on multi-agent consensus theory, distributed coordinated control can realize the autonomous cooperative control among multiple agents in the system only through the limited data transmission between adjacent agents, eliminates the centralized control center, and has high robustness in dealing with communication changes or faults. Therefore, the distributed coordinated control has high flexibility and scalability and is an effective method to solve the above problems of centralized control, and the distributed coordinated control has been widely used in the field of power system. In addition, a large number of studies have shown that the use of grouping control for multiple energy storage units is conducive to improve the capacity utilization of the energy storage units and prolong the service life of battery energy storage. Therefore, a grouped consensus power allocation method for multiple energy storage units is proposed herein to solve the above problems.

SUMMARY OF THE DISCLOSURE

An objective of the disclosure is to provide a grouped consensus power allocation method for multiple energy storage units to solve the problems existing in the related art.

In order to achieve the above purpose, the disclosure provides a grouped consensus power allocation method for multiple energy storage units, including:

determining a grouped coordination control strategy of multiple battery energy storage units (BESUs), wherein a battery energy storage system (BESS) includes the multiple BESUs, the multiple BESUs are communicated through a communication topology, and the multiple BESUs are grouped into a charging group C_Z and a discharging group D_Z; and performing power coordination control between the charging group C_Z and the discharging group D_Z based on a distributed algorithm, and determining whether a switching between the charging group C_Z and the discharging group D_Z is triggered during an operation process of the BESS.

In an embodiment, the determining whether a switching between the charging group C_Z and the discharging group D_Z is triggered includes: switching charging/discharging states of the charging group C_Z and the discharging group D_Z, in response to a state of charge (SOC) of any one group of the BESUs reaches an upper threshold or a lower threshold.

In an embodiment, the charging group C_Z and the discharging group D_Z are defined as mutual dual groups. A BESU i of the charging group C_Z and a BESU i^∘ of the discharging group D_Z are mutual dual units when the BESU i and the BESU i^∘ meet a direct intergroup physical communication condition, where i ε {1, 2 . . . , n}, i^∘ ε {1, 2 . . . , n}, and n is the number of the BESUs in the charging group C_Z.

In an embodiment, the performing power coordination control between the charging group C_Z and the discharging group D_Z based on a distributed algorithm includes:

step S31, initializing the BESS; and

step S32, determining an adjustment strategy of an adjacency matrix A, including:

determining, by the BESU i, whether a dual group is needed to participate in a process of the iterating, and adjusting the adjacency matrix A through a formula (22) to achieve intergroup coordination in response to the dual group is needed to participate in the process of the iterating, where the formula (22) is expressed as follows:

$\begin{array}{cc}\{\begin{array}{cc}{a}_{i{i}^o}=1& x_i\left(k_c\right)\notin \left[{\underset{¯}{x}}_i,{\overline{x}}_i\right]\\ a_{i{i}^o}=0& x_i\left(k_c\right)\in \left[{\underset{¯}{x}}_i,{\overline{x}}_i\right]\end{array},& \left(22\right)\end{array}$

where x_i and x_i respectively represent a lower limit and an upper limit of a state constraint of the BESU i as a i-th node, x_i(k_c) represents a state of a consensus variable of the i-th node at a time k_c, a_ii∘ represents an element in the adjacency matrix A, a_ii∘=1 represents that the BESU i exchanges data with the BESU i^∘ being a dual unit of the BESU and a_ii∘=0 represents there is no data exchange;

determining whether the BESU i with an initial power of P_bi⁰=P_BS^ref(t) exits a communication iteration of a current control cycle, and adjusting the BESS as per a formula (25) in response to it is determined that the BESU i with the initial power exits the communication iteration, otherwise continuing to iterating without the adjusting until reaching a convergence accuracy, where the formula (25) is expressed as follows:

$\begin{array}{cc}\{\begin{array}{cc}{x}_j^{\prime }\left(k\right)=x_j\left(k\right)+\frac{x_i\left(k\right)w_i}{d_i\left(k-1\right)w_j}& \\ x_i^{\prime }\left(k\right)=0& \forall j\in N_i\left(k-1\right)\\ a_{\mathrm{ij}}=0& \\ d_j^{\prime }\left(k\right)=d_j\left(k\right)-1& \end{array},& \left(25\right)\end{array}$

where a_ij represents an element in the adjacency matrix A representing a communication network topology, x_i(k) represents a state of the i-th node at a time k, x_j(k) represents a state of a j-th node at the time k, w_i represents a weight coefficient of the BESU i, w_j represents a weight coefficient of a BESU j, N_i(k−1) represents an adjacent node set of the BESU i before exits, x′_j(k) represents a state quantity of the j-th node in N_i(k−1) after the BESU i exits, and d′_hd(k) represents the number of the adjacent nodes of the j-th node in N_i(k−1) after the BESU i exits; and

calculating a power P_bi(t) of each of the multiple BESUs as per a formula (28) in response to the convergence accuracy is reached, where the formula (28) is expressed as follows:

$\begin{array}{cc}\{\begin{array}{c}P\left(t\right)=W{X}^{*}\\ P\left(t\right)={\left[P_{b1}\left(t\right),P_{b2}\left(t\right),\cdots ,P_{\mathrm{bi}}\left(t\right)\cdots ,P_{\mathrm{bn}}\left(t\right)\right]}^T\\ P_B\left(t\right)={\sum }_{i=1}^n{P}_{\mathrm{bi}}\left(t\right)\end{array},& \left(28\right)\end{array}$

where X* represents a convergence state value meeting power constraints of all the n number of nodes, W=diag(w₁, w₂, . . . , w_n) represents a weight matrix of the BESS, P_B(t) represents a sum of powers of all the multiple BESUs, and P_bi(t) represents the power of the BESU i.

In an embodiment, the S31 includes:

1) defining a total output power of the BESS as P_BS^ref(t) and an initial power for iterating of the BESU i as P_bi⁰; where one of the BESUs of the charging group C_Z receives the total output power P_BS^ref(t) when P_BS^ref(t)≤0, and one of the BESUs of the discharging group D_Z receives the total output power P_BS^ref(t) when P_BS^ref(t)>0, thereby, for the BESU i received the total output power P_BS^ref(t), the initial power P_bi⁰ is expressed as follows:

$\begin{array}{cc}{P}_{\mathrm{bi}}^0=\{\begin{array}{ccc}{P}_{\mathrm{BS}}^{\mathrm{ref}}\left(t\right)& ,& i\in C_Z\mathrm{and}{P}_{\mathrm{BS}}^{\mathrm{ref}}\left(t\right)\le 0\\ P_{\mathrm{BS}}^{\mathrm{ref}}\left(t\right)& ,& i\in D_Z\mathrm{and}{P}_{\mathrm{BS}}^{\mathrm{ref}}\left(t\right)>0\\ 0& ,& \mathrm{other}\mathrm{conditions}\end{array},& \left(19\right)\end{array}$

where the BESU i belongs only to the charging group C_Z or the discharging group D_Z at a time t, and the initial power P_bi⁰ is 0 for the other BESUs that do not receive the total output power P_BS^ref(t).

2) defining an initial state quantity x_i⁰ for iterating, as per a formula (20) expressed as follows:
x_i⁰=w_i⁻¹P_bi⁰  (20),

where P_bi⁰ represents the initial power for iterating of each control cycle of the BESU i, and x_i⁰ represents an initial value of the consensus variable of each control cycle of the BESU i.

Therefore, the initializing the BESS for iterating is expressed as follows:

$\begin{array}{cc}\{\begin{array}{c}{P}_0\left(t\right)={\left[P_{b1}^0,P_{b2}^0,\cdots ,P_{\mathrm{bn}}^0\right]}^T\\ X_0=W^{-1}{P}_0\left(t\right)\end{array}.& \left(21\right)\end{array}$

In an embodiment, before the step S32, the method further includes: obtaining the adjacency matrix A, including:

making the charging group C_Z and the discharging group D_Z of the BESS at an initial time of each control cycle do not exchange information, that is, the BESU i and the dual unit BESU i^∘ thereof correspond to a_ii∘=0, while the charging group C_Z and the discharging group D_Z each are intragroup communicated, and a corresponding adjacency matrix is an initial adjacency matrix A₀ expressed as follows:

$A_0=\left[\begin{array}{cccccccc}0& 1& 0& 1& 0& 0& 0& 0\\ 1& 0& 1& 0& 0& 0& 0& 0\\ 0& 1& 0& 1& 0& 0& 0& 0\\ 1& 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& 0& 1\\ 0& 0& 0& 0& 1& 0& 1& 0\\ 0& 0& 0& 0& 0& 1& 0& 1\\ 0& 0& 0& 0& 1& 0& 1& 0\end{array}\right],\mathrm{and}$

expressing the adjacency matrix A of the BESS at the initial time of iterating of each control cycle as A=A₀, where the adjacency matrix A changes with iteration times k.

In an embodiment, the determining, by the BESU i, whether a dual group is needed to participate in a process of iterating includes:

defining the k_c as a positive integer, substituting X₀ into a formula (12) for iterating:
X(k+1)=M(k)X(k)+W⁻¹U(k)  (12),

where M (k) represents a non-negative row random matrix, U(k) represents a decoupling control quantity; X₀=[x₁⁰,x₂⁰, . . . ,x_n⁰]^T, where x_n⁰ is a state of an n-th BESU at the initial time; X(k)=[x₁(k),x₂(k), . . . ,x_n(k)]^T, where x_n(k) is a state of the n-th BESU at the time k; and X(k+1)=[x₁(k+1),x₂(k+1), . . . ,x_n(k+1)]^T, where x_n(k+1) is a state of the n-th BESU at a time k+1;

when a k_c-th time of iterating is ended, determining, by a BESU i, the total output power P_BS^ref(t) of the BESS at this time exceeds a total response capacity of the group where the BESU i is located in response to the state x_i(k_c) of the BESU i still does not meet the state constraint, requesting the dual group to participate in consumption and letting a_ii∘=1; otherwise letting a_ii∘=0 and continuing to iterating until reaching the convergence accuracy.

In an embodiment, the method further includes: setting the state constraint, including:

setting a state of the consensus variable of the i-th node at the time k as x_i(k) in considering of the BESS with n numbers of agent nodes as x_i(k); where the state constraint of the consensus variable of the i-th node is expressed as follows:
x_i(k)∈[x_i,x_i],

where [x,x] represents an intersection of constraints of all the n number of nodes, [x,x]=∩_i=1ⁿ[x_i,x_i], and x_i represents a power related quantity;

introducing a variable r_ij(k) on a basis of an ordinary weighted-consensus algorithm to realize that a distributed power coordinated allocation result meets power constraints of the respective BESUs, where an iterating and updating method of the consensus variable is expressed by a formula (1) as follows:
x_i(k+1)=x_i(k)+εw_i⁻¹Σ_j∈Ni(k)r_ij(k)a_ij·(s_j(k)−x_i(k)),∀i∈_n  (1),

where ε represents a convergence coefficient and ε>0; for an undirected graph, a_ij is expressed as a_ij=1 when the i-th node and the j-th node have direct data exchange, otherwise a_ij=0; N_i represents an adjacent node set of the i-th node, depending on the communication network topology and changing with a change of the communication network topology; _n represents a set of the all the n number of nodes of the BESS; and r_ij(k) represents an artificially introduced variable.

In an embodiment, the method further includes: obtaining the weight coefficient by a formula (3) expressed as follows:

$\begin{array}{cc}{w}_i=\frac{w_i^{\prime }n\epsilon }{{\min }_{\forall i\in {\mathbb{Z}}_n}\left\{w_i^{\prime }\right\}},& \left(3\right)\end{array}$

where w′_i is expressed as

$w_i^{\prime }=\{\begin{array}{c}{E}_{\mathrm{ba},i}·{\mathrm{SOC}}_B^{\max }-E_{b,i}\left(t\right),P_{\mathrm{BS}}^{\mathrm{ref}}\left(t\right)\le 0\\ E_{b,i}\left(t\right)-E_{\mathrm{ba},i}·{\mathrm{SOC}}_B^{\min },P_{\mathrm{BS}}^{\mathrm{ref}}\left(t\right)>0\end{array},$
E_ba,i represents a capacity of the BESU i, E_b,i(t) represents a residual capacity at the current time of the BESU i, n is the number of the BESUs, SOC_B^max represents an upper limit of a state of charge of each of the BESUs, and SOC_B^min represents a lower limit of the state of charge of each of the BESUs.

In an embodiment, the method further includes obtaining the introduced variable r_ij(k), including:

adjusting a constraint range by:

defining Δ_i^u(k) as a distance between the state x_i(k) and the upper limit x_i and defining Δ_i^l(k) as a distance between the state x_i(k) and the lower limit x_i, where Δ_i^u(k) and Δ_i^l(k) are respectively expressed as follows:

$\begin{array}{cc}\{\begin{array}{c}{\Delta }_i^u\left(k\right)=❘x_i\left(k\right)-{\overline{x}}_i❘\\ {\Delta }_i^l\left(k\right)=❘x_i\left(k\right)-{\underset{\_}{x}}_i❘\end{array},& \left(4\right)\end{array}$

adjusting the upper limit and the lower limit during iterating as per a formula (5) and a formula (6) in order to make a state outside of the constraint range be transferred to be within the constraint range in a process of iterating, where the formula (5) and the formula (6) are respectively expressed as follows:

$\begin{array}{cc}{{\overline{x}}_i}^{\prime }=\{\begin{array}{cc}{\overline{x}}_i-\tau {\Delta }_i^u\left(k\right)& x_i\left(k\right)>{\overline{x}}_i\\ {\overline{x}}_i& x_i\left(k\right)\le {\overline{x}}_i\end{array},& \left(5\right)\end{array}$ $\begin{array}{cc}{{\underset{\_}{x}}_i}^{\prime }=\{\begin{array}{cc}{\underset{\_}{x}}_i+{\mathrm{\tau \Delta }}_i^l\left(k\right)& x_i\left(k\right)<{\underset{\_}{x}}_i\\ {\underset{\_}{x}}_i& x_i\left(k\right)\ge {\underset{¯}{x}}_i\end{array},& \left(6\right)\end{array}$

the formula (5) and the formula (6) each represent that a constraint boundary is adjusted according to a distance between a current state and the constraint boundary, so that the state outside the constraint range is transferred to be within the state constraint in the process of iterating, and a value range of τ is:

$0<\tau <\min \left(\frac{❘{\overline{x}}_i❘}{{\Delta }_i^u\left(k\right)},\frac{❘{\underset{\_}{x}}_i❘}{{\Delta }_i^l\left(k\right)}\right);$

calculating intermediate variables as per a formula (7) expressed as follows:

$\begin{array}{cc}\{\begin{array}{c}\Delta {\overline{x}}_i\left(k\right)=❘x_i\left(k\right)-{{\overline{x}}_i}^{\prime }❘\\ \Delta {\underset{\_}{x}}_i\left(k\right)=❘x_i\left(k\right)-{{\underset{\_}{x}}_i}^{\prime }❘\end{array},& \left(7\right)\end{array}$

calculating the variable r_ij(k), including: sending, by the i-th node, the intermediate variables (Δx_i(k),Δx_i(k)) and the state x_i(k) to adjacent nodes, and receiving intermediate variables and a state of adjacent j-th node for calculating the variable r_ij(k) as per a formula (8) expressed as follows:

$\begin{array}{cc}{r}_{\mathrm{ij}}\left(k\right)=\{\begin{array}{cc}\min \left\{1,\frac{\Delta {\overline{x}}_i\left(k\right)}{❘{\sigma }_{\mathrm{ij}}\left(k\right)❘},\frac{\Delta {\underset{\_}{x}}_j\left(k\right)}{❘{\sigma }_{\mathrm{ij}}\left(k\right)❘}\right\}& {\sigma }_{\mathrm{ij}}\left(k\right)>0\\ \min \left\{1,\frac{\Delta {\overline{x}}_j\left(k\right)}{❘{\sigma }_{\mathrm{ij}}\left(k\right)❘},\frac{\Delta {\underset{¯}{x}}_i\left(k\right)}{❘{\sigma }_{\mathrm{ij}}\left(k\right)❘}\right\}& {\sigma }_{\mathrm{ij}}\left(k\right)<0\\ 1& {\sigma }_{\mathrm{ij}}\left(k\right)=0\end{array},& \left(8\right)\end{array}$

where σ_ij(k)=x_j(k)−x_i(k), j ε N_i(k); r_ij(k) is configured to prevent x_i(k) from exceeding the constraint boundaries in the process of iterating when transferring along a direction approaching the constraint boundaries to thereby ensure that a convergence result meets the state constraint.

The embodiments of the disclosure may mainly have beneficial effects as follows.

The BESS includes the multiple BESUs. The multiple BESUs are set to be communicated through the communication topology. The multiple BESUs are divided into the charging group C_Z and the discharging group D_Z, so that switching times of charging/discharging states of the multiple BESUs are reduced, and the service life of energy storage operation is prolonged. Using the distributed algorithm, the grouping distributed coordinated control of the BESS can reduce the switching times of charging/discharging states of the multiple BESUs and prolong the service life of energy storage operation. In addition, it can ensure that the power allocation results meet the power constraints of the multiple BESUs, reduce the times of low-power operation of the multiple BESUs and improve the efficiency of energy storage operation.

BRIEF DESCRIPTION OF DRAWINGS

The accompanying drawings, which are incorporated in and constitute a part of the disclosure, are intended to provide a further understanding of the disclosure. The illustrative embodiments of the disclosure and descriptions thereof are used to explain the disclosure and do not constitute undue restrictions on the disclosure. In the accompanying drawings:

FIG. 1 illustrates a schematic flowchart of a battery energy storage system (BESS) control according to an embodiment of the disclosure;

FIG. 2a illustrates a schematic diagram of an iterative convergence process of a classical weighted-consensus algorithm;

FIG. 2b illustrates a schematic diagram of an iterative convergence process of a target algorithm according to an embodiment of the disclosure;

FIG. 3 illustrates a schematic diagram showing iterative convergence processes obtained by battery energy storage units (BESUs) of the BESS using an adjustment strategy of an adjacency matrix A according to an embodiment of the disclosure;

FIG. 4 illustrates a schematic structural diagram of a grid-connected microgrid with wind charge storage and optical charge storage according to an embodiment of the disclosure;

FIG. 5 illustrates a schematic diagram showing regulation demand power curves of an energy storage system in a dispatching cycle of the grid-connected microgrid according to an embodiment of the disclosure;

FIG. 6 illustrates a schematic diagram showing power curves of a group Z1 and a group Z2 under a control strategy of according to an embodiment of the disclosure;

FIG. 7 illustrates a comparison diagram of grid-connected power volatility under different allocation strategies;

FIG. 8 illustrates a schematic diagram showing state of charge (SOC) change curves of the BESUs within 24 hours;

FIG. 9 illustrates a comparison diagram of switching times of charging/discharging states with different distribution strategies; and

FIG. 10 illustrates a topological diagram of the BESUs of the BESS of the disclosure.

DETAILED DESCRIPTION OF EMBODIMENTS

It should be noted that embodiments in the disclosure and features in the embodiments can be combined with each other without conflict. Hereafter, the disclosure will be described in detail below with reference to the accompanying drawings and in combination with the embodiments.

It should be noted that the steps shown in the flowcharts of the accompanying drawings may be exemplarily performed in a computer system with a set of computer executable instructions, and although a logical order is shown in the flowcharts, in some cases, the steps shown or described may be performed in a different order than here.

As shown in FIGS. 1-10, this embodiment provides a grouped consensus power allocation method for multiple energy storage units, including:

determining a grouped coordination control strategy of multiple battery energy storage units (BESUs); and

performing power coordination control by using a weighted-consensus algorithm based on distributed model predictive control (DMPC) and considering consensus variable constraints.

In an illustrated embodiment, the determining a grouped coordination control strategy of multiple BESUs may include steps as follows.

Due to the great uncertainty, randomness and volatility of new energy power generation, to effectively alleviate adverse impacts of the new energy power generation on a power system, it is necessary for an energy storage system to frequently switch between charging/discharging states, which accelerates the decline of energy storage life to a certain extent.

In order to reduce the switching times of charging/discharging states of the multiple BESUs, the multiple BESUs in a battery energy storage system (BESS) are divided into two groups Z1 and Z2. When one group is a charging group and the other group is a discharging group, the multiple BESUs are communicated with each other through a communication topology. The charging group is defined as C_Z and the discharging group is defined as D_Z. For the convenience of subsequent expression, it is defined the charging group C_Z and the discharging group D_Z as mutual dual groups. When a BESU i of the charging group C_Z and a BESU i^∘ of the discharging group D_Z meet a direct intergroup physical communication condition (as shown by broken lines in FIG. 10), the BESU i and the BESU i^∘ are mutual dual units, where i ε {1, 2 . . . , n}, and i^∘ε{1, 2 . . . , n}, and {1, 2 . . . , n} respectively represent serial numbers of the multiple BESUs. During operation, when a state of charge (SOC) of any one group of the BESUs reaches an upper threshold or a lower threshold, charging/discharging states of the charging group C_Z and the discharging group D_Z are switched, that is, the original charging group C_Z is converted to the discharging group D_Z and the original discharging group D_Z is converted to the charging group C_Z.

Research shows that the energy conversion efficiency of the BESU is closely related to the operating power. When the operating power per unit value is small, the loss of electrical equipment is large, which leads to the low energy conversion efficiency of the BESUs. When the operating power per unit value is large, the energy conversion efficiency is high. Therefore, the unnecessary energy loss caused by the low efficiency operation of the BESUs can be reduced, so as to improve the energy conversion efficiency.

As shown in FIG. 4, a distributed BESS is composed of multiple BESUs and sparse communication networks. Each of the multiple BESUs may be an agent with certain communication, calculation and control capabilities. In order to realize distributed coordination of the output of each of the multiple BESUs, ensure that power instructions of the multiple BESUs meet power constraints and a principle of “capable people do more work”, improve consensus of the multiple BESUs and improve a convergence speed of a distributed algorithm, the disclosure proposes a weighted-consensus algorithm based on DMPC considering consensus variable constraints. In this situation, each of the multiple BESUs is performed a limited data exchange with adjacent BESUs without a central node or a “leader” node, that is, an adaptive power allocation between the multiple BESUs can be completed in a very short time. Specifically, the algorithm is described as follows.

Step S1, considering the consensus variable constraints.

A distributed BESS with n numbers of agent nodes is considered. A state of a consensus variable of an i-th node at a time k is defined as x_i(k). A state constraint the consensus variable of the i-th node is expressed as follows:
x_i(k)ε[x_i,x_i],

where x_i and x_i respectively represent a lower limit and an upper limit of the state constraint of the BESU i as an i-th node, [x, x] represents an intersection of constraints of all the n numbers of nodes, [x, x]=∩_i=1ⁿ[x_i, x_i], and x_i represents a power related quantity. Specifically, x_i is a physical quantity related to weight and power, see a formula (1a) for details, and the upper and lower limits of the state constraint can be calculated according to a formula (1b) as follows.

In order to realize that coordinated allocation of distributed power meets power constraints of the multiple BESUs, a variable r_ij(k) is first introduced based on a classical weighted-consensus algorithm (also referred to as an ordinary weighted-consensus algorithm). An updating method of iterating is expressed by a formula (1) as follows:
x_i(k+1)=x_i(k)+εw_i⁻¹Σ_jεNi(k)r_ij(k)a_ij·(s_j(k)−x_i(k)),∀iε_n  (1),

where ε represents a convergence coefficient and ε>0; a_ij represents an element in an adjacency matrix A representing a communication network topology, for an undirected graph, a_ij=1 when the i-th node and a j-th node have direct data exchange, otherwise a_ij=0; N_i represents an adjacent node set of the i-th node, depending on the communication network topology and changing with a change of the communication network topology; _n represents a set of the all then number of nodes of the BESS; and w_i represents a weight coefficient of the BESU i. The weight coefficient is related to the SOC of the BESU. In the illustrated embodiment, a relationship among the weight coefficient w_i, the consensus variable x_i and a power P_b,i of the BESU i is established, as shown in the formula (1a). A calculation formula of an upper threshold and a lower threshold of a state of the consensus variable x_i of the BESU i in an operation cycle t is shown in a formula (1b). Specifically, the formula (1a) and the formula (1b) are respectively expressed as follows:

$\begin{array}{cc}{P}_{\mathrm{bi}}=w_i{x}_i,& \left(1a\right)\end{array}$ $\begin{array}{cc}\{\begin{array}{c}{\overline{x}}_i\left(t\right)=w_i^{-1}{P}_{b,i,\max }^{\mathrm{dis}}\\ {\underset{\_}{x}}_i\left(t\right)=w_i^{-1}{P}_{b,i,\max }^{\mathrm{ch}}\end{array},& \left(1b\right)\end{array}$

where P_b,i,max^ch and P_b,i,max^dis respectively represents a maximum allowable discharge power and a maximum allowable charging power in the operation cycle t of the BESU i.

In order to realize that the powers of the multiple BESUs meet the principle of “capable people do more work”, the weight coefficient is defined according to a capacity, a residual capacity and other information of the BESU as follows:

$\begin{array}{cc}{w}_i^{\prime }=\{\begin{array}{c}{E}_{\mathrm{ba},i}·{\mathrm{SOC}}_B^{\max }-E_{b,i}\left(t\right),P_{\mathrm{BS}}^{\mathrm{ref}}\left(t\right)\le 0\\ E_{b,i}\left(t\right)-E_{\mathrm{ba},i}·{\mathrm{SOC}}_B^{\min },P_{\mathrm{BS}}^{\mathrm{ref}}\left(t\right)>0\end{array},& \left(2\right)\end{array}$

where E_ba,i represents a capacity of the BESU i, and E_b,i(t) represents a residual capacity of the BESU i at a current time.

Due to the large battery capacity, when the BESU is discharged (charged), if the residual capacity of the BESU is high (low), and a value of w′_i will be large, while when the convergence coefficient ε has been given, the weight coefficient that is too large will lead to a slow convergence speed. It is shown that when the convergence coefficient ε is given, the weight coefficient of the BESS should satisfy

$\min \left\{w_i\right\}>\frac{n\epsilon }{2},i\in {\mathbb{Z}}_n$
to ensure a convergence of the model, and when min{w_i}=nε is satisfied, the BESS has a faster convergence speed and a better stability margin. Therefore, the weight coefficient w_i is defined as follows:

$\begin{array}{cc}{w}_i=\frac{w_i^{\prime }n\epsilon }{{\min }_{\forall i\in {\mathbb{Z}}_n}\left\{w_i^{\prime }\right\}},& \left(3\right)\end{array}$

where w′_i is expressed as

$w_i^{\prime }=\{\begin{array}{c}{E}_{\mathrm{ba},i}·{\mathrm{SOC}}_B^{\max }-E_{b,i}\left(t\right),P_{BS}^{\mathrm{ref}}\left(t\right)\le 0\\ E_{b,i}\left(t\right)-E_{\mathrm{ba},i}·{\mathrm{SOC}}_B^{\min },P_{BS}^{\mathrm{ref}}\left(t\right)>0\end{array},$
E_ba,i represents the capacity of the BESU i, E_b,i(t) represents the residual capacity at the current time of the BESU i, ε represents the convergence coefficient, n is the number of the BESUs, SOC_B^max represents an upper limit of a state of charge of each of the BESUs, and SOC_B^min represents a lower limit of the state of charge of each of the BESUs.

In an illustrated embodiment, r_ij(k) in the formula (1) is an artificially defined/introduced variable, which is the key to realize the state constraint of the consensus variable. Calculation steps of the variable r_ij(k) are described as follows.

S11, adjusting a restraint range.

A distance between the state x_i(k) and the upper limit x_i is defined as Δ_i^u(k), a distance between the state x_i(k) and the lower limit x_i is defined as Δ_i^l(k), and the distances Δ_i^u(k) and Δ_i^l(k) are respectively expressed as the following formula (4):

$\begin{array}{cc}\{\begin{array}{c}{\Delta }_i^u\left(k\right)=❘x_i\left(k\right)-{\overline{x}}_i❘\\ {\Delta }_i^l\left(k\right)=❘x_i\left(k\right)-{\underset{\_}{x}}_i❘\end{array}.& \left(4\right)\end{array}$

In order to make a state outside of the constraint range be transferred to be within the constraint range in the process of iterating, the upper and lower limits of the state during iterating are adjusted as per formulas (5) and (6) expressed as follows:

$\begin{array}{cc}{\overline{x}}_i^{\prime }=\{\begin{array}{cc}{\overline{x}}_i-{\mathrm{\tau \Delta }}_i^u\left(k\right)& x_i\left(k\right)>{\overline{x}}_i\\ {\overline{x}}_i& x_i\left(k\right)\le {\overline{x}}_i\end{array},& \left(5\right)\end{array}$ $\begin{array}{cc}{\underset{¯}{x}}_i^{\prime }=\{\begin{array}{cc}{\underset{¯}{x}}_i+{\mathrm{\tau \Delta }}_i^l\left(k\right)& x_i\left(k\right)<{\underset{¯}{x}}_i\\ {\underset{¯}{x}}_i& x_i\left(k\right)\ge {\underset{¯}{x}}_i\end{array},& \left(6\right)\end{array}$

where the formulas (5) and (6) each represent that a constraint boundary is adjusted according to a distance between a current state and the constraint boundary, so that the state outside the constraint range is transferred to be within the state constraint in the process of iterating. A value range of τ is:

$0<\tau <\min \left(\frac{❘{\overline{x}}_i❘}{{\Delta }_i^u\left(k\right)},\frac{❘{\underset{¯}{x}}_i❘}{{\Delta }_i^l\left(k\right)}\right).$

S12, calculating intermediate variables as per a formula (7) expressed as follows:

$\begin{array}{cc}\{\begin{array}{c}\Delta {\overline{x}}_i\left(k\right)=❘x_i\left(k\right)-{\overline{x}}_i^{\prime }❘\\ \Delta {\underset{\_}{x}}_i\left(k\right)=❘x_i\left(k\right)-{\underset{\_}{x}}_i^{\prime }❘\end{array}.& \left(7\right)\end{array}$

S13, calculating the variable r_ij(k).

Specifically, the intermediate variables (Δx_i(k),Δx_i(k)) and the state x_i(k) are sent to adjacent nodes by the node i, and intermediate variables and a state of adjacent j-th node are received by the i-th node for calculating the variable r_ij(k) as per a formula (8) expressed as follows:

$\begin{array}{cc}{r}_{\mathrm{ij}}\left(k\right)=\{\begin{array}{cc}\min \left\{1,\frac{\Delta {\stackrel{\_}{x}}_i\left(k\right)}{❘{\sigma }_{\mathrm{ij}}\left(k\right)❘},\frac{\Delta {\underset{\_}{x}}_j\left(k\right)}{❘{\sigma }_{\mathrm{ij}}\left(k\right)❘}\right\}& {\sigma }_{\mathrm{ij}}\left(k\right)>0\\ \min \left\{1,\frac{\Delta {\stackrel{\_}{x}}_j\left(k\right)}{❘{\sigma }_{\mathrm{ij}}\left(k\right)❘},\frac{\Delta {\underset{¯}{x}}_i\left(k\right)}{❘{\sigma }_{\mathrm{ij}}\left(k\right)❘}\right\}& {\sigma }_{\mathrm{ij}}\left(k\right)<0\\ 1& {\sigma }_{\mathrm{ij}}\left(k\right)=0\end{array},& \left(8\right)\end{array}$

where σ_ij(k)=x_j(k)−x_i(k), j ε N_i(k); x_i(k) is configured to prevent x_i(k) from exceeding the constraint boundaries thereof in the process of iterating when transferring along a direction approaching the constraint boundaries to thereby ensure that a convergence result meets the state constraint.

In an illustrated embodiment, C_ij(k)=εw_i⁻¹r_ij(k), M_ij(k)=C_ij(k)a_ij, and M_ii(k)=1−Σ_jεNi(k) M_ij(k), (∀i≠j) are defined. It can be known that ε>0, ∀i ε _n, and w_i≥nε is satisfied from preset conditions, therefore, a formula (9) is satisfied for ∀i,j ε _n, C_ij(k), M_ij(k), and M_ii(k) as follows:

$\begin{array}{cc}\{\begin{array}{c}{C}_{\mathrm{ij}}\left(k\right)\ge 0,{\sum }_{j\in N_i\left(k\right)}{C}_{\mathrm{ij}}\left(k\right)<1\\ M_{\mathrm{ij}}\left(k\right)\ge 0,M_{\mathrm{ii}}\left(k\right)>0,{\sum }_{j=1}^n{M}_{\mathrm{ij}}\left(k\right)=1\end{array}.& \left(9\right)\end{array}$

Therefore, the formula (1) can be expressed equivalently as:
x_i(k+1)=M_ii(k)x_i(k)+Σ_jεNi(k)M_ij(k)x_j(k)  (10).

A consensus control rule of the BESS is expressed in a matrix form as follows:
X(k+1)=M(k)X(k)  (11).

From the formula (9), X(k)=[x₁(k),x₂(k), . . . ,x_n(k)]^T. From the formula (9), M(k) in the formula (11) is a non-negative row random matrix, its diagonal element M_ii(k) is a state conversion term, and its non-diagonal element M_ij(k) is a dynamic coupling term. The j-th node has a direct impact on the i-th node when M_ijk)≠0 is satisfied, otherwise, there is no direct impact.

S2, improving the convergence speed based on the DMPC.

Improving consensus iterative speed of the distributed coordination power allocation is conducive to reducing the response delay of the BESS and improving the effect of energy storage regulation. The introduction of the DMPC into the consensus algorithm can significantly improve the convergence speed. In this embodiment, the formula (11) is further improved based on the DMPC, a decoupling control quantity U(k) is first introduced into on a basis of the formula (11) expressed as follows:
X(k+1)=M(k)X(k)+W⁻¹U(k)  (12),

where W=w diag(w₁, w₂, . . . , w_n) represents a weight matrix of the BESS; X₀=[w₁⁰,x₂⁰, . . . ,x_n⁰]^T, where x_n⁰ represents a state of a n-th BESU at an initial time, X(k)=[x₁(k),x₂(k), . . . ,x_n(k)]^T, where x_n(k) is a state of the n-th BESU at a time k; and X (k+1)=[x₁(k+1),x₂(k+1), . . . ,x_n(k+1)]^T, where x_n(k+1) is a state of an n-th BESU at a time k+1. The current time is set to be the time k, a prediction time domain is defined as N_p, a control time domain is defined as N_c, and N_p≥N_c, and the decoupling control quantity U(k) outside the control time domain N_c remain unchanged. The control variable of the i-th node only acts on itself, and an updating formula of the i-th node is expressed as follows:
x_i(k+1)=M_ii(k)x_i(k)+Σ_j∈Ni(k)M_ij(k)x_j(k)+w_i⁻¹u_i(k)  (12a).

Taking the current time k as a starting point, the state vector X_i(k) of the i-th node can be predicted in the prediction time domain based on the formula (12a) as a formula (13):
X_i(k)=F_i(k)x_i(k)+Σ_j∈Ni(k)F_ij(k)X_j(k)+G_i(k)U_i(k)  (13),

where

$X_i\left(k\right)=\left[\begin{array}{c}{x}_i\left(k+1/k\right)\\ x_i\left(k+2/k\right)\\ ⋮\\ x_i\left(k+N_p/k\right)\end{array}\right],$ $F_{\mathrm{ij}}\left(k\right)=\left[\begin{array}{cccc}{M}_{ij}\left(k\right)& 0& \dots & 0\\ M_{\mathrm{ii}}\left(k\right)M_{\mathrm{ij}}\left(k\right)& M_{ij}\left(k\right)& \dots & 0\\ ⋮& \ddots & \dots & ⋮\\ M_{\mathrm{ii}}^{N_p-1}\left(k\right)M_{\mathrm{ij}}\left(k\right)& M_{\mathrm{ii}}^{N_p-2}\left(k\right)M_{\mathrm{ij}}\left(k\right)& \dots & M_{ij}\left(k\right)\end{array}\right],$ $F_i\left(k\right)=\left[\begin{array}{c}{M}_{\mathrm{ii}}\left(k\right)\\ M_{\mathrm{ii}}^2\left(k\right)\\ ⋮\\ M_{\mathrm{ii}}^{N_p}\left(k\right)\end{array}\right],$ $G_i\left(k\right)=\left[\begin{array}{cccc}{w}_i^{-1}& & & \\ M_{\mathrm{ii}}\left(k\right)w_i^{-1}& w_i^{-1}& & \\ ⋮& \dots & \ddots & \\ M_{ii}^{N_P-2}\left(k\right)w_i^{-1}& M_{ii}^{N_P-2}\left(k\right)w_i^{-1}& \dots & {\sum }_{c=0}^{N_P-N_c}{M}_{ii}^c\left(k\right)w_i^{-1}\end{array}\right],$ $\mathrm{and}{X}_j\left(k\right)=x_j\left(k\right)I_{N_p},$
where I_Np represents a column vector whose elements are all 1.

A quadratic optimization objective function is defined as a formula (14):
minJ_i(k)=∥X_i(k)−H_i(k)∥_Q²+∥U_i(k)∥_R²  (14)

where Q=qI_Np, R=rI_Nc, q>0, r>0; H_i(k) represents a reference vector of the consensus variable, H_i(k) and U_i(k) are calculated as follows.

1) When x_i(k) ε [x_i,x_i] is satisfied,

$\begin{array}{cc}{H}_i\left(k\right)\mathrm{is}\mathrm{expressed}\mathrm{as}{H}_i\left(k\right)=\{\begin{array}{cc}{\overline{x}}_i^{\prime }{I}_{N_p}& x_i\left(k\right)>{\overline{x}}_i\\ {\underset{\_}{x}}_i^{\prime }{I}_{N_p}& x_i\left(k\right)<{\underset{\_}{x}}_i\end{array};& \left(15\right)\end{array}$

for the formula (14), a partial derivative of U_i(k) is calculated and the partial derivative is made to be 0, so as to obtain an optimal decoupling control quantity, which is expressed as a formula (16):

$\begin{array}{cc}\{\begin{array}{c}{U}_i\left(k\right)={\Psi }_i\left(k\right)\left[H_i\left(k\right)-F_i\left(k\right)-{\sum }_{j\in N_i\left(k\right)}{F}_{\mathrm{ij}}\left(k\right)X_j\left(k\right)\right]\\ {\Psi }_i\left(k\right)={\left(G_i^T\left(k\right)Q{G}_i\left(k\right)+R\right)}^{-1}\left(G_i^T\left(k\right)Q\right)\end{array}.& \left(16\right)\end{array}$

2) When x₁(k) ε [x_i,x_i] is satisfied,
U_i(k) is expressed as U_i(k)=0·I_Np  (17),

so far, the calculation of U_i(k) is completed, and then a first element u_i(k) in U_i(k) is taken and corrected for a second time. A correction formula is shown in a formula (18) expressed as:
u′_i(k)=Σ_j∈Nia_ij(u_i(k)−u_j(k))  (18).

The corrected control vector, i.e., the decoupling control quantity, U_i(k)=[u′₁(k),u′₂(k), . . . , u′_n(k)]^T, is input into the formula (12) to complete the updating of the BESS. It can be seen that the above derivation process only needs to transfer a small amount of information between adjacent nodes, which is a completely distributed process.

S3, implementing of the distributed algorithm

To realize the grouped distributed coordination control of the BESS using the distributed algorithm, iterative initialization and the adjustment strategy of the adjacency matrix A are needed. Iterative initialization at the initial time of each control cycle is a necessary link to realize the distributed coordination control of the BESS, and the adjustment strategy of the adjacency matrix A is designed to realize the grouping coordinated control and improve the energy conversion efficiency of the BESU.

Step S31, initializing the BESS.

Specifically, after receiving a total load, i.e., a total output power, the distributed BESS needs to be initialized for the consensus iteration according to information of the received total output power, and then the consensus iteration can be carried out. P_bi⁰ is defined as an initial power of the iterating in each control cycle of the BESU i and x_i⁰ is an initial value of the consensus variable of the iterating in each control cycle of the BESU i. In order to reduce the communication cost between a superior energy management center and the BESS, it is assumed that the superior energy management center randomly sends the total output power P_BS^ref(t) of the BESS to one BESU of the charging group C_Z and one BESU of the discharging group D_Z at a beginning of each control cycle.

1) defining an initial power P_bi⁰ for iterating.

Specifically, the total output power of the BESS is defined as P_BS^ref(t), one of the BESUs of the charging group C_Z receives the total output power P_BS^ref(t) when P_BS^ref(t)≤0, and one of the BESUs of the discharging group D_Z receives the total output power P_BS^ref(t) when P_BS^ref(t)>0. Therefore, for the BESU i received the total output power P_BS^ref(t), the initial power P_bi⁰ is expressed as follows:

$\begin{array}{cc}{P}_{bi}^0=\{\begin{array}{cc}{P}_{BS}^{\mathrm{ref}}\left(t\right),& i\in C_Z\mathrm{and}{P}_{BS}^{\mathrm{ref}}\left(t\right)\le 0\\ P_{BS}^{\mathrm{ref}}\left(t\right),& i\in D_Z\mathrm{and}{P}_{BS}^{ref}\left(t\right)>0\\ 0,& \mathrm{other}\mathrm{conditions}\end{array},& \left(19\right)\end{array}$

where the BESU i belongs only to the charging group C_Z or the discharging group D_Z at the time t, and the initial power P_bi⁰ is 0 for the other BESUs that does not receive the total output power P_BS^ref(t).

2) defining the initial state quantity x_i⁰ for iterating, as per a formula (20) expressed as follows:
x_i⁰=w_i⁻¹P_bi⁰  (20).

After the BESU i determines its own value of the initial power P_bi⁰, the initial state quantity for iterating is calculated according to an inverse transformation of the formula (1a) as follows:
x_i⁰=w_i⁻¹P_bi⁰  (20a).

Therefore, the initialization of the distributed BESS is expressed as:

$\begin{array}{cc}\{\begin{array}{c}{P}_0\left(t\right)={\left[P_{b1}^0,P_{b2}^0,\dots ,P_{bn}^0\right]}^T\\ X_0=W^{-1}{P}_0\left(t\right)\end{array}.& \left(21\right)\end{array}$

Step S32, determining the adjustment strategy of the adjacency matrix A.

Specifically, a condition that the charging group C_Z and the discharging group D_Z of the BESS at an initial time of each control cycle are set to do not exchange information is set, that is, the BESU i and the dual unit BESU i⁰ thereof correspond to a_ii∘=0, while the charging group C_Z and the discharging group D_Z each are intragroup communicated as C_Z(D_Z), and a corresponding adjacency matrix is an initial adjacency matrix A₀ expressed as follows:

$A_0=\left[\begin{array}{cccccccc}0& 1& 0& 1& 0& 0& 0& 0\\ 1& 0& 1& 0& 0& 0& 0& 0\\ 0& 1& 0& 1& 0& 0& 0& 0\\ 1& 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& 0& 1\\ 0& 0& 0& 0& 1& 0& 1& 0\\ 0& 0& 0& 0& 0& 1& 0& 1\\ 0& 0& 0& 0& 1& 0& 1& 0\end{array}\right].$

The adjacency matrix A of the BESS at the initial time of iterating of each control cycle is expressed as A=A₀, and the adjacency matrix A changes with iteration times k. The following adjustment strategy of the adjacency matrix A is designed to solve intergroup coordination and avoid the low energy conversion efficiency caused by the low operating power of the BESUs.

1) Responding to the Intergroup Coordination.

Specifically, k_c is defined as a positive integer and X₀ is substituted into the formula (12) for iterating. When a k_c-th time of iterating is ended, the BESU i can determine the total output power P_BS^ref(t) of the BESS at this time exceeds a total response capacity of the group where the BESU i is located in response to the state x_i(k_c) of the BESU i still does not meet the state constraint, and the dual group is requested to participate in consumption and let a_ii∘=1. Otherwise let a_ii∘=0 and continue to iterating until reaching the convergence accuracy. The adjacency matrix A is adjusted according to a formula (22) to achieve the intergroup coordination:

$\begin{array}{cc}\{\begin{array}{cc}{a}_{{\mathrm{ii}}^o}=1& x_i\left(k_c\right)\notin \left[{\underset{¯}{x}}_i,{\overline{x}}_i\right]\\ a_{{\mathrm{ii}}^o}=0& x_i\left(k_c\right)\in \left[{\underset{\_}{x}}_i,{\overline{x}}_i\right]\end{array},& \left(22\right)\end{array}$

where x_i and x_i respectively represent a lower limit and an upper limit of a state constraint of the BESU i as an i-th node, x_i(k_c) represents a state of a consensus variable of the i-th node at a time k_c, a_ii∘ represents an element in the adjacency matrix A, a_ii∘=1 represents that the BESU i exchanges data with the BESU i^∘ being a dual unit of the BESU i, and a_ii∘=0 represents there is no data exchange.

2) Improving the Energy Conversion Efficiency of the BESUs.

For the BESU i with an initial power of P_bi⁰=P_BS^ref(t), the BESU i does not exchange information with its adjacent nodes when a formula (23) is satisfied, and only the BESU i bears by the BESU i, that is, P_bi(t)=P_BS^ref(t). Otherwise, continue iterating. The formula (23) is expressed as follows:
|P_Bi⁰|≤2ζP_bn,i  (23),

where ζ ε (0,1), P_bn,i is a rated power of the BESU i.

For ∀i ε _n, when the k-th iteration result has reached a certain convergence accuracy and the power P_bi^k corresponding to x_i(k) (see the formula (1a)) satisfies the formula (24), then the BESU i exits the communication iteration of a current control cycle. In this situation, the BESS needs to be adjusted as shown in a formula (25). Otherwise, no adjustment will be made and the iteration will continue until reaching the convergence accuracy. The formula (24) and the formula (25) are respectively expressed as follows:

$\begin{array}{cc}0<❘P_{\mathrm{bi}}^k❘\le \zeta P_{\mathrm{bn},i},& \left(24\right)\end{array}$ $\begin{array}{cc}\{\begin{array}{c}{x}_j^{\prime }\left(k\right)=x_j\left(k\right)+\frac{x_i\left(k\right)w_i}{d_i\left(k-1\right)w_j}\\ x_i^{\prime }\left(k\right)=0\\ a_{\mathrm{ij}}=0\\ d_j^{\prime }\left(k\right)=d_j\left(k\right)-1\end{array}\forall j\in N_i\left(k-1\right),& \left(25\right)\end{array}$

where a_ij represents an element in the adjacency matrix A representing a communication network topology, x_i(k) represents a state of the i-th node at the time k, x_j(k) represents a state of a j-th node at the time k, w_i represents a weight coefficient of the BESU i, w_j represents a weight coefficient of a BESU j, N_i(k−1) represents an adjacent node set of the BESU i before exits, x′_j(k) represents a state quantity of the j-th node in N_i(k−1) after the BESU i exits, and d′_j(k) represents the number of adjacent nodes of the j-th node in N_i(k−1) after the BESU i exits. As can be seen from the formula (25), before and after the i-th node exits, there are:
Σ_{j∈Ni(k−1)∪{i}}w_jx′_j(k)=Σ_{jÅNi(k−1)∪{i}}w_jx_j(k)  (26).

Therefore, the consensus algorithm proposed by the disclosure can ensure that the k-th iteration result of the BESS satisfies a formula (27), so as to ensure that a sum of final power allocations of the multiple BESUs is equal to the total output power P_BS^ref(t) of the BESS. The formula (27) is expressed as follows:
w_ix_i(0)=w_ix_i(k)  (27).

To sum up, at the beginning of each control cycle, each BESU in the BESS is initialized according to the above method, the initial state value X₀ is substituted into the formula (12), the convergence accuracy is set, and a convergence state value X* that satisfies the power constraints of all the n number of nodes can be obtained according to the above method and strategy. Then, the power P_bi(t) of each BESU is calculated by using a formula (28) as follows:

$\begin{array}{cc}\{\begin{array}{c}P\left(t\right)={\mathrm{WX}}^{*}\\ P\left(t\right)={\left[P_{b1}\left(t\right),P_{b2}\left(t\right),\dots ,P_{\mathrm{bn}}\left(t\right)\right]}^T\\ P_B\left(t\right)={\sum }_{i=1}^n{P}_{\mathrm{bi}}\left(t\right)\end{array},& \left(28\right)\end{array}$

where X* represents the convergence state value that satisfies the power constraints of all the n number of nodes, W=w diag(w₁, w₂, . . . , w_n) represents the weight matrix of the BESS, P_B(t) represents a sum of the powers of all the n number of BESUs, and the power of the BESU i.

The power allocation of each BESU in the BESS in this control cycle is completed through the distributed algorithm. When the SOC of any group of BESUs reaches the set SOC upper or lower limits for switching charging/discharging states during operation, the state switching of the charge/discharging group is triggered, that is, the original charging group C_Z is changed into the discharging group D_Z, and the original discharging group D_Z is changed in the charging group C_Z.

Verification of the Feasibility of the Method

Specifically, the advantages of the proposed method are demonstrated from two aspects including algorithm performance and control effect. The parameters of the BESS and the algorithm are shown in Table 1 as follows.

TABLE 1
Energy storage system parameters	Algorithm parameters
	Numerical		Numerical
Parameters	value	Parameters	value
BESU quantity	n = 8	Convergence	ε = 0.1
		coefficient
BESU rated power/(kW)	Pbn = 10	Convergence	10−4
		accuracy
BESU rated capacity/	Eba = 40	Prediction time	Np = 5
(kW · h)		domain Np
Upper and lower limits	0.9/0.1	Control time domain	Nc = 3
of SOC SOCminmax		Nc
Upper and lower thresholds	0.88/0.122	ζ	ζ = 0.2
of SOC of charging/
discharging states switching

1) Algorithm Performance

The experiment may include two parts. Firstly, an energy storage system with four energy storage units is taken as an example, the classical weighted-consensus algorithm and the algorithm proposed in this paper are compared and analyzed. Secondly, another energy storage system with 8 energy storage units is taken as an example, the adjustment strategy of the adjacency matrix proposed in the S32 is analyzed to verify the advantages of the algorithm proposed in this paper.

The energy storage system with four energy storage units (where the communication topology is ring-shaped) as an example, a discharge weight matrix of the time t is defined as W=diag(1, 2, 3, 4), and the total output power P_BS^ref(t)=30 kW is sent to No. 3 energy storage unit by the superior energy management center. The iterative initialization of the classical weighted-consensus algorithm and the algorithm in this paper are P₀=[0,0,30,0]^T and X₀=[0,0,10,0]^T. FIG. 2a illustrates an iterative convergence process of the classical weighted-consensus algorithm, and the convergence result is P(t)=[3,6,9,12]^T. FIG. 2b illustrates an iterative process of the algorithm proposed in this paper, and the power allocation result is: P(t)=[3.3333, 6.6667, 10, 10]^T.

The adjustment strategy of the adjacency matrix A is analyzed by using the battery energy storage unit system with eight energy storage units. The eight energy storage units are divided into two groups Z1 and Z2, the group Z1 is composed of energy storage units from No. 1 through No. 4, the group Z2 is composed of energy storage units from No. 5 through No. 8, and the group Z1 is the charging group and the group Z2 is the discharging group. The communication topology of the energy storage system is shown in FIG. 10, where ζ=0.2. The discharge weight matrix of the energy storage system at the time t as: W=diag(1, 1.5, 2, 4, 5, 5.5, 6, 7), and the P_BS^ref(t)=50 kW is sent to No. 5 energy storage unit in response to a demand by the microgrid center. Iteration initialization is expressed as: P₀=[0,0,0,0,50,0,0,0]^T and X₀=[0,0,0,0,10,0,0,0]^T. The iterative convergence process is shown in FIG. 3.

2) Control effect of BESS grouping strategy

Specifically, the control effect of this method is analyzed and compared with the control effect of traditional non-grouping strategies (mainly including maximum power allocation strategy and proportional allocation strategy) through this experiment. A grid-connected microgrid with wind charge storage and solar charge storage is built as shown in FIG. 4. In the microgrid model, the parameters of the distributed energy storage system are shown in the Table 1, and the communication topology is shown in FIG. 10. In the energy storage system, The BESUs from No. 1 through No. 4 form the group Z1, and the BESUs from No. 5 through No. 8 form the group Z2. The source load data comes from an experimental microgrid platform, and a sampling period is 5 minutes. In this scenario, due to the large vitality of new energy and load power, in order to reduce the adverse impact of microgrid on large power grid, it is necessary to use energy storage to smooth the grid-connected tie-line power. The group Z1 is set as the charging group and the group Z2 as the discharging group at the initial time, and the initial SOC of the BESUs from No. 1 through No. 8 is (0.2, 0.25, 0.3, 0.35, 0.65, 0.7, 0.75, 0.8). A curve of the total power P_BS^ref(t) of the energy storage system in a scheduling period of the microgrid to cycle is shown in FIG. 5.

Under the grouping control strategy in this paper, power curves of the group Z1 and the group Z2 are shown in FIG. 6. In order to demonstrate the smoothing effect of the grouping control strategy in this paper on the grid-connected power fluctuation, the restriction condition of the microgrid grid-connected power fluctuation rate γ is γ≤7% within 30 minutes. Under this restriction condition, the strategy proposed in the disclosure, the traditional non-grouping strategy and traditional DBESS (double battery energy storage system) strategy are statistically compared in one scheduling period, and statistical results are shown in FIG. 7. In addition, SOC change curves of respective BESUs within 24 hours under the strategy proposed in the disclosure is shown in FIG. 8.

In order to demonstrate the effectiveness of the proposed strategy in reducing the switching times of charging/discharging states of BESUs, the change of power symbols of BESUs at adjacent moments is defined as one switching of charging/discharging states, and the switching times of charging/discharging states of each BESU in one scheduling period under the proposed strategy, the maximum power allocation strategy in the non-grouping strategy and the proportional allocation strategy are counted respectively. The statistical results are shown in FIG. 9. In addition, in order to verify the effectiveness of the proposed strategy in this paper to improve the efficiency, the number of low-power operations of all BESUs under different strategies is counted here. The operating state of the BESU with operating power 0<|P_b(t)|≤0.2 P_bn is defined as a low-power and low-efficiency operation with reference to relevant research. The times of BESU operating power less than 0.2 P_bn under different strategies are compared, and the proportion of low-power operation is calculated. The statistical results of the low-power operation times per 24 hours of the BESU under different allocation strategies are shown in Table 2:

	TABLE 2
	Low-power operation	Total operation
	times	times	Low-power
	Total times of	Average	Total times of	Average	operation
Power allocation strategy	eight BESUs	times	eight BESUs	times	ratio
Proposed strategy	23	2.875	1301	162.625	1.77%
Non-	Maximum power	64	8	983	122.875	6.51%
strategy	allocation
grouping	Proportional	779	97.375	2304	288	33.81%
	allocation

Through the above two experiments, this study draws the following conclusions:

1) Firstly, because the classical consensus algorithm does not consider the constraints, the phenomenon that the power of the energy storage unit with larger weight exceeds the limit occurs, while the consensus algorithm proposed in this paper considers the power constraints and improves the convergence speed of the algorithm based on DMPC. Therefore, the power allocation results of the algorithm in this paper meet the constraints and the principle of “capable people do more work”, and the convergence speed of the algorithm in this paper is faster than that of the classical consensus algorithm. In addition, the adjacency matrix adjustment strategy proposed in this paper can realize the intergroup coordination of the energy storage units. When the total regulation demand is greater than the regulation capacity of a single group, the unit of the dual group also participate in the regulation, so as to give full play to the regulation capacity of the energy storage system, overcome the shortcomings of insufficient power capacity of the energy storage system under the traditional DBESS grouping control strategy, and avoid multiple energy storage units running at low-power conditions at the same time, to thereby improve the energy conversion efficiency of the BESUs.

2) FIG. 6 illustrates power curves of the group Z1 and the group Z2. At about 13:00, the states of the two groups of batteries switch, that is, the Z1 group changes from a charging state to a discharging state, and the Z2 group changes from a discharging state to a charging state. For the color overlapping part in FIG. 6, due to the violent fluctuation of new energy power in these periods, the response demand for BESS is greater than the rated power of a single group, and both groups of batteries will participate in the response. From the volatility shown in FIG. 7, it can be seen that the smoothing results of the grouping strategy in this paper meet the volatility constraints and can effectively smooth the grid-connected power fluctuation. FIG. 8 illustrates SOC change curves of respective energy storage units within 24 hours. Based on comprehensive analysis of FIG. 6, FIG. 7 and FIG. 8, it can be seen that the grouping strategy in this paper has advantages as follows. Firstly, the capacity utilization of respective BESUs is improved. Secondly, the BESUs in the same group can be in the same charging/discharging state for a long time, so that the orderliness of the charging/discharging behavior of each BESU is improved, and the frequent conversions of the charging/discharging state in the operation process of the BESUs is avoided. Thirdly, the SOC states of the BESUs in the same group tend to be consistent by adopting consensus control, so that individual BESUs are prevented from being unable to participate in response when entering extreme conditions, and the power response capability of the energy storage system is ensured.

The statistical results in FIG. 9 show that the average charging/discharging conversion times of energy storage units in the proposed strategy, the maximum power allocation strategy and the proportional allocation strategy are 51, 76 and 141 times respectively in a scheduling period. The switching times of charging/discharging states of the proposed strategy are 67% of the maximum power allocation strategy and 36% of the proportional allocation strategy. Moreover, the switching times of charging/discharging states of each BESU in the proposed strategy are relatively uniform, which is helpful to prevent the rapid decline of the life of individual units and to prolong the cycle life of BESS. The experimental results shown in the Table 2 show that the BESU operates at low power less than 3 times on average in 288 control cycles a day under the proposed strategy, and the low power operation rate is 4.74% less than that of the maximum power allocation strategy and 32.04% less than that of the proportional allocation strategy. The control strategy proposed in this paper optimizes the operating state of BESU, reduces the power loss of BESU, and improves the energy conversion efficiency of BESU.

Through this research, this paper finds that the following effects can be achieved by grouping and distributed consensus control of multiple energy storage units. Firstly, the switching times of charging/discharging states of energy storage units are reduced and the service life of energy storage is prolonged. Secondly, the times of low-power operation of energy storage unit are reduced and the efficiency of energy storage operation is improved. Thirdly, the consensus of each unit in the energy storage system is improved and the regulation capability of the energy storage system is improved. Fourthly, the distributed control of energy storage system is realized and the advantages of distributed control are fully played. It provides a reference idea for the application of energy storage technology. In addition, the proposed weighted-consensus algorithm based on DMPC and consensus variable constraints effectively improves the convergence speed of the algorithm and ensures that the power allocation results meet the power constraints of each energy storage unit.

The above mentioned is only preferred specific embodiments of the disclosure, but the scope of protection of the disclosure is not limited to this. Any changes or substitutions that can be easily thought of by those skilled in the related art within the technical scope disclosed in the disclosure shall be covered by the scope of protection of the disclosure. Therefore, the scope of protection of the disclosure should be subject to the scope of protection of the claims.

Claims

1. A grouped consensus power allocation method for multiple energy storage units, comprising: { a ii o = 1 x i ( k c ) ∉ [ x _ i, x _ i ] a ii o = 0 x i ( k c ) ∈ [ x _ i, x _ i ], ( 22 ) { x j ′ ( k ) = x j ( k ) + x i ( k ) ⁢ w i d i ( k - 1 ) ⁢ w j x i ′ ( k ) = 0 a ij = 0 d j ′ ( k ) = d j ( k ) - 1 ⁢ ∀ j ∈ N i ( k - 1 ) ( 25 ) { P ⁡ ( t ) = WX * P ⁡ ( t ) = [ P b ⁢ 1 ( t ), P b ⁢ 2 ( t ), L, P bi ( t ) ⁢ L, P bn ( t ) ] T P B ( t ) = ∑ i = 1 n ⁢ P bi ( t ) ( 28 )

determining a grouped coordination control strategy of multiple battery energy storage units (BESUs), wherein a battery energy storage system (BESS) comprises the multiple BESUs, the multiple BESUs are communicated through a communication topology, and the multiple BESUs are grouped into a charging group CZ and a discharging group DZ; and

performing power coordination control between the charging group CZ and the discharging group DZ based on a distributed algorithm, and determining whether a switching between the charging group CZ and the discharging group DZ is triggered during an operation process of the BESS;

wherein the performing power coordination control between the charging group CZ and the discharging group DZ based on a distributed algorithm comprises:

step S31, initializing the BESS; and

step S32, determining an adjustment strategy of an adjacency matrix A, comprising: determining, by a BESU i, whether a dual group is needed to participate in a process of iterating, and adjusting the adjacency matrix A through a formula (22) to achieve intergroup coordination in response to the dual group is needed to participate, wherein the formula (22) is expressed as follows:

where xi and xi respectively represent a lower limit and an upper limit of a state constraint of the BESU i as an i-th node, xi(kc) represents a state of a consensus variable of the i-th node at a time kc, aii∘ represents an element in the adjacency matrix A, aii∘=1 represents that the BESU i exchanges data with the BESU i∘ being a dual unit of the BESU i, otherwise, there is no data exchange; determining whether the BESU i with an initial power of Pbi0=PBSref(t) exits a communication iteration of a current control cycle, and adjusting the BESS as per a formula (25) in response to it is determined that the BESU i with the initial power exits the communication iteration, otherwise continuing to iterating without the adjusting until reaching a convergence accuracy, wherein the formula (25) is expressed as follows:

where aij represents an element in the adjacency matrix A representing a communication network topology, xi(k) represents a state of the i-th node at a time k, xj(k) represents a state of a j-th node at the time k, wi represents a weight coefficient of the BESU i, wj represents a weight coefficient of a BESU j, Ni(k−1) represents an adjacent node set of the BESU i before exits, xj′(k) represents a state quantity of the j-th node in Ni(k−1) after the BESU i exits, and dj′(k) represents the number of adjacent nodes of the j-th node in Ni(k−1) after the BESU i exits; and calculating a power Pbi(t) of each of the multiple BESUs as per a formula (28) in response to the convergence accuracy is reached, wherein the formula (28) is expressed as follows:

where X* represents a convergence state value meeting power constraints of all the n number of nodes, W=diag(w1, w2,..., wn) represents a weight matrix of the BESS, PB(t) represents a sum of powers of all the multiple BESUs, and Pbi(t) represents the power of the BESU i.

2. The method according to claim 1, wherein the determining whether a switching between the charging group CZ and the discharging group DZ is triggered comprises:

switching charging/discharging states of the charging group CZ and the discharging group DZ, in response to a state of charge (SOC) of any one group of the BESUs reaches an upper threshold or a lower threshold.

3. The method according to claim 2, wherein the charging group CZ and the discharging group DZ are defined as mutual dual groups;

wherein a BESU i of the charging group CZ and a BESU i∘ of the discharging group DZ are mutual dual units when the BESU i and the BESU i∘ meet a direct intergroup physical communication condition, where i∈{1, 2..., n}, i∘∈{1, 2..., n}, and n is the number of the BESUs in the charging group CZ.

4. The method according to claim 1, wherein the step S31 comprises: P bi 0 = { P BS ref ( t ), i ∈ C Z ⁢ and ⁢ P BS ref ( t ) ≤ 0 P BS ref ( t ), i ∈ D Z ⁢ and ⁢ P BS ref ( t ) > 0 0, other ⁢ conditions, ( 19 ) { P 0 ( t ) = [ P b ⁢ 1 0, P b ⁢ 2 0, …, P bn 0 ] T X 0 = W - 1 ⁢ P 0 ( t ). ( 21 )

1) defining a total output power of the BESS as PBSref(t) and an initial power for iterating of the BESU i as Pbi0; wherein one of the BESUs of the charging group CZ receives the total output power PBSref(t) when PBSref(t)≤0, and one of the BESUs of the discharging group DZ receives the total output power PBSref(t) when PBSref(t)>0, thereby, for the BESU i received the total output power PBSref(t), the initial power Pbi0 is expressed as follows:

where the BESU i belongs only to the charging group CZ or the discharging group DZ at a time t, and the initial power Pbi0 is 0 for the other BESUs that do not receive the total output power PBSref(t);

2) Defining an initial state quantity for iterating, as per a formula (20) expressed as follows: xi0=wi−1Pbi0  (20),

where Pbi0 represents the initial power for iterating of each control cycle of the BESU i, and xi0 represents an initial value of the consensus variable of each control cycle of the BESU i;

wherein the initializing the BESS is expressed as follows:

5. The method according to claim 1, wherein before the step S32, the method further comprises: A 0 = [ 0 1 0 1 0 0 0 0 1 0 1 0 0 0 0 0 0 1 0 1 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 1 0 1 0 0 0 0 0 0 1 0 1 0 0 0 0 1 0 1 0 ],

obtaining the adjacency matrix A, comprising: making the charging group CZ and the discharging group DZ of the BESS at an initial time of each control cycle do not exchange information, that is, the BESU i and the dual unit BESU i∘ thereof correspond to aii∘=0, while the charging group CZ and the discharging group DZ each are intragroup communicated, and a corresponding adjacency matrix is an initial adjacency matrix A0 expressed as follows:

expressing the adjacency matrix A of the BESS at the initial time of iterating of each control cycle as A=A0, wherein the adjacency matrix A changes with iteration times k.

6. The method according to claim 1, wherein the determining, by a BESU i, whether a dual group is needed to participate in a process of iterating comprises:

defining the kc as a positive integer, substituting X0 into a formula (12) for iterating: X(k+1)=M(k)X(k)+W−1U(k)  (12)

where M(k) represents a non-negative row random matrix, U(k) represents a decoupling control quantity; wherein X0=[x10, x20,..., xn0]T, where xn0 is a state of an n-th BESU at the initial time; X(k)=[x1(k), x2(k),..., xn(k)]T, where xn(k) is a state of the n-th BESU at the time k; and X (k+1)=[x1(k+1), x2(k+1),..., xn(k+1)]T, where xn(k+1) is a state of the n-th BESU at a time k+1;

when a kc-th time of iterating is ended, determining, by a BESU i, the total output power PBSref(t) of the BESS at this time exceeds a total response capacity of the group where the BESU i is located in response to the state xi(kc) of the BESU i still does not meet the state constraint, requesting the dual group to participate in consumption and letting aii∘=1; otherwise letting aii∘=0 and continuing to iterating until reaching the convergence accuracy.

7. The method according to claim 6, further comprising: x i ( k + 1 ) = x i ( k ) + ε ⁢ w i - 1 ⁢ ∑ j ∈ N i ( k ) r ij ( k ) ⁢ a ij · ( x j ( k ) - x i ( k ) ), ∀ i ∈ n, ( 1 )

setting the state constraint, comprising: setting a state of the consensus variable of the i-th node at the time k as xi(k) in considering of the BESS with n numbers of agent nodes; wherein the state constraint of the consensus variable of the i-th node is expressed as follows: xi(k)∈[xi,xi],

where [x, x] represents an intersection of constraints of all the n number of nodes, [x,x]=Ii=1n[xi,xi], and xi represents a power related quantity;

introducing a variable rij(k) on a basis of an ordinary weighted-consensus algorithm to realize that a distributed power coordinated allocation result meets power constraints of the respective BESUs, wherein an iterating and updating method of the consensus variable is expressed by a formula (1) as follows:

where ε represents a convergence coefficient and ε>0; for an undirected graph, aij=1 when the i-th node and the j-th node have direct data exchange, otherwise aij=0; Ni represents an adjacent node set of the i-th node, depending on the communication network topology and changing with a change of the communication network topology; cn represents a set of the all the n number of nodes of the BESS; and rij(k) represents an artificially introduced variable.

8. The method according to claim 7, further comprising: w i = w i ′ ⁢ n ⁢ ε min ∀ i ∈ ℤ n { w i ′ } ( 3 ) w i ′ = { E ba, i · SOC B max - E b, i ( t ), P B ( t ) ≤ 0 E b, i ( t ) - E ba, i · SOC B min, P B ( t ) > 0, Eba,i represents a capacity of the BESU i, Eb,i(t) represents a residual capacity at current time of the BESU i, n is the number of the BESUs, SOCBmax represents an upper limit of a state of charge of each of the BESUs, and SOCBmin represents a lower limit of the state of charge of each of the BESUs.

obtaining the weight coefficient by a formula (3) expressed as follows:

where wi′ is expressed as

9. The method according to claim 7, further comprising: { Δ i u ( k ) = ❘ "\[LeftBracketingBar]" x i ( k ) - x _ i ❘ "\[RightBracketingBar]" Δ i l ( k ) = ❘ "\[LeftBracketingBar]" x i ( k ) - x _ i ❘ "\[RightBracketingBar]", ( 4 ) x _ i ′ = { x _ i - τ ⁢ Δ i u ( k ) x i ( k ) > x _ i x _ i x i ( k ) ≤ x _ i, ( 5 ) x _ i ′ = { x _ i + τ ⁢ Δ i l ( k ) x i ( k ) < x _ i x _ i x i ( k ) ≥ x _ i, ( 6 ) 0 < τ < min ⁡ ( ❘ "\[LeftBracketingBar]" x _ i ❘ "\[RightBracketingBar]" Δ i u ( k ), ❘ "\[LeftBracketingBar]" x _ i ❘ "\[RightBracketingBar]" Δ i l ( k ) ); { Δ ⁢ x _ i ( k ) = ❘ "\[LeftBracketingBar]" x i ( k ) - x _ i ′ ❘ "\[RightBracketingBar]" Δ ⁢ x _ i ( k ) = ❘ "\[LeftBracketingBar]" x i ( k ) - x _ i ′ ❘ "\[RightBracketingBar]"; ( 7 ) r ij ( k ) = { min ⁢ { 1, Δ ⁢ x _ i ( k ) ❘ "\[LeftBracketingBar]" σ ij ( k ) ❘ "\[RightBracketingBar]", Δ ⁢ x _ j ( k ) ❘ "\[LeftBracketingBar]" σ ij ( k ) ❘ "\[RightBracketingBar]" } σ ij ( k ) > 0 min ⁢ { 1, Δ ⁢ x _ j ( k ) ❘ "\[LeftBracketingBar]" σ ij ( k ) ❘ "\[RightBracketingBar]", Δ ⁢ x _ i ( k ) ❘ "\[LeftBracketingBar]" σ ij ( k ) ❘ "\[RightBracketingBar]" } σ ij ( k ) < 0 1 σ ij ( k ) = 0 ( 8 )

obtaining the introduced variable rij(k), comprising: adjusting a constraint range by: defining Δiu(k) as a distance between the state xi(k) and the upper limit xi and defining Δil (k) as a distance between the state xi(k) and the lower limit xi, wherein Δiu(k) and Δil(k) are respectively expressed as follows:

adjusting the upper limit and the lower limit during iterating as per a formula (5) and a formula (6) in order to make a state outside of the constraint range be transferred to be within the constraint range in a process of iterating, wherein the formula (5) and the formula (6) are respectively expressed as follows:

wherein the formula (5) and the formula (6) each represent that a constraint boundary is adjusted according to a distance between a current state and the constraint boundary, so that the state outside the constraint range is transferred to be within the state constraint in the process of iterating, and a value range of τ is:

calculating intermediate variables as per a formula (7) expressed as follows:

calculating the variable rij(k), comprising: sending, by the i-th node, the intermediate variables (Δxi(k), Δxi(k)) and the state xi(k) to adjacent nodes, and receiving intermediate variables and a state of adjacent j-th node for calculating the variable rij (k) as per a formula (8) expressed as follows:

where σij(k)=xj(k)−xi(k), j∈Ni(k); rij(k) is configured to prevent xi(k) from exceeding the constraint boundaries in the process of iterating when transferring along a direction approaching the constraint boundaries to thereby ensure that a convergence result meets the state constraint.