Cutting-planes application procedure for a specific problem

Question

Sort of following up with this question. I reformulated another model to make it convex and possibly solve it with some cut generation method. I would like to double-check whether I am doing it correctly. Below is the nonlinear model due to $1/h_p$. The variable $h_p$ is non-negative continuous and satisfies $H^- \leq h_p \leq H^+$. The parameters: $f_p,B_p,\tau_p,M,H^-,H^+$ are non-negative real numbers. The set $P$ is polynomially-sized.

\begin{alignat}2\min &\quad \mathbf{C} = \sum_{p\in P}\frac{h_p}{2}\tag1\\\text{s.t.}&\quad \frac{f_ph_p}{30}\leq B_p \qquad \forall p\in P\tag2\\&\quad\sum_{p\in P}\frac{\tau_p}{h_p} \leq M\tag3\\&\quad h_p\in \mathbb{R}^+, H^- \leq h_p \leq H^+.\end{alignat}

Theorem 1: Assume $\phi_p\left(h_p\right)=\frac{\tau_p}{h_p}$. Then, $\phi_p\left(h_p\right)$ is convex in $h_p$ under the domains $h_p,\tau_p\in\mathbb{R}^+$.

Proof: Showing the second derivative of $\phi_p\left(h_p\right)$ with respect to $h_p$ being non-negative will prove the convexity. Since $\frac{d^2 \phi_p\left(h_p\right)}{dh_p^2}=\frac{2\tau_p}{h_p^3}\geq 0$ in the domains $h_p,\tau_p\in\mathbb{R}^+$, my proof is complete.

I will introduce $(4)$ to represent the new definition assuming $\phi_p\equiv\phi_p\left(h_p\right)$. Due to Theorem 1, I will say $a_p+b_p h_p$ supports $\phi_p\left(h_p\right)$ at $h_p=\tilde{h}_p$, where $a_p=\phi_p\left(\tilde{h}_p\right)-b\tilde{h}_p$ and $b_p=\frac{d\phi_p\left(\tilde{h}_p\right)}{d\tilde{h}_p}$. So, if I introduce the cut constraint $(5)$ to the problem $(1)-(2), ~(4)$ iteratively, I expect to solution to converge to optimality.

\begin{alignat}2\sum_{p\in P}\phi_p&\leq M\,\tag4\\\phi_p &\geq a_p+b_p h_p.\tag5\end{alignat}

Here is my planned solution procedure. Solve $(1)-(2),~(4)$ with non-negativity conditions with $h_p$ boundary and $\phi_p\geq 0 $. The solution is $h_p=H^-~\forall p\in P$ due to minimization assuming $30B_p/f_p\geq H^-$. If $(3)$ is not satisfied with this solution, feed the solution $H^-$ into $\tilde{h}_p$, introduce $(5)$, resolve the problem...

I am sort of stuck with the solution procedure as I cannot really determine the termination criterion. I know I find a lower bound to $\mathbf{C}$ by iteratively solving. But, how can I calculate the upper bound (call it $\mathbf{\hat{C}}$)? If I could understand it, I would say while the gap between the bounds are less than a satisfactory ratio, keep adding cuts; terminate when gap satisfies the ratio.

Is there any better approach to attack solving this problem?

Answer 1

If $|P|$ is not too large, you could try an integer programming formulation. Fix an integer $N>1$ (which will control the granularity of the approximation) and let $\Delta=\frac{H^+ - H^-}{N}$. For each $p\in P$ and each $n\in \lbrace 0,\dots, N\rbrace$, introduce variable $t_{p,n}\in [0,1]$. Now add the constraints$$\sum_{n=0}^N t_{p,n} = 1\quad \forall p$$and $$h_p = H^- + \sum_{n=0}^N t_{p,n} \cdot (n\Delta).$$Also force $\lbrace t_{p,0},\dots, t_{p,N}\rbrace$ to be a type 2 special ordered set (SOS2), which is what makes the problem an integer program. You will need an integer programming solver that understands the SOS2 designation (or else you will need to enforce it with an alternative formulation, which is possible). In essence, what this does is discretize the interval $[H^-,H^+]$ and express each $h_p$ as a convex combination of two adjacent grid points in the interval.

Now we address your constraint (3). Replace the left side of (3) with $$\sum_{p\in P}\sum_{n=0}^N \left(\frac{\tau_p}{H^- + n\Delta}\right) t_{p,n}.$$ This interpolates each $\phi_p()$ linearly, so it's an approximation.

Assuming your solver can solve the approximate model, calculate the discrepancy if any in (3). If you can live with it, great. If not, you can try tightening the approximation of each $h_p$ in the vicinity of its "optimal" value. I posed the initial approximation as a uniform grid for simplicity, but in fact you can space the grid points any way you like. So to refine your solution, you can increase the density of grid points near the values the solver chose for the $h_p$ and optionally reduce the density elsewhere to conserve model size.