Improved Approximation Guarantees for Power Scheduling Problems With Sum-of-Squares Constraints

Management of power systems involves operational routines concerned with the optimal allocation/dispatch of limited resources (power/generation) under various types of objectives and constraints, including total fuel cost, total profit, emissions, active power loss, and voltage magnitude deviation, to name a few. These problems are conventionally modeled as quadratic programs (possibly with discrete variables), which are in general hard to solve exactly or approximately. In this paper, we focus on a class of scheduling problems where the energy constraints are modeled as a sum of squares of linear functions. Such constraints have been studied in a number of works, including electric vehicle charging in smart grids (Yu and Chau 2013; Karapetyan et al. 2018; Elbassioni, Karapetyan, and Nguyen 2019; Khonji et al. 2019); welfare maximization in gas supply networks (Klimm et al. 2022), and scheduling of tasks on machines (Wierman, Andrew, and Tang 2012; Bansal, Kimbrel, and Pruhs 2004).

In (Yu and Chau 2013), the authors investigated the so-called Complex-Demand Knapsack Problem (CKP), formulated to optimize the distribution of power among consumers within Alternating Current (AC) electrical systems. In CKP, power (representing user demand) is expressed in complex numbers (Grainger and Stevenson 1994), where the real part corresponds to active power, the imaginary part stands for reactive power, and the magnitude of the complex number measures the apparent power. Each customer possesses a specific power demand and receives a certain utility if their demand is fully met. However, due to physical limitations on the magnitude of the total power supply available in the system, only a subset of users can be selected for power allocation. Mathematically, the power constraint is formulated as a quadratic function, which is the sum of two squares of linear functions.

The sum-of-squares constraint has also appeared in the context of gas supply networks (Klimm et al. 2022). As noted in (Klimm et al. 2022), a gas pipeline network can be modeled as a directed path graph with different entry and exit nodes. Transportation requests are characterized by their entry and exit nodes and the amounts of gas to be transported. The gas flow on an edge from node $u$ to a node $v$ is modeled by the Weymouth equations (Weymouth 1912), describing the difference of the squared pressures between these nodes. For the operation of gas networks, one needs to choose a subset of transportation requests to be served such that the maximum difference of the squared pressures in the network does not exceed a given threshold. This pressure constraint can then be formulated as a sum of squares of linear functions.

Mathematically, the above two problems can be cast as a binary program that seeks to maximize a linear objective function subject to a sum-of-squares constraint. This problem covers the classical Knapsack problem as a special case and is thus $\mathrm{NP}$-hard. The first algorithmic result was given in (Woeginger 2000), where it was proved that the problem is strongly $\mathrm{NP}$-hard, and thus does not admit a fully polynomial time approximation scheme (FPTAS). Yu and Chau (Yu and Chau 2013) developed a $(0.5-\epsilon)$-approximation algorithm, which was subsequently improved to a PTAS in (Chau, Elbassioni, and Khonji 2016), the best possible result one can hope for the problem. In (Elbassioni and Nguyen 2017) the authors further extended this result to the case where the constraint contains a constant number of squares and devised a ($\frac{1}{e}-\epsilon$)-approximation algorithm for variant with a submodular objective. Recently, (Klimm et al. 2022) studied the version with an unbounded number of squares and presented a constant $0.618$-approximation algorithm for a linear objective and ruled out the existence of an approximation factor of $0.99$. These results rely on the so-called pipage rounding technique. One interesting question raised here is whether this rounding technique can be extended to nonlinear objectives such as quadratic and submodular functions, which are frequently encountered in real-world economic contexts.

The minimization version of CKP ($\mathsf{MinCKP}$) with a linear objective and a sum-of-two-squares constraint has been explored in (Elbassioni, Karapetyan, and Nguyen 2019). As noted therein, this version constitutes a special case of the Economic Dispatch problem (a.k.a, Unit Commitment, Generation Dispatch Control), which is principal to power systems engineering and is responsible for minimizing generation costs while maintaining the balance between the net supply and demand (Wood and Wollenberg 2012). Unlike the maximization version, however, $\mathsf{MinCKP}$’s natural relaxation is non-convex and hence harder to approximate. As of now, the best-known approximation was the quasi-PTAS (QPTAS) attained in (Elbassioni, Karapetyan, and Nguyen 2019) via a geometric approach, whereas the existence of a PTAS remained an open question.

Taking a step forward, we present two novel approximations for $\mathsf{MinCKP}$ and a generalized variant of CKP. More concretely, the contributions of this study are as follows:

