## 2 Experimental Procedures and Research Methods
### Design of BPNN Optimized by QPSO  
The initial weights and thresholds of BPNN training were randomly selected and easily fell into the local optimal solution, resulting in poor prediction performance of BPNN. Therefore, the QPSO algorithm was used to optimize the weights and thresholds of BPNN, which is called the BPNN-QPSO neural network model. In the QPSO algorithm, the particle swarm is an \(ND\times NP\) two-dimensional array, denoted as \((\text{O}_{id,jp})_{ND\times NP}\). The position vector and velocity vector of the \(id\)-dimensional particle in the \(z\)th update are \(x_{id}^{z}=(x_{id,1}^{z},x_{id,2}^{z},\cdots x_{id,NP}^{z}\) and \(\mathbf{v}_{id}^{z}=(v_{id,1}^{z},v_{id,2}^{z},\cdots v_{id,NP}^{z})\), respectively. The optimal value of the particle after the \(z\)th update is called the best personal position (\(\textit{Pbest}_{id,jp}^{z}\)). If the best personal position is also the optimal value in the particle swarm, it is also called the best global position (\(\textit{Gbest}_{id,jp}^{z}\)). The flow chart of the optimization algorithm is shown in Figure 5.  
* Step 1: Extract the weight and threshold of BPNN, initialize the particle swarm, position vector (\(\mathbf{x}_{id}\)), and velocity vector (\(\mathbf{v}_{id}\)). * Step 2: Calculate the fitness of the initial particles (_fitness_) The standard error (\(\sigma\)) is selected as fitness value (_fitness_), as shown in Equation (8), to evaluate the performance of BPNN with different weights and thresholds, i.e., to evaluate the best personal position and the global position of particles [37]. Where \(\varepsilon_{k}\) is the relative error. \[\sigma=\sqrt{\frac{\varepsilon_{1}^{2}+\varepsilon_{2}^{2}+\varepsilon_{3}^{2 }+\cdots\cdots\varepsilon_{k}^{2}}{N}},\ \varepsilon_{k}=\left|\frac{d_{k}-net_{k}}{net_{k}}\right|\times 1 00(\%)\] (8)  
Figure 4: The structure of BPNN. * Step 3: Compare the best personal position (\(\textit{Pbest}_{id,jp}^{z+1}\)) of the particle in the (\(z+1\))th update with the best personal position (\(\textit{Pbest}_{id,jp}^{z}\)) of the \(z\)th update. If (\(z+1\))th _fitness_ is better than \(z\)th _fitness_, then (\(\textit{Pbest}_{id,jp}^{z}\)) is replaced by (\(\textit{Pbest}_{id,jp}^{z+1}\)). * Step 4: Compare the best global position (\(\textit{Gbest}_{id,jp}^{z+1}\)) of the particle in the (\(z+1\))th update with the best global position (\(\textit{Gbest}_{id,jp}^{z}\)) of the \(z\)th update. If (\(z+1\))th _fitness_ is better than \(z\)th _fitness_, then (\(\textit{Gbest}_{id,jp}^{z}\)) is replaced by (\(\textit{Gbest}_{id,jp}^{z+1}\)). * Step 5: Update \(\textit{Mbest}_{jp}^{z}\), \(\vec{v}_{id}^{z}\) and \(\vec{x}_{id}^{z}\) of particles. \(\textit{Mbest}_{jp}^{z}\) is the average value of \(\textit{Pbest}_{id,jp}^{z}\) of all particles, which is the center of the best personal position of all particles, called the mean best position, as shown in Equation (9). The particles update their velocity vector and position vector according to Equations (10)-(12) in the QPSO algorithm: \[\textit{Mbest}_{jp}^{z}=\left(\frac{1}{ND}\sum_{id=1}^{ND}\textit{Pbest}_{id,1 }^{z},\frac{1}{ND}\sum_{id=1}^{ND}\textit{Pbest}_{id,2}^{z},\cdots\frac{1}{ND }\sum_{id=1}^{ND}\textit{Pbest}_{id,N_{p}}^{z}\right)\] (9) \[\left\{ \begin{array}{c}v_{id}^{z+1}=\varphi\cdot\textit{Pbest}_{id}^{z}+(1 -\varphi)\textit{Gbest}_{id}^{z}\\ \textit{Pbest}_{id}^{z}=\left(\textit{Pbest}_{id,1}^{z},\textit{Pbest}_{id,2 }^{z},\cdots\textit{Pbest}_{id,NP}^{z}\right)\\ \textit{Gbest}_{id}^{z}=\left(\textit{Gbest}_{id,1}^{z},\textit{Gbest}_{id,2 }^{z},\cdots\textit{Gbest}_{id,NP}^{z}\right)\end{array} \right.\] (10) \[x_{id}^{z+1}=v_{id}^{z+1}+rand(u)\times\alpha\Big{|}\textit{Mbest}_{jp}^{z}-x_ {id}^{z}\Big{|}\ln\left(\frac{1}{u}\right)\] (11) \[rand(u)=\left\{ \begin{array}{c}+1,u\geq 0.5\\ -1,u<0.5\end{array} \right.\] (12) where \(\varphi\) and \(u\) are a sequence of random numbers on the uniform distribution within (0,1), and \(\alpha\) is a contraction factor.