4.2.3. Updating Policies

Step 1. Before the update operation, individuals are selected based on the elite strategy to increase the running speed of the algorithm. That is, the fitness of individuals generated in the population is first sorted from largest to smallest, and the top 50% of individuals are retained.

Step 2. DIWPSO is used for updating. Although the standard PSO has a fast convergence speed in the early stage, it is slow in later stages and easily converges locally. Therefore, this algorithm is improved from the perspective of the inertia weight. The inertia weight *ω* indicates the extent to which the original speed is retained; if *ω* is larger, the global search ability is stronger, and if *ω* is small, the local search ability is strong.

The update strategy is as follows: in the position vector, Equations (26) and (27) are used to update the particle velocity and position for the continuous factor *visk*(*t*):

$$v\_{id}(\tau+1) = v\_{id}(\tau)\omega + c\_1 r\_1(\tau)[p\_{id}(\tau) - x\_{id}(\tau)] + c\_2 r\_2(\tau)[\mathbf{g}\_d(\tau) - x\_{id}(\tau)] \tag{26}$$

$$\mathbf{x}\_{id}(\tau+1) = \mathbf{x}\_{id}(\tau) + \upsilon\_{id}(\tau+1) \tag{27}$$

$$
\omega = \omega\_{\rm min} + (\omega\_{\rm max} - \omega\_{\rm min}) \ast e^{-\frac{\tau}{\hbar \max}} + \sigma \ast \text{beta} \,\text{arcnd}(p, q) \tag{28}
$$

where *τ* represents the current iteration number; *τmax* represents the maximum number of iterations; *ωmax* represents the maximum inertia weight, which is set to 0.9; *ωmin* represents the minimum inertia weight, which is 0.1; *σ* is the inertia adjustment factor, which is 0.1; *p* = 1, *q* = 3; *c*<sup>1</sup> and *c*<sup>2</sup> are learning factors; *r*<sup>1</sup> and *r*<sup>2</sup> are uniform random numbers between [0,1]; *xid*(*τ*) and *vid*(*τ*) represent the position and velocity of the *d* dimension elements, respectively, of the *i* particle after the *τ* iteration; *pid*(*τ*) represents the individual optimal position of the *i* particle in the *d* dimension; and *gd*(*τ*) represents the global optimal position of all particles in the *d* dimension.