• •

  First, we provide a PTAS for $\mathsf{MinCKP}$, which is a notable improvement upon the state-of-the-art result of a QPTAS provided in (Elbassioni, Karapetyan, and Nguyen 2019). At the heart of the approximation scheme is a polynomial time algorithm for solving the (non-convex) relaxation problem to optimality, by exploiting the special structure of optimal solutions. Note that a PTAS is the best possible result one can hope for, unless $\mathrm{P}=\mathrm{NP}$, as one can use a technique similar to that in (Woeginger 2000)¹¹1In (Woeginger 2000) the author provided a reduction from a variant of Partition to disprove the existence of an FPTAS for the CKP problem. With some minor modifications, this reduction can also be established for the minimization version of CKP. to rule out the non-existence of an FPTAS.

• •

  Next, we show that any $\alpha$-approximation algorithm for the relaxation problem maximizing a (non-monotone) submodular quadratic function over a general (non-negative) convex quadratic constraint can lead to an $(\frac{\alpha\phi^{2}}{2})$-approximation algorithm for the discrete case, where $\phi=\frac{2}{1+\sqrt{5}}$ is the reverse golden ratio. The idea is to generalize the technique of (Klimm et al. 2021) developed for a similar problem with a linear objective: given a fractional solution, round it in such a way that the resulting discrete solution remains feasible and the objective value does not decrease. Since in the studied problem both the objective and constraint are quadratic, we have to develop a more sophisticated rounding technique than that derived in (Klimm et al. 2022) for the linear objective case.

In this section, we formalize the mathematical models of the two power scheduling problems under study. The first one, termed $\mathsf{MinCKP}$, is the minimization version of CKP (Yu and Chau 2013) and is defined as

where ${\bf p},{\bf q}\in\mathbb{Q}_{+}^{n}$, ${\bf c}\in\mathbb{Q}^{n}$ and $C>0$. We may assume, w.l.o.g., that ${\bf c}>0$ (if $c_{i}\leq 0$ for some $i\in[n]$, the corresponding variable $x_{i}$ can be set to $1$). Also, it is assumed that ${\bf p},{\bf q}$ are non-negative, otherwise the problem doesn’t admit an approximation algorithm with finite ratio (Chau, Elbassioni, and Khonji 2016). In $\mathsf{MinCKP}$, the scalar $C$ denotes the aggregate apparent power demand to be met, while the column vectors ${\bf p}^{\top}$ and ${\bf q}^{\top}$ capture the active and reactive power outputs of the contracted generation units. Accordingly, the vector ${\bf c}$ represents the costs associated with dispatching these units. These costs may quantify the expenses carried by the load-serving entity (operator of the grid), including the utilization of a generation source, or alternatively, the amount paid due to CO2 emissions.

The second problem, referred to as $\mathsf{MaxSUB}$, generalizes the one studied by (Elbassioni and Nguyen 2017; Klimm et al. 2022) and takes the following form:

where $K\leq n$, ${\bf p}_{k}\in\mathbb{Q}_{+}^{n}$ and the objective is submodular, i.e., the matrix $A$ has non-positive off-diagonal entries. Let ${\bf a}$ be the vector of diagonal entries of $A$. One can assume, w.l.o.g., that ${\bf a}$ is non-negative, since if $a_{i}<0$ then setting $x_{i}=0$ does not affect the optimal solution. The constraint (2) can be written in the form of ${\bf x}^{T}Q{\bf x}\leq C$, for some positive semi-definite (PSD) matrix $Q\in\mathbb{Q}_{+}^{n\times n}$. Since ${\bf x}\in\{0,1\}^{n}$, one can write the objective function in the equivalent form of $f({\bf x}):={\bf x}^{\top}A^{*}{\bf x}+{\bf a}^{\top}{\bf x}$, where $A^{*}$ is obtained from $A$ by setting all the diagonal entries to $0$, without changing optimal solutions. It is not hard to verify that in this case the (discrete) function $f({\bf x})$ (over $\{0,1\}^{n}$) satisfies submodularity, i.e., $f({\bf x}^{S})+f({\bf x}^{V})\geq f({\bf x}^{S\cup V})+f({\bf x}^{S\cup V})$, for every two sets $S,V\subseteq\{1,2,\ldots,n\}$, where ${\bf x}^{S}$ denotes the binary vector with $x_{i}^{S}=1$ if and only if $i\in S$.

