## 3 Modeling
### Additive Three-Stage Network DEA Model
 \tag{12}\]  
The above fractional model (9) and (10) can be transformed into the following linear model (12)and (13) based on Equation (11), the result is as follows:  
\[\max{PSI}=\sum_{\delta_{i}=1}^{D_{i}}{\gamma_{\delta_{i}}}{z_{\delta_{i}k}}+\sum _{\delta_{i}=1}^{D_{i}}{\gamma_{\delta_{i}}}{z_{\delta_{i}k}}+\sum_{r=1}^{\delta _{i}}{\theta_{i}}{y_{\delta_{r}k}}\]  
By simplifying Equation (11), we can get:  
\[\sum_{i_{1}=1}^{m_{1}}{\pi_{i_{1}}}{x_{i_{1}k}}+\sum_{i_{2}=1}^{m _{2}}{\pi_{i_{2}}}{x_{i_{2}k}}+\sum_{\delta_{r_{2}=1}^{\delta_{2}}}{\theta^{2} }{{}_{r_{2}}y^{2}}{{}_{r_{2}k}}+\sum_{\delta_{r_{2}=1}^{\delta_{2}}}{\theta^{2} }{{}_{r_{2}}y^{2}}{{}_{r_{2}k}}+\sum_{\delta_{r_{2}=1}^{\delta_{2}}}{\gamma_{ \delta_{r_{2}}}}{z_{\delta_{r_{2}k}}}=1 \tag{112}\] \[\sum_{\delta_{i}=1}^{D_{i}}{\gamma_{\delta_{i}}}{z_{\delta_{i}k}} -\sum_{i_{2}=1}^{m_{1}}{\pi_{i_{2}}}{x_{i_{2}j}}-\sum_{\delta_{r_{2}=1}^{ \delta_{2}}}{\theta^{2}}{{}_{r_{2}}y^{2}}{{}_{r_{2}j}}\leq 0\] \[\sum_{\delta_{i}=1}^{\delta_{2}}{\gamma_{\delta_{i}}}{z_{\delta_{ i}d_{j}}}-\sum_{i_{2}=1}^{m_{1}}{\pi_{i_{2}}}{x_{i_{2}j}}-\sum_{\delta_{i}=1}^{D_{i}}{ \gamma_{\delta_{i}}}{z_{\delta_{i}d_{j}}}-\sum_{\delta_{r_{2}=1}^{\delta_{2}}} {\theta^{2}}{{}_{r_{2}}y^{2}}{{}_{r_{2}j}}\leq 0\] \[\sum_{r=1}^{z}{\theta_{r}}{y_{rj}}-\sum_{r_{2}=1}^{\delta_{2}}{ \theta^{2}}{{}_{r_{2}}y^{2}}{{}_{r_{2}j}}-\sum_{\delta_{i}=1}^{D_{i}}{\gamma_ {\delta_{i}}}{z_{\delta_{i}d_{j}}}\leq 0 \tag{13}\]  
Boundary constraints in Equation (12-13):  
\[{\gamma_{\delta_{1}}}{y_{\delta_{1}}}{z_{\delta_{1}}}{z_{\delta_{1}}}{z_{ \delta_{1}}}{z_{\delta_{1}}}{z_{\delta_{1}}}{z_{\delta_{1}}}{z_{\delta_{1}}}{z _{\delta_{1}}}{z_{\delta_{1}}}{z_{\delta_{1}}}{z_{\delta_{1}}}\]  
By simplifying Equation (6) based on Equation (11), we can also get the new \(T\)I (\(T\)I\({}_{h}\)):  
\[T_{h}=\sum_{i_{1}=1}^{m_{1}}{\pi_{i_{1}}}{x_{i_{1}k}}+\sum_{i_{2 }=1}^{m_{2}}{\pi_{i_{2}}}{x_{i_{2}k}}+\sum_{\delta_{r_{2}=1}^{\delta_{2}}}{ \theta^{2}}{{}_{r_{2}}y^{2}}{{}_{r_{2}k}}+\sum_{r_{2}=1}^{\delta_{2}}{\theta^{ 2}}{{}_{r_{2}}y^{2}}{{}_{r_{2}k}}+\sum_{r_{2}=1}^{\delta_{2}}{\theta^{2}}{{}_ {r_{2}}y^{2}}{{}_{r_{2}k}} \tag{14}\] \[+\sum_{\delta_{r_{2}=1}^{\delta_{2}}}{\gamma_{\delta_{r_{2}}}}{z_{ \delta_{2}d_{2}k}}+\sum_{\delta_{r_{2}=1}^{\delta_{2}}}{\gamma_{\delta_{r_{2}}}} {z_{\delta_{2}d_{2}k}}\]  
Once we obtain an optimal solution to linear model (12) and (13), the PSI for the three-stage network processes can be calculated accordingly. However, because of the linear systems Equations (12) and (13) has multiple unknown variables, so the linear model may have multiple optimal solutions, the respective \(E_{p}\) of the sub-DMUs may not be unique. Therefore, we can choose a set of multipliers based on three sub-DMUs that produce the highest first or second or third stage \(E_{p}\) while maintaining the PSI of the whole network processes. In addition, minimizing the efficiency gap between the three sub-DMUs in the three stages is also a standard strategy for choosing multipliers [27]. There is the environmental efficiency in stage 1, 2, and 3, respectively (\(E_{p}\) in stage 1, 2, 3): simplify by the Charnes-Cooper transformation, we also can get Equations (13) and (15) in Stage 1 sub-DMU:  
\[ \begin{gathered}{E_{1}}=\max{\sum_{\delta_{i}=1}^{D_{i}}{\gamma_{ \delta_{i}}}{z_{\delta_{i}k}}}\\ \sum_{i_{2}=1}^{m_{1}}{\pi_{i_{1}}}{x_{i_{1}k}}+\sum_{\tau_{1}=1}^ {\delta_{1}}{\theta^{2}}{{}_{r_{1}}y^{2}}{{}_{r_{1}k}}=1\end{gathered}  \tag{15}\]  
As above, we also can get Equations (13) and (16) in Stage 2 sub-DMU:  
\[ \begin{gathered}{E_{2}}=\max{\sum_{\delta_{i}=1}^{D_{i}}{\gamma_{ \delta_{i}}}{z_{\delta_{i}k}}}\\ \sum_{i_{2}=1}^{m_{2}}{\pi_{i_{2}}}{x_{i_{2}k}}+\sum_{\tau_{2}=1}^ {\delta_{2}}{\theta^{2}}{{}_{r_{2}}y^{2}}{{}_{r_{2}k}}+\sum_{\delta_{1}=1}^{D_{ 1}}{\gamma_{\delta_{1}}}{z_{\delta_{1}k}}{\gamma_{\delta_{1}}}{z_{\delta_{1}k}}{ \gamma_{\delta_{1}}}\end{gathered}