## 2 The TPUMSTSR Under RLV Model
### Solution Method/Algorithm  
The solution process for this model is illustrated in Figure 5.  
**(1)**: Stage One: Solving the Optimal Power Flow Model of the Main Grid  
Each grid-connected CS and independently dispatched power source calculates its own upper and lower power output limits for the current time period. Those limits are then transferred to the main grid. After receiving the feasible power region constraints from all participants, the main grid dispatch center solves for the minimum required power generation for each grid-connected CS and IPP based on the predicted load demand. This information is then provided to the corresponding CSs and IPPs. **(2)**: Stage Two: Solving the Total PR Revenue \(C_{t}^{\text{TN,D}}\)  
Since CA is essentially the distribution of the compensation revenue, in terms of the time sequence, it is possible to first compute the peak-shaving compensation revenue and then calculate the CA for the units that do not participate in peak shaving. Therefore, in Stage Two, the internal operating model of the MCS without CA is solved, and the computed \(C_{t}^{\text{TN,D}}\) is fixed. In other words, \(C_{t}^{\text{TN,D}}\) can initially be approximated as a constant representing the PR revenue and then incorporated into Stage Three for further optimization. Figure 5: Computation flowchart for economic operation of three-stage MCS considering CA. After solving Stage Two, the nonlinear objective involving cross-calculations of fractional and multiplicative variables in the original CA is partially linearized, significantly reducing the solution difficulty and computation time. **(3)**: Stage Three: Solving the Economic Dispatch Model of CSs Considering CA  
After introducing CA, the PR compensation revenue \(C_{t}^{\text{TN,D}}\) is fixed as a constant in Stage Two. Therefore, the CA calculation only needs to address the fractional nonlinearity of the ratio of the adjusted output of a single unit to the total adjusted output of all units in the system. To handle the fractional nonlinearity in the objective function, this paper employs second-order cone programming (SOCP) for relaxation. Since the fractional nonlinearity in the CA calculation for thermal power and renewable energy units only differs in the numerator (adjusted output), the following example uses the CA for thermal power units not participating in PR to demonstrate SOCP relaxation:  
\[f_{i,u,t}^{\text{US}}=C_{t}^{\text{TN,D}}\frac{P_{i,u,t}^{\text{TPU,C}}}{ \sum\limits_{g\in G}P_{g,t}^{\text{C}}} \tag{33}\]  
Let \(P_{i,u,t}^{\text{TPU,C}}=A\sum\limits_{g\in G}P_{g,t}^{\text{C}}\), \(A^{2}=\alpha\), and \((\sum\limits_{g\in G}P_{g,t}^{\text{C}})^{2}=\beta\). The original nonlinear CA problem can be described as follows:  
\[f_{i,u,t}^{\text{US}}=C_{t}^{\text{TN,D}}A \tag{34}\]  
Additionally, the following second-order cone and variable bound constraints are introduced:  
\[\small \begin{array}{l}\small\begin{array}{l}\small\begin{array}{l}(P_{ i,u,t}^{\text{TPU,C}})^{2}\leq\alpha\beta\\ A^{2}\leq\alpha\\ (\sum\limits_{g\in G}P_{g,t}^{\text{C}})^{2}\leq\beta\\ \sum\limits_{g\in G}P_{g,t}^{\text{C}}\leq\sum\limits_{g\in G}P_{g,t}^{\text{C}} \leq\overline{\sum\limits_{g\in G}P_{g,t}^{\text{C}}}\\ P_{i,u,t}^{\text{TPU,C}}\leq P_{i,u,t}^{\text{TPU,C}}\leq P_{i,u,t}^{\text{TPU,C}}\leq P_{i,u,t}^{\text{TPU,C}}\end{array}\end{array}\end{array}  \tag{35}\]  
After relaxation, the linear CA makes the solution of the optimal dispatch scheme for MCSs convenient. To minimize the error introduced by \(C_{t}^{\text{TN,D}}\) in Stage Two, the PR compensation revenue \(f_{i,u,t}^{\text{VD,III}}\) from Stage Three is summed to obtain the total compensation revenue \(C_{t}^{\text{TN,D,III}}\), which is then compared with the fixed \(C_{t}^{\text{TN,D}}\) from Stage Two. If \(C_{t}^{\text{TN,D,III}}\geq C_{t}^{\text{TN,D}}\), the constant \(C_{t}^{\text{TN,D}}\) used in Stage Three is updated to \(C_{t}^{\text{TN,D,III}}\), and the Stage Three optimization model is solved again until \(C_{t}^{\text{TN,D,III}}=C_{t}^{\text{TN,D}}\).