In this section, we present a PTAS for $\mathsf{MinCKP}$. As previously remarked, the relaxation with continuous variables in $[0,1]$ is non-convex and thus cannot be solved by existing techniques for convex (quadratic) programs. Interestingly, by exploiting the structure of a class of the relaxation’s optimal solutions, and carefully analyzing an equivalent dual program, we can reduce the space of candidates so that an optimal solution can be found efficiently. It can then be used to derive an integer solution with a bounded guarantee.

Denote by $\mathsf{NLP}$ the (non-convex) relaxation of $\mathsf{MinCKP}$. It might be the case that the optimum solution ${\bf x}$ for $\mathsf{NLP}$ is attained at the boundary of the feasible set, and thus has irrational components. Hence, by saying “an optimal solution” we essentially mean a rational solution $\widehat{}{\bf x}\in[0,1]^{n}$ with ${\bf c}^{\top}\widehat{}{\bf x}\leq{\bf c}^{\top}{\bf x}^{*}+\delta$ and $({\bf p}^{\top}\widehat{}{\bf x})^{2}+({\bf q}^{\top}\widehat{}{\bf x})^{2}\geq C-\delta$, for an arbitrary small $\delta>0$, where ${\bf x}^{*}$ is an optimum solution for $\mathsf{NLP}$. For a given $\delta>0$, we say that $\widehat{}{\bf x}$ is $\delta$-feasible for $\mathsf{NLP}$ if it satisfies the second inequality, and $\delta$-optimal if it satisfies both inequalities. At the heart of our method is the following key result.

We assume that $\|{\bf p}\|_{1}^{2}+\|{\bf q}\|_{1}^{2}\geq C$, so that the problem is feasible. We also assume that ${\bf p}$ and ${\bf q}$ are linearly independent vectors; otherwise, the problem simplifies to a linear program. We will prove that, given a fixed positive value $T$, the problem of finding a $\delta$-feasible solution, if it exists, of value exactly $T$ can be achieved in polynomial time. We denote this problem by $\mathsf{NLP}[T]$. As the objective function and the constraint (2) are both monotone²²2Note that, given a feasible solution to $\mathsf{NLP}$ with value ${\bf c}^{\top}{\bf x}<T$, we can increase the components of ${\bf x}$, one by one, without violating the inequality (2), until either ${\bf c}^{\top}{\bf x}$ becomes equal to $T$, or $x_{i}$ becomes equal to $1$ for all $i$; the latter case cannot happen unless $T\geq\|{\bf c}\|_{1}$. Thus, checking if the optimum value of $\mathsf{NLP}$ is at most $T$ is equivalent to checking if $\mathsf{NLP}[T]$ is feasible., this result combined with a binary search (via Algorithm 1) yields an efficient algorithm for finding a $\delta$-optimal solution to $\mathsf{NLP}$, and hence proving Theorem 2.

We first give a structural property of an optimal solution of $\mathsf{NLP}$, which can help to check the feasibility of $\mathsf{NLP}[T]$ efficiently. Let ${\bf x}^{*}$ be an optimal solution to $\mathsf{NLP}$. We will use ${\bf x}^{*}$ as a parameter in analyzing the structure of an optimal solution. Let $\alpha={\bf p}^{\top}{\bf x}^{*}$ and $\beta={\bf q}^{\top}{\bf x}^{*}$. Consider the following linear program:

Note that ${\mathsf{LP}}$ and $\mathsf{NLP}$ have the same optimal solution. Indeed, it can be seen that every feasible solution to ${\mathsf{LP}}$ is also feasible to $\mathsf{NLP}$. Since both programs have the same objective, every optimal solution to ${\mathsf{LP}}$ must also be optimal to $\mathsf{NLP}$. In what follows we make use of the dual of ${\mathsf{LP}}$ in analyzing the structure of such a solution.

The dual of (${\mathsf{LP}}$) is

Since the primal is feasible with the primal solution ${\bf x}^{*}$, then the dual is also feasible and let $({\bf z}^{*},{\bf y}^{*})$ be an optimal (basic feasible) dual solution³³3Note that the linear independence of ${\bf p}$ and ${\bf q}$ guarantees that the optimum of the dual is attained at a vertex. corresponding to the primal solution ${\bf x}^{*}$. The optimum values of the primal-dual pair coincide by the strong duality criterion. Moreover, by the complementary slackness condition, it holds that, for all $i\in[n]$,

(i)

$x^{*}_{i}=0~{}\text{or}~{}p_{i}\cdot z^{*}_{1}+q_{i}\cdot z^{*}_{2}-y^{*}_{i}-c_{i}=0$,

