Recently, I've been working on some two-stage stochastic programming problems. Due to the presence of integer second-stage variables in the model, I've preferred to use the Progressive Hedging Algorithm (PHA). I've implemented my algorithms in GAMS and run them successfully. However, a common issue between all problems was that occasionally PHA underestimated the optimal solution of the original problem (obtained via solving the deterministic equivalent problem).
I've found out that on some problems the method is highly sensitive to the penalty term, $r$, in the objective function depicted below (i.e., the one penalizing the distance between the current solution from the previous nonanticipative solution). Please note that $\rho$ and $\hat{x}$ are parameters determined iteratively within the PHA. Also, $c^Tx_k$ and $q^T_ky_k$ are the first-stage and second-stage costs, respectively.
\begin{align}\min\qquad&\sum_{k=1}^{K} p_{k}\left[c^Tx_k+q^T_ky_k+\rho^{v,T}_{k}(x_k-\hat{x}_k)+\dfrac{r}{2}\lVert x_k-\hat{x}_k\rVert^2\right]\\\text{s.t.}\qquad&Ax_k = b,\\\qquad&Wy_k=h_k-T_kx_x,\\\qquad&x_k \ge 0,\,\exists x_k\in \mathbb{Z}^+,\\\qquad&y_k \ge 0, \,\exists y_k\in \mathbb{Z}^+.\end{align}
Unfortunately, tuning $r$ could take considerable time and be a pain. I was wondering if there is any way to tune $r$ and find the optimal or near-optimal solutions with less hassle. I would really appreciate it if you could share your experience and tricks to overcome this difficulty.
Update #1: The problems under study are completely linear. In other words, they are two-stage stochastic mixed-integer programming (with integer first- and second-stage variables).
Update #2: As my first-stage variables are binary, I'm able to completely linearize the term $\lVert x_k-\hat{x}_k\rVert^2$. Therefore, the problem solved in each iteration of PHA is completely linear. In addition, it is separable under index $k$ for different scenarios.
Update #3: The progressive hedging algorithm was implemented as discussed here on page 257.
Step 0. Suppose some nonanticipative $x_0$, some initial multiplier $\rho_0$, and $r > 0$ . Let $ν = 0$. Go to Step 1.
Step 1. Let $(x_k^{v+1},y_k^{v+1}) \; \forall k = 1,...,K$ solve the problem above. Let $\hat{x}^{v+1} = (\hat{x}^{v+1,1},...,\hat{x}^{v+1,K})^T$ where $\hat{x}^{v+1,k} = \sum_{l=1}^{K} p_lx^{v+1,l} \; \forall k = 1,...,K$.
Step 2. Let $\rho^{v+1}=\rho^{v} + r(x^{v+1,k}-\hat{x}^{v+1})$. If $\hat{x}^{v+1}=\hat{x}^{v}$ and $\rho^{v+1}=\rho^{v}$, then, stop; $\hat{x}^{v}$ and $\rho^{v}$ are optimal. Otherwise, let $v = v+1$ and go to Step 1.
