Note that the constant term of \(R(1 - \sum _{k=1}^{K} \lambda _k)\)  is eliminated because it does not affect convexity. In problem (REF ), constraints (REF ) and () are coupling constraints, while (), (), and () are coupling variables. To handle this problem, we utilize the decomposition method to relax constraints (REF ), (), and () [1]}. Hence, the partial Lagrangian function of problem (REF ) can be expressed as
\(\begin{split}\mathcal {L}^{^{\prime }}(\mathbf {\tau },\mathbf {\varphi },\mathbf {e},\mathbf {\mu },\mathbf {\xi },\nu ) & = \sum _{k=1}^K \lambda _k R_k(\tau _k,\varphi _k,e_k) \\& - \sum _{k=1}^{K}{\mu _k \bigg ( e_k - \zeta P_t h_k \sum _{n=0}^{k-1}\tau _n \bigg )} \\& - \sum _{k=1}^K \xi _k(\varphi _k - \tau _k) - \nu \bigg ( \sum _{k=0}^{K}\tau _k - T \bigg ),\\\end{split}\) 