(ii)

$x^{*}_{i}=1~{}\text{or}~{}y^{*}_{i}=0.$

The complementary slackness criterion implies that the primal solution ${\bf x}^{*}$ satisfies the following conditions

Note that $x^{*}_{i}$ can take any value in $[0,1]$ when the equality $p_{i}\cdot z^{*}_{1}+q_{i}\cdot z^{*}_{2}=c_{i}$ happens.

The following observation is due to the linear independence of ${\bf p}$ and ${\bf q}$ and the fact that $({\bf z}^{*},{\bf y}^{*})$ is assumed to be a basic feasible solution to the dual.

Indeed, since $({\bf z}^{*},{\bf y}^{*})$ is a basic feasible solution to the dual, there must exist $n+2$ linearly independent constraints among the constraints (3)-(4) that are tight at $({\bf z}^{*},{\bf y}^{*})$. This implies that we must have two distinct indices $j,k\in[n]$ such that the two inequalities (3) and (4) hold as equalities for $i=j,k$ and such that the four obtained equations are linearly independent. This leads to the claim in the observation.

The above observation implies that $z^{*}_{1},z^{*}_{2}$ are uniquely and fully determined by a $2\times 2$ linear system obtained by selecting two indices $j,k\in[n]$ such that the two vectors $(p_{j},p_{k})$ and $(q_{j},q_{k})$ are linearly independent. Once the solutions $z_{1}^{*},z_{2}^{*}$ are obtained, one can classify the values of variables in the primal solution ${\bf x}$ as follows

• •

  If $p_{i}\cdot z^{*}_{1}+q_{i}\cdot z^{*}_{2}<c_{i}$, then $x_{i}=0$,

• •

  If $p_{i}\cdot z^{*}_{1}+q_{i}\cdot z^{*}_{2}>c_{i}$, then $x_{i}=1$,

• •

  If $p_{i}\cdot z^{*}_{1}+q_{i}\cdot z^{*}_{2}=c_{i}$, then $x_{i}\in[0,1]$.

In other words, based on the values of $z^{*}_{1},z^{*}_{2}$, we know that some of the variables in the primal solution must be equal to $1$ or $0$, but possibly not all. Define $\Delta_{1},\Delta_{0}$ to be the sets containing indices $i$ with $x^{*}_{i}=1$ and $x^{*}_{i}=0$, respectively, and let $\Delta=[n]\setminus\Delta_{1}\cup\Delta_{0}$. Then the values of $x^{*}_{i}$’s, $i\in\Delta$, must be the optimal to the $\mathsf{NLP}$ program ($\mathsf{NLP}[\Delta]$), parameterized by $\Delta$:

where ${\bf p}(\Delta_{1})={\sum}_{i\in\Delta_{1}}p_{i},{\bf q}(\Delta_{1})={\sum}_{i\in\Delta_{1}}q_{i},{\bf c}(\Delta_{1})={\sum}_{i\in\Delta_{1}}c_{i}$. It follows that the problem of checking if $\mathsf{NLP}[T]$ has a feasible solution is equivalent to checking if there is a feasible solution to $\mathsf{NLP}[\Delta]$ with $\sum_{i\in\Delta}c_{i}x_{i}+{\bf c}(\Delta_{1})=T$. By the definition of $\Delta$, it follows that

Based on this, the problem is now to check if the following nonlinear system is feasible with variables $\xi=\sum_{i\in\Delta}p_{i}x_{i}\in[0,\overline{\xi}],~{}\zeta=\sum_{i\in\Delta}q_{i}x_{i}\in[0,\overline{\zeta}]$:

where $\overline{\xi}={\bf p}(\Delta)$, and $\overline{\zeta}={\bf q}(\Delta)$ are upper bounds on the values of $\xi$ and $\zeta$, respectively.

Since we assume ${\bf c}>0$, it follows that at least one of the dual solutions $z_{1}^{*},z_{2}^{*}$ must be positive. Without loss of generality, we assume $z_{2}^{*}>0$. To check the feasibility of (5), one can compute $\zeta$ as a function of $\xi$ from the first equation, and then plug it into the second one, leading to a quadratic inequality in one variable $\xi$:

from which one can easily get the solution of $\xi$ as

