1. Introduction
In the background of the Energy Internet and the “dual carbon” goals, breaking the isolated system of traditional energy systems and constructing a multi-energy microgrid (MEMG) with multiple energy couplings has become a primary measure for governments and enterprises worldwide to promote energy transition [
1]. Meanwhile, as an important energy technology with the capability of energy time shifting and spatiotemporal decoupling, energy storage plays a significant role in providing flexible regulation during the planning and construction of MEMGs [
2]. Therefore, it is of great significance to explore the planning and construction strategies of MEMGs supported by energy storage resources [
3].
Currently, many studies have focused on MEMG planning methods that take energy storage resources into account [
4]. For example, a capacity planning model for MEMG energy conversion equipment and energy storage devices considering dynamic energy conversion efficiency coefficients is proposed in [
5]. A capacity planning method for MEMG equipment is proposed in [
6], which considers the investment constraints of energy conversion and storage devices to reduce the planning cost of MEMG. The energy hub modeling method of MEMG is improved using directed acyclic graph theory in [
7], and based on this, a low-carbon capacity planning model for MEMG energy stations with carbon capture devices and energy storage devices is proposed, aiming to minimize the total planning and operating cost of MEMG. An optimal planning method for hydrogen-based rural MEMG is proposed in [
8], which considers multiple time scales and multi-energy storage devices, aiming to achieve economic and reliable carbon neutral energy supply in rural communities. A device location and capacity optimization method for MEMG is proposed in [
9], which considers various energy facilities such as heat pumps, electrical energy storage, and thermal energy storage, aiming to reduce the overall planning cost of MEMG. In addition to the aforementioned physical energy storage, some studies have also explored planning methods for MEMG considering virtual energy storage. For example, the impact of demand response on the economic efficiency of MEMG equipment planning is analyzed in [
10], and an optimal planning model for MEMG that considers demand response as virtual energy storage is proposed to reduce MEMG’s total planning costs and carbon emissions. A two-stage robust optimization planning model for MEMG considering demand response is constructed in [
11], and the optimal equipment configuration strategy for MEMG is obtained by iteratively optimizing the upper-level planning model and the lower-level operation model. A device capacity planning model for MEMG that considers virtual heat storage in heating networks is proposed in [
12], and it is verified through case studies that virtual heat storage in heating networks can improve the economic efficiency of MEMG’s planning while improving the consumption level of renewable energy. A coordinated siting and sizing model for MEMG considering generalized energy storage is established in [
13], where the physical energy storage device and virtual heat energy storage are simultaneously considered in this paper to reduce the total planning and operating costs and carbon emissions.
The aforementioned studies focus on scenarios where energy storage is independently used by individual MEMGs. In these cases, the investment and operational costs of energy storage are relatively high, and the utilization rate of the equipment is relatively low, hindering the further development and application of energy storage in MEMGs. Shared energy storage (SES) is a new business model for energy storage that has emerged in the context of the current sharing economy, whose core idea is to temporarily transfer the usage rights of idle energy storage resources from energy storage owners to lessees for a certain fee. Under the background of the increasing renewable energy penetration, sharing energy storage among multiple MEMGs can effectively improve the economic efficiency of MEMG operations, the renewable energy consumption rate, and the utilization rate of energy storage devices. Since Qinghai province initiated China’s first SES pilot in 2019, SES has gradually become a hot research topic in China’s energy storage field in recent years. Against this background, some scholars have explored optimization planning methods for SES, aiming to improve the utilization rate and economic efficiency of energy storage devices. For example, a bi-level planning model for SES operators (SESOs) and users is constructed in [
14], aiming to increase the revenue for both parties. An SES capacity planning model for wind farms and photovoltaic power stations in distribution networks is constructed in [
15]. A business model in which SESOs provide ancillary services to renewable energy generators is proposed in [
16], aiming to enhance the investment enthusiasm of SESOs. The abovementioned studies focus on large-scale centralized SES planning methods, some studies have also addressed optimization planning methods for distributed SES on the user side. For instance, an optimal configuration model for user-side SES under a peer-to-peer energy sharing framework is constructed in [
17]. An optimal configuration method for electricity retailers is proposed in [
18], which aims to reduce electricity procurement costs for retailers by analyzing the investment benefits of SES. An optimal configuration model for SES in communities with multiple buildings is proposed in [
19], aiming to reduce the SES planning costs for communities. A community-oriented SES planning framework is proposed in [
20], and the effectiveness of SES in reducing energy storage configuration costs is verified by case studies. In the context of the Energy Internet, some scholars have also explored planning methods for MEMGs considering SES. For example, a capacity planning model for SES jointly built by multiple industrial MEMGs is constructed in [
21], and the impact of SES on the economic efficiency of industrial MEMGs’ planning under different capacity ratios is investigated. An SES planning and scheduling model for MEMG is established in [
22], and a non-dominated sorting equilibrium optimizer algorithm is introduced to prevent the Pareto solution set from getting trapped in local optima and to guarantee the successful execution of the proposed benefit distribution mechanism.
The aforementioned studies have made significant contributions to the development of MEMG and SES, but the following shortcomings remain: (1) Most existing research focuses on the independent planning method of MEMG or SESO, with less attention paid to the collaborative planning of multiple MEMG alliances and SESO. In fact, the SES leasing and scheduling strategies of MEMG alliances are based on SESO’s SES pricing, while the SES leasing strategies of MEMG alliances also affect SESO’s physical storage planning, scheduling, and pricing strategies. Therefore, it is urgent to study the collaborative planning strategies of SESO and MEMG alliances. (2) Existing research on SES capacity planning mostly does not consider the impact of storage lifespan degradation. In reality, ignoring the degradation of storage lifespan may accelerate the depreciation of equipment residual value, thereby indirectly increasing SESO’s planning costs. Thus, it is crucial to account for the impact of storage lifespan degradation in SESO’s planning process.
Hence, a joint planning method of SESO and MEMG alliances based on a dynamic game with perfect information is proposed in this paper to address the abovementioned issues. The contributions of this paper are summarized as follows:
- (1)
An upper-level model for energy storage capacity configuration and pricing planning of SESO is proposed, aiming to maximize the total planning and operational income of SESO. Energy storage operation constraints that consider storage lifespan degradation are introduced to improve the lifespan of energy storage devices, thereby reducing SESO’s energy storage configuration costs and subsequently increasing SESO’s total planning income.
- (2)
A lower-level model for the optimal configuration of the MEMG alliance considering SES is proposed, aiming to minimize the total planning and operational costs of the MEMG alliance. In the lower-level model, each MEMG independently invests in ECDs, while jointly leasing and utilizing SES. The planning and operational costs of SES are allocated using the Shapley value method. Compared with the planning model where each MEMG independently invests in energy storage devices, the proposed lower-level model can achieve lower total planning and operational costs for the MEMG alliance.
- (3)
A solving algorithm based on the dynamic game theory with perfect information and the backward induction method is proposed to obtain the Nash equilibrium solution of the proposed bi-level optimization models. A dynamic non-cooperative game framework with a leader–follower hierarchical structure is constructed for the proposed bi-level optimization models. On this basis, the bi-level optimization model is iteratively solved using the dynamic game theory with perfect information, and the subgame perfect Nash equilibrium solution is obtained using the backward induction method. The most beneficial collaborative planning strategy for both SESO and the MEMG alliance can be obtained through the proposed solving algorithm.
2. Upper-Level Model for Energy Storage Capacity Configuration and Pricing Strategy Planning of SESO
The upper-level model is from the perspective of SESO, and reasonable strategies for SES capacity configuration electricity, charging and discharging, electricity purchasing and selling with the upper-level grid, and SES rental service pricing are formulated based on the SES leasing demand of MEMGs. It should be noted that the charging and discharging power of SES is transmitted through the distribution network, so SESO needs to pay a certain transmission fee to the distribution network.
The objective function of the upper-level model is to maximize SESO’s net income, which is the SES service income over the planning period minus the physical energy storage investment cost, physical energy storage maintenance cost, and electricity trading cost with the grid company, as shown in Equation (1). It should be noted that the occurrence probability of scenario
s, i.e.,
xs, is introduced in the following Equation (1) to realize the modeling of the stochastic programming method, thereby dealing with uncertainty in the planning process:
where
Y,
S, and
T are the number of planning years, typical planning scenarios, and scheduling time, respectively;
is the discount coefficient;
and
are the rental unit prices for the energy capacity and power capacity of the SES service used by MEMGs’ alliance, respectively;
and
are the energy capacity and power capacity of SES rented by MEMGs’ alliance from SESO at year
y, respectively;
and
are the unit investment costs of the energy capacity and power capacity of the physical energy storage devices invested by SESO, respectively;
and
are the energy capacity and power capacity of the physical energy storage devices invested by SESO, respectively;
is the unit maintenance cost of the physical energy storage devices;
and
are the charging and discharging power of SESO’s physical energy storage devices at time
t in scenario
s of year
y, respectively;
is the interval between adjacent scheduling times;
and
are the unit price of electricity sales and repurchases by power grids at time
t in scenario
s of year
y, respectively; and
and
are the power purchased from and sold to the grid company by SESO at time
t in scenario
s of year
y, respectively.
In the upper-level model, a lifespan degradation model for physical energy storage devices is constructed based on the rainflow counting method, aiming to extend the lifespan of SES by limiting the daily charge–discharge cycle count, thereby indirectly reducing SESO’s SES investment costs. The rainflow counting method is a commonly used analysis method in the field of material fatigue life research and is widely applied in engineering practice to analyze the relationship between discharge depth and lifespan of SES [
23]. The rainflow counting method counts cycles based on the nonlinear relationship between the state of charge (SOC) of SES and scheduling times, thereby reflecting the relationship between discharge depth and lifespan of SES. Specifically, the steps to calculate the charge–discharge cycle count and discharge depth of SES within a scheduling period using the rainflow counting method are as follows:
First, the relationship curve between SOC and scheduling times of SES is rotated 90 degrees clockwise, assuming that the virtual rainflow starts to flow down from the initial point. Then, the rainflow drips vertically when it reaches each “eave”, i.e., the peaks and valleys of SOC, until it falls to the next roof that is larger than the “peak eave” or smaller than the “valley eave”, forming a charge–discharge cycle and recording the SOC values at the start and end points of the rainflow path. Next, the rainflow continues to flow down from the last charge–discharge cycle until it reaches the next “eave”, repeating the abovementioned steps until it flows to the end of the scheduling period. Finally, based on each rainflow path and its start and end SOC values, the discharge depth of each rainflow path is calculated.
Taking the relationship curve between SOC and scheduling times of SES shown in
Figure 1 as an example, it can be concluded according to the above steps that the first charge–discharge cycle is B-C-B′, with a discharge depth of 16%; the second charge–discharge cycle is E-F-E′, with a discharge depth of 10%; the third partial charge–discharge cycle is A-B-B′-D, with a discharge depth of 58%; and the fourth partial charge–discharge cycle is D-E-E′-G, with a discharge depth of 40%.
Based on the abovementioned principles, it is evident that the process of calculating discharge depth and charge–discharge cycle count using the rainflow counting method involves solving nonconvex and nonlinear problems, such as extracting extreme values, which are challenging to express analytically. Therefore, the discharge level of ESD at each scheduling time is discretized in this paper, and equivalently converted into the maximum charge–discharge cycle count at a 100% discharge depth. By limiting the sum of the maximum charge–discharge cycle counts at each scheduling time within the scheduling period, a corresponding charge–discharge strategy of SES is formulated to reduce the degradation of SES lifespan. It should be noted that the depth of discharge of the ESD at each time step is normalized to a 100% depth of discharge. The reason is that the lifespan of the storage at a 100% depth of discharge is a standard for comparison. By converting the depth of discharge at each time step to a 100% depth of discharge, the following equations can be utilized to calculate the lifespan of the ESD at the current depth of discharge. The specific modeling steps are given below.
First, according to industry test data on the discharge depth and maximum charge–discharge cycle count of SES, a power function is used for curve fitting [
24], deriving the relationship between the discharge depth and the maximum charge–discharge cycle count at the time, expressed as
where
is the actual discharge depth value of SES at time
t;
is the number of charge–discharge cycles that SES can perform at discharge depth
until reaching the maximum lifespan;
is the number of charge–discharge cycles that SES can perform at a 100% discharge depth until reaching the maximum lifespan; and
is the fitting parameter.
Then,
is converted to the equivalent maximum cycle count
at a 100% discharge depth.
Equation (3) is a nonconvex and nonlinear power function, making it difficult to construct a linear optimization model. Therefore, the power function is piecewise linearized in this section, as shown in Equation (4), with the addition of two constraints shown in Equations (5) and (6), as follows:
where
K is the number of linear segments;
and
are the slope and intercept of the
kth segment, respectively;
is an auxiliary binary variable with a value of 1 indicating that the discharge depth of SES at time
t in scenario
s of year
y falls within the
kth segment (thus,
=
), and a value of 0 indicating that the discharge depth does not fall within the
kth segment (thus,
= 0);
and
are the upper and lower bounds of the discharge depth for the
kth segment, respectively; and
is an auxiliary binary variable reflecting the charge–discharge cycle status of SES at time
t in scenario
s of year
y, with a value of 1 indicating the completion of a charge–discharge cycle at that time. It should be noted that, for Equation (4), according to the constraints shown in Equation (5), when
= 0,
= 0, and when
= 1,
=
. Therefore, Equation (4) can be equivalently relaxed to Equation (7).
Since the discharge depth
of SES is a decision variable influenced by the energy storage charge–discharge strategy, it is subject to constraints related to SOC and the charge–discharge cycle status, as shown in (8) and (9). More specifically, Equation (8) represents the relationship between the discharge depth and SOC of SES; Equation (9) represents the relationship between the charge–discharge cycle status and the charge–discharge state, which indicates that a charge–discharge cycle is completed if and only if SES transitions from a discharge state to a charge state, as follows:
where
is the remaining energy capacity of SES at time
t − 1 in scenario
s of year
y, and
and
are auxiliary binary variables reflecting the charging state of SES at times
t and
t − 1, respectively. For Equation (8), the big-M method can be employed to equivalently relax it into the following Equations (10) and (11):
where
M is a positive real number with a sufficiently large numerical value.
Based on the aforementioned modeling, as long as the sum of the maximum charge–discharge cycle counts at each scheduling time is limited to within the maximum daily average charge–discharge cycle count derived from the expected lifespan of SES, SES can be effectively ensured to reach its expected lifespan, as shown in the following Equation (12):
where
is the maximum daily average charge–discharge cycle count derived from the expected lifespan of SES, where
,
is the expected lifespan of SES.
For the nonlinear term in Equation (12), it can be converted to Equation (13) by multiplying both sides of the inequality by
. Then, the product of binary and continuous variables, i.e.,
, can be handled using the big-M method, thereby equivalently converting Equation (12) into a linear constraint.
At this point, the energy storage lifespan degradation model based on the rainflow counting method has been constructed. Incorporating the basic operational constraints of SES shown in Equations (14)–(17) [
25], they together form the operational constraints of SES considering lifespan degradation. More specifically, Equation (14) represents the constraint relationship between SOC and the charge–discharge power of SES; Equation (15) represents the constraint relationship between SOC and the invested energy capacity of SES; and Equations (16) and (17) represent the constraint relationships between the charge–discharge power and the invested power capacity of SES, as follows:
where
is the remaining energy capacity of SES at time
t in scenario
s of year
y;
and
are the charging and discharging loss coefficients of SES, respectively; and
and
are the upper limits of the energy capacity and power capacity investment for SES, respectively.
The constraints of the upper-level model also include the upper and lower bound constraints of the rental pricing for the energy capacity and power capacity of SES, as shown in Equations (18) and (19), respectively; additionally, the constraint shown in Equation (20) represents that the SES charge–discharge demand of the MEMG alliance needs to be met, as follows:
where
and
are the upper and lower bounds of the rental pricing for the energy capacity of SES set by SESO, respectively;
and
are the upper and lower bounds of the rental pricing for the power capacity of SES set by SESO, respectively;
N is the number of MEMGs; and
and
are the charging and discharging power of the
nth MEMG at time
t in scenario
s of year
y, respectively.
3. Lower-Level Model for Optimal Configuration of MEMGs’ Alliance Considering SES
The lower-level model is designed from the perspective of MEMGs. Each MEMG independently plans its own ECDs based on its load demand and formulates its energy scheduling strategy. The MEMG alliance also develops a reasonable SES rental strategy. The following reasonable assumptions are made for the MEMG alliance: (1) The MEMG alliance needs to elect a representative to summarize the daily energy storage charge and discharge plans of each MEMG and notify SESO. (2) The MEMG alliance can rent SES services either annually or per planning period, based on negotiations between the MEMG alliance and SESO. In this paper, it is assumed that MEMGs rent SES services annually.
Each MEMG’s energy supply sources include the electric power company (EPC), natural gas company (NGC), and distributed renewable energy sources (RES); the ECDs include electric heaters (EH), gas boilers (GB), and combined heat and power units (CHP); and the load types of MEMG include electric load and thermal load. Assuming that
N MEMGs form an alliance to jointly use SES services, the objective function of the lower-level model, i.e.,
, is to minimize the sum of the life cycle costs of each MEMG in the alliance, as shown in the following Equation (21):
where
is the residual value of the ECDs invested by the
nth MEMG in year
y [
26];
,
, and
are the capacities of the newly built EH, GB, and CHP by the
nth MEMG in year
y, respectively;
,
, and
are the unit capacity investment costs of EH, GB, and CHP, respectively;
is the unit price of gas sold by the NGC at time
t in scenario
s of year
y; and
and
are the power purchased from the EPC and the volume of gas purchased from the NGC by the
nth MEMG at time
t in scenario
s of year
y, respectively.
The constraints of the lower-level model include the power balance constraints shown in Equations (22)–(24), the SES operational constraints shown in Equations (25)–(28), the ECD operational power upper limit constraint shown in Equation (29), and the energy conversion constraints shown in Equations (30)–(32), as follows:
where
and
are the electric and thermal load powers of the
nth MEMG at time
t in scenario
s of year
y, respectively;
and
are the electric powers supplied by the CHP and the RES of the
nth MEMG at time
t in scenario
s of year
y, respectively;
is the electric power consumed by the EH of the
nth MEMG at time
t in scenario
s of year
y;
,
, and
are the thermal powers supplied by the CHP, EH, and GB of the
nth MEMG at time
t in scenario
s of year
y, respectively;
and
are the natural gas consumption by the CHP and GB of the
nth MEMG at time
t in scenario
s of year
y, respectively;
and
are the remaining energy capacities of SES leased by the
nth MEMG at time
t in scenario
s of year
y, respectively;
is the output power of ECDs, where
E denotes CHP, GB, and EH;
is the newly installed capacity of
E of the
nth MEMG in year
x;
and
are the energy capacity and power capacity leased by the
nth MEMG in year
y for SES, respectively;
,
,
, and
are the electrical-to-thermal conversion coefficient of CHP, electrical efficiency coefficient of CHP, thermal-to-electrical conversion coefficient of EH, and gas-to-thermal conversion coefficient of GB, respectively.
It should be noted that, after the lower-level model determines the ECD capacity planning and SES service leasing strategy for MEMGs within the planning period, the SES service leasing costs generated by the alliance during the planning period need to be reasonably allocated to maintain the stability of the alliance. The cooperative game theory points out that the differences in marginal benefits among alliance members will result in varying bargaining power during the cost-sharing process. Therefore, the marginal benefit contribution index is introduced in this section to assess the position of each MEMG in SES service cost-sharing negotiations. The Shapley value method achieves a fair distribution of SES service leasing costs generated by the alliance by measuring each MEMG’s marginal benefit contribution in all potential sub-alliances. Hence, the SES service cost allocation result for the
ith MEMG in the coalition
N, i.e.,
, can be obtained by the Shapley value method, which is expressed as
where
S is the potential sub-alliances in the MEMGs alliance
N; |
S| is the number of MEMGs in
S;
and
are the alliance income of the sub-alliance
S and the sub-alliance
S −
i (which excludes the
ith MEMG), respectively; and
−
is the marginal benefit brought to the sub-alliance
S by the joining of the
ith MEMG.
4. Solving Algorithm for Bi-Level Optimization Models Based on the Dynamic Game Theory with Perfect Information and the Backward Induction Method
For the previously constructed bi-level optimization model, the SES leasing and scheduling strategy of the MEMG alliance is based on SESO’s SES pricing strategy. Therefore, SESO’s SES pricing strategy indirectly affects each MEMG’s ECD planning and scheduling strategy. Additionally, the SES leasing strategy of the MEMG alliance will also impact SESO’s physical energy storage planning, scheduling, and pricing strategies. In summary, this process conforms to a dynamic non-cooperative game framework with a leader–follower hierarchical structure. Furthermore, the SES pricing strategy transmitted from the upper-level model to the lower-level model and the SES leasing and scheduling strategy feedback from the lower-level model are mutually public information. Therefore, the proposed joint planning bi-level optimization problem for SESO and MEMGs is a dynamic game with perfect information. In this context, a bi-level optimization model solving algorithm based on the dynamic game theory with perfect information and the backward induction method is proposed in this paper. More specifically, the bi-level optimization model is iteratively solved using the dynamic game theory with perfect information, and the subgame perfect Nash equilibrium solution of the bi-level optimization model is obtained using the backward induction method. The specific steps are as follows.
Let SESO be the leader and the MEMG alliance be the follower in the bi-level model, establishing a dynamic game with a perfect information framework, including the following elements: participants {SESO, MEMGs}, strategies {, , , }, and incomes {, }. ={, } represents SESO’s SES energy capacity and power capacity pricing strategies; = {, , , , , } (y = 1, 2, …, Y; s = 1, 2, …, S; t = 1, 2, …, T) represents SESO’s planning and scheduling strategies; = {, , , } (n = 1, 2, …, N; y = 1, 2, …, Y; s = 1, 2, …, S; t = 1, 2, …, T) represents MEMGs’ SES service usage and charge–discharge strategies; = {, , , , , , , } represents MEMGs’ planning and scheduling strategies; and and represent the incomes for SESO and MEMGs, respectively.
When the follower makes an optimal response based on the leader’s strategy, and the leader accepts this response without further changing its strategy, the dynamic game with perfect information reaches Nash equilibrium. Assuming that NE* = {
,
,
,
} is the Nash equilibrium solution of the proposed bi-level model in this paper, it should satisfy the following conditions:
According to the dynamic games theory with perfect information, the actions of participants are sequential, with the leader and follower making strategies in an interdependent manner in turn. The original dynamic game with perfect information can be broken down into a series of smaller subgame problems, making it easier to intuitively describe the sequence of actions, strategy information, and incomes of the game participants. Therefore, a game tree is used in this paper to illustrate the dynamic game process with perfect information. In the game tree, nodes representing sequential game actions are termed as nodes of the game tree, including initial nodes, decision nodes, and terminal nodes. It is evident that decision nodes can have multiple direct successor nodes, indicating that the participant at these decision nodes can choose from multiple strategies; terminal nodes indicate the end of the game and, thus, have no direct successors. Tracing back from a terminal node to its direct predecessor and repeating this process to the initial node, a path in the game tree can be obtained, representing a combination of action plans among the participants.
The game tree of the dynamic game with perfect information for the proposed bi-level optimization model is shown in
Figure 2. In
Figure 2, decision node 1 represents SESO setting the energy capacity and power capacity pricing strategies of SES; based on the strategy from decision node 1, decision node 2 represents the MEMG alliance formulating the SES service usage and charge–discharge strategies, as well as the planning and scheduling strategies for each MEMG; decision node 3 represents SESO receiving the response from the MEMG alliance and then accordingly formulating its planning and scheduling strategies; the terminal node represents the final incomes for SESO and MEMGs based on the strategies from decision nodes 1–3.
N1,
N2, and
N3 denote the number of subgame nodes for decision nodes 1, 2, and 3, respectively.
Based on the aforementioned concepts, a subgame is defined as follows: Starting from a decision node in the game tree and including all subsequent nodes of that decision node, every decision node in the game tree can serve as an initial node, and together with all its subsequent nodes, they form a subgame. For the dynamic game with a perfect information problem of the proposed bi-level optimization model, if the strategy combination {
,
,
,
} is the Nash equilibrium of the original dynamic game with a perfect information problem, and it is also the subgame Nash equilibrium for each subgame, then this strategy combination is called the subgame perfect Nash equilibrium (SPNE) of the proposed model [
27]. Reference [
28] has proven that a dynamic game with perfect information has a subgame perfect Nash equilibrium, and it has also proven that the backward induction method can be used to find the subgame perfect Nash equilibrium in a dynamic game with perfect information. Therefore, the backward induction method is employed in this paper to obtain the subgame perfect Nash equilibrium [
29]. The backward induction method-based solving algorithm for the proposed bi-level optimization model is shown in Algorithm 1. The process of backward induction from a decision node to its direct predecessor can be solved using Gurobi 11.0 or CPLEX 12.9 commercial solvers.
Algorithm 1 Backward induction method-based solving algorithm for the proposed bi-level optimization model |
1: Input: Parameters of the proposed bi-level optimization model in the first stage. |
2: Initialize: i=1, j=1, k=1, S1flag_1=0, S1flag_2=0, S1flag _3=0, enter the initial node. |
3: Repeat |
4: Enter decision node 1; the upper-level model provides the strategy , and then transmits it to decision node 2. |
5: Repeat |
6: Enter decision node 2; the lower-level model provides the strategies and , and then transmits it to decision node 3. |
7: Repeat |
8: Enter decision node 3; the upper-level model provides the strategy , and then transmits the strategies {, , , } to the terminal node. |
9: Enter the terminal node, solving the upper- and lower-level models proposed in Section 2 and Section 3 to obtain {, }. |
10: Backward induction to decision node 2, and according to (34), determine whether {, , , } is the Nash equilibrium solution of the subgame tree starting at {, }. If it is, set S1flag_3=1; otherwise, set k=k + 1. |
11: Until k > N3 or S1flag_3=1 |
12: Backward induction to decision node 1, and according to (34), determine whether {, , , } is the Nash equilibrium solution of the subgame tree starting at . If it is, set S1flag_2=1; otherwise, set j=j + 1. |
13: Until j > N2 or S1flag_2=1 |
14: Backward induction to initial node 1, and according to (34), determine whether {, , , } is the Nash equilibrium solution of the initial game tree. If it is, set SPNE*={, , , }, S1flag_1=1; otherwise, set i=i + 1. |
15: Until i > N1 or S1flag_1=1 |
16: Output: SPNE* |
17: End |