*ω* in the update strategy is an improved strategy for the dynamic adjustment of inertia weight [33], and the exponential function is used to control the change in inertia weight *ω*. With an increase in the number of iterations, *e*−*τ*/*τmax* decreases nonlinearly; thus, *ω* can ensure the breadth of global search in the early stage and gradually decrease in the later stage to improve the ability of the local search and ensure its accuracy. Betarnd is a random number generator in MATLAB that can generate random numbers in line with the beta distribution. In addition, an inertia adjustment factor *σ* was added to control the deviation of the inertia weight, to make the adjustment more reasonable.

Step 3. Particle evaluation. To avoid generating infeasible particle positions and excessive velocities during the iteration, they must be within the corresponding limits.

$$\upsilon\_{\rm id} = \begin{cases} \upsilon\_{\rm max}, \upsilon\_{\rm id} > \upsilon\_{\rm max} \\ \upsilon\_{\rm min}, \upsilon\_{\rm id} < \upsilon\_{\rm min} \end{cases}, \chi\_{\rm id} = \begin{cases} \chi\_{\rm max}, \chi\_{\rm id} > \chi\_{\rm max} \\ \chi\_{\rm min}, \chi\_{\rm id} < \chi\_{\rm min} \end{cases} \tag{29}$$

Step 4. Particle adjustment. Since the fitness function is designed with the belief that larger is better, the individual *Pi*(*τ*) and the global *G*(*τ*) optimums are updated by calculating the fitness of the particles.

Step 4.1. For the individual optimum *Pi*(*τ*), if *Fitness*[*Yi*(*τ*)] > *Fitness*[*Pi*(*τ* − 1)], update *Pi*(*τ*) = *Yi*(*τ*); otherwise, maintain the original value.

Step 4.2. For the global optimum *G*(*τ*), if *Fitness*[*Pi*(*τ*)] > *Fitness*[*G*(*τ* − 1)], update *G*(*τ*) = *Pi*(*τ*); otherwise, maintain the original value, namely *G*(*τ*) = *G*(*τ* − 1).

Step 5. Premature particle determination. To judge the convergence degree of the particles, the population fitness variance [34] was introduced as the judgment mechanism of particle prematurity.

$$\delta^2 = \frac{1}{N} \sum\_{i=1}^{N} \left( \frac{f\_i - f\_{avg}}{f} \right)^2 \tag{30}$$

$$f = \begin{cases} \max|f\_i - f\_{avg}|, \max|f\_i - f\_{avg}| > 1\\ 1, \text{ otherwise} \end{cases} \tag{31}$$

where *δ*<sup>2</sup> is the variance of population fitness; the larger *δ*<sup>2</sup> is, the better the population diversity, and vice versa. *fi* is the fitness of the *i* particle; *favg* is the average fitness of the population, and *f* is the normalization factor, which limits the size of *δ*2.

A population fitness judgment threshold *δ*<sup>2</sup> *<sup>T</sup>* is selected for premature judgment: when *δ*<sup>2</sup> < *δ*<sup>2</sup> *<sup>T</sup>*, the particle enters premature convergence. *<sup>δ</sup>*<sup>2</sup> *<sup>T</sup>* is generally much smaller than the fitness variance of the initial population; *δ*<sup>2</sup> *<sup>T</sup>* = 0.001 is taken here.

Step 6. The mutation operation exists to improve the ability of the algorithm to jump out of premature convergence, ensure the diversity of the population, and keep the algorithm from falling into local convergence in the later stage to stop searching for a better solution. The mutation mechanism of the differential evolution algorithm is used to mutate the identified premature particles.

$$V\_i(\tau+1) = \mathbf{x}\_{r1}(\tau) + \eta \left[\mathbf{x}\_{r2}(\tau) - \mathbf{x}\_{r3}(\tau)\right] \tag{32}$$

*r*1,*r*2,*r*3 ∈ (1, 2, . . . , *N*) is a random number and *r*1 = *r*2 = *r*3 = *i*, and *η* is a scaling factor adjusted by adaptive strategy:

$$
\eta = \eta\_{\text{max}} - \tau (\eta\_{\text{max}} - \eta\_{\text{min}}) / \tau\_{\text{max}} \tag{33}
$$

where *η*max and *η*min are the upper and lower limits of the scaling factor, respectively.

#### *4.3. Overall Process Framework of the Algorithm*

The algorithm designed in this study is a two-layer GADS/DIWPSO hybrid algorithm. In the project scheduling problem, GADS is first used to initialize the feasible strategy and introduce it into lower-level planning. Then, DIWPSO is used to find the corresponding optimal solution of resource provisioning and the input to the upper planning is returned. Then, GADS is used to decode and generate the current optimal solution. This process is repeated until the upper optimal solution satisfies the stop condition. Through this dynamic interaction, the Stackelberg-Nash equilibrium strategy of the MRCPSP-MPMSP ensemble system is finally obtained.

The flow chart of this hybrid algorithm is shown in Figure 4, where the left part is the flow of GADS solving the upper-level project scheduling problem and the right part is the flow of DIWPSO solving the lower-level resource supply problem.

**Figure 4.** Flow chart of the hybrid algorithm.

#### **5. Practical Application**

The practical application and calculation test of a dam project verified the practicability and effectiveness of the proposed optimization method and provided decisionmaking guidance.

#### *5.1. Project Description*

In this study, a large hydropower project located in southeast China was considered as an application example. The project had a variety of hydraulic structures such as river dams, flood discharge structures, and hydraulic power generation systems. The river dam was a concrete double-curvature arch dam with a height of 610 m.

The concrete double-curvature arch dam construction project, which consists of 17 engineering activities, is the most important part. A flowchart is shown in Figure 5. Each activity has several optional modes, and each mode has a certain duration and resource demand. At the construction site, there are two large-scale resource demand points to allocate resources for each activity within the project, and the three resources required by the demand points are supplied by an external resource supplier with four resource supply points.

**Figure 5.** Construction flow chart of a concrete double-curvature arch dam.

#### *5.2. Data Collection and Setting*

#### 5.2.1. Project Scheduling Data Processing

To collect relevant data for this practical application, we conducted interviews and surveys with relevant construction companies. The construction process of a concrete double-curvature arch dam can be divided into 17 activities, among which there are three types of common resources. Table 2 shows the activities in which each demand point is responsible for providing resources, and the other necessary data are shown in Table 3.

**Table 2.** Demand point-project activity mapping table.


**Table 3.** Details on other parameters.


According to the preliminary data collected, the data of each activity in the project were processed in detail; specifically, uncertain variables were expressed in the form of random variables. The detailed processing data are shown in Table 4. In addition, the project planning period and available budget are *D* = 52 and *B* = 8510, respectively, the indirect cost of each time period is *c*<sup>0</sup> = 5.8, and the storage capacity of each period is *IC* = 300. The weights of the objective functions in the upper model were set to *μ*<sup>1</sup> = *μ*<sup>2</sup> = 0.5.


**Table 4.** Concrete double-curvature arch dam project activity details.