where $\xi_{1},\xi_{2}$ are the two (possibly identical) roots of the quadratic equation obtained by setting the inequality in (6) to an equality⁴⁴4Note that it is possible that the quadratic equation has no real solution in which case, we set $\xi_{1}=-\infty$ and $\xi_{2}=+\infty$; this means that our assumption about the feasibility of $\mathsf{NLP}[T]$ is wrong.. Note that, depending on the signs of $\xi_{1}$ and $\xi_{2}$, $\Omega(\xi)$ consists of at most two intervals on the real line. Using each interval in $\Omega(\xi)$ and the equation in (5), the remaining variables $x^{*}_{i}$’s, $i\in\Delta$, can be recovered by solving at most 2 linear programs, one for each choice of $I\in\Omega(\xi)$:

Note that, due to the possible non-rationality of $\xi_{1}$ and $\xi_{2}$, we may need to work with $\delta^{\prime}$-approximations of $\xi_{1}$ and $\xi_{2}$, i.e., compute numbers $\tilde{\xi}_{1}$ and $\tilde{\xi}_{2}$ such that $|\xi_{1}-\tilde{\xi}_{1}|\leq\delta^{\prime}$ and $|\xi_{2}-\tilde{\xi}_{2}|\leq\delta^{\prime}$. Such approximations can be found in time $\mathsf{poly}(L,\log\frac{1}{\delta^{\prime}})$, and by choosing $\delta^{\prime}$ sufficiently small, we can guarantee that (6) is satisfied within an additive error of $\delta$, thus yielding a $\delta$-feasible solution for $\mathsf{NLP}[T]$.   ❑ 

Let $z^{*}$ and $z^{\prime}$ be the optimal values for $\mathsf{NLP}[S_{0}^{*},S_{1}^{*}]$ and ${\mathsf{LP}}^{\prime}[S_{0}^{*},S_{1}^{*}]$, respectively. Then, as $x^{*}$ is both $\delta$-optimal for $\mathsf{NLP}[S_{0}^{*},S_{1}^{*}]$ and feasible for ${\mathsf{LP}}^{\prime}[S_{0}^{*},S_{1}^{*}]$, we must have $z^{\prime}\leq z^{*}+\delta$. Hence, there exists a basic optimal solution of ${\mathsf{LP}}^{\prime}[S_{0}^{*},S_{1}^{*}]$, with at most $2$ fractional components (which can be found efficiently), which is also a $\delta$-optimal solution of $\mathsf{NLP}[S_{0}^{*},S_{1}^{*}]$. By rounding each of these two fractional components to $1$, we obtain a $\delta$-feasible solution $\widehat{}{\bf x}$ to $\mathsf{MinCKP}$ of value ${\bf c}^{\top}\widehat{}{\bf x}$ at most

if we choose $\delta=\epsilon\cdot 2^{-2L-1}\leq\epsilon\cdot\min_{j\in[n]}c_{j}$. Observe that, for this choice of $\delta$, $({\bf p}^{\top}\widehat{}{\bf x})^{2}+({\bf q}^{\top}\widehat{}{\bf x})^{2}\geq C-\delta$ implies that $({\bf p}^{\top}\widehat{}{\bf x})^{2}+({\bf q}^{\top}\widehat{}{\bf x})^{2}\geq C$, since $C-(({\bf p}^{\top}\widehat{}{\bf x})^{2}+({\bf q}^{\top}\widehat{}{\bf x})^{2})$ is a rational number of value at least $2^{-2L}>\delta$. We conclude that $\widehat{}{\bf x}$ is feasible solution to $\mathsf{MinCKP}$ of value at most $(1+O(\epsilon))\mathrm{OPT}$. The running time is $n^{\frac{2}{\epsilon}+O(1)}L^{O(1)}$.   ❑ 

In this section we provide a pipage rounding procedure for MaxSUB, given a (fractional) solution of its relaxation (denoted as $\mathcal{F}$, obtained by changing binary variables to continuous ones). Our aim is to prove the following result.

To prove the Theorem 3, we construct an algorithm, which is formally presented in Algorithm 3, and prove its correctness. Before doing so, (inspired by (Klimm et al. 2022)) let us introduce two relaxations of $\mathsf{MaxSUB}$, which are needed in the performance analysis of the algorithm later on. In both relaxations, we replace the binary variables ${\bf x}\in\{0,1\}^{n}$ by ${\bf x}\in[0,1]^{n}$. Recall that $f({\bf x})={\bf x}^{\top}A^{*}{\bf x}+{\bf a}^{\top}{\bf x}$, define

where ${\bf q}$ is the vector containing all the diagonal entries of $Q$. It is worth noting that, because of the non-negativity of $Q$, $\mathcal{F}$ and $\mathcal{F}_{1}$ are equivalent.