# Continuous Time Multi-stage Stochastic Reserve and Unit Commitment

Kári Hreinsson, Bita Analui and Anna Scaglione  ECEE, Arizona State University
<EMAIL_ADDRESS>
###### Abstract

In this paper we introduce a continuous time multi stage stochastic
optimization for scheduling generating units, their commitment, reserve
capacities and their continuous time generation profiles in the day-ahead
wholesale electricity market. Our formulation approximates the solution of a
variational problem, in which the balance, generation capacity and ramping
constraints are in continuous time. Due to the greater accuracy of our
representation of ramping events this approach improves the system reliability
and lowers the real-time cost.

###### Index Terms:

Continuous time unit commitment, multistage stochastic programming, stochastic
unit commitment.

This work was funded in part by the Advanced Research Projects Agency-Energy
(ARPA-E), U.S. Department of Energy, under Award Number DWS1103. The views and
opinions of authors expressed herein do not necessarily state or reflect those
of the United States Government or any agency thereof.

##  I Introduction and Motivation

There is a vast literature dealing with stochastic versions of the unit
commitment (UC) and of the security constrained unit commitment (SCUC)
problems. Currently, uncertainty in power systems is managed by scheduling
reserve capacity in advance to compensate for errors in (net-)load forecasts.

Stochastic Unit Commitment (SUC) appears to be the appropriate tool to capture
the exogenous uncertainties directly into the decision process. In fact, SUC
ensures that a feasible solution exists for all considered scenarios, and that
the expected cost is minimized. The most commonly studied versions of SUC
formulations are two-stage SUC problems, pioneered by [1]. The multi-stage
stochastic unit commitment (MSUC) is an extension of the two-stage SUC where
decisions are taken sequentially at certain time instants. While the MSUC
formulation allows for smoother boundary conditions on the commitment
variables, it has received less attention than the two-stage SUC because the
complexity is often prohibitive. First in [2] and then in a series of follow
up work (see e.g. [3], [4], [5], [6] and the references therein) many authors
worked on curbing the MSUC computational complexity. It is natural to use the
more accurate representation of the uncertainties in SUC to optimally oversee
a more economic commitment of reserves. The authors of [7] have proposed a
scenario based approach for a short-term scheduling problem in the presence of
wind power generation. More specifically, their proposed approach focuses on
the lack of accurate prediction for wind power generation and the allocation
of reserve capacities from Hydro-Thermal generation units to secure the system
against immense variability of wind power generation. In a similar stream, our
work provides a novel solution to overcome the existing challenges in
committing units and scheduling reserves, but our focus is also on improving
the accuracy beyond the conventional stochastic optimization frameworks cited
above.

We observe that, while representing the decision uncertainty more accurately
leads to greater reliability, one of the aspects that the works cited above
ignore is the inadequacy in representing the continuous time scenarios and all
of the literature cited above fail in capturing the continuous time nature of
some actions, such as ramping needs and inter-hourly variations. For example,
in locations where solar power is significant, early morning and evening hours
have pronounced curvatures, due to the mismatch between demand peaks and the
solar production.

The main contribution of this paper is to show a viable numerical
approximation of the stochastic variational problem that power system
operators try to address in their day-ahead commitment and reserve decisions.
We leverage our prior work in [8] focused on the deterministic continuous time
unit commitment. More specifically, we introduce the continuous time multi-
stage reserve and unit commitment (CT-MSRUC) problem, in which an underlying
load scenario tree is used to decide the baseline day-ahead dispatch,
commitment and reserve capacity, considering continuous time generation
trajectories as part of the decision variables. Our approach allows to better
schedule the reserve capacity and power, not only because the formulation
accounts for the future expected real-time cost in dispatching them, but also
because it provides a more accurate representation of the possible future
inter-hourly ramping needs.

The paper is organized as follows. In Section II we introduce our
approximation of continuous-time trajectories and scenario trees of net-load.
The formulation of the CT-MSRUC follows in Section III. The formulation is
further shown to reduce to a framework similar to the conventional discrete
time representation of the generation trajectory, by relaxing specific
constraints associated with the time-continuity and using a first order
polynomial. This relaxation (which we refer to as the discrete-time MSRUC —
DT-MSRUC) is used as a benchmark to assess the performance of the proposed
framework in Section IV, where we show that the CT-MSRUC leads to greater
system reliability and lowers operating cost.

##  II General Polynomial Interpolation of Electric Net-Load Process on a
Scenario Tree

To approximate continuous net-load trajectories, we utilize the Bernstein
polynomial basis:

$b_{k,n}(x)=\binom{n}{k}x^{k}(1-x)^{n-k}\Pi(x),\quad k\in[0,n]$ (1)

where $\Pi(x)$ is the rectangular function and $\Pi(x)=1$ for $0\leq x\leq 1$.
We define the vector of polynomials of degree $n$ as $\mathbf{b}_{n}(x)$. The
Bernstein polynomial operator $\mathfrak{B}_{n}$ takes a function
$f(x)\in[0,1]$ and maps it into an $n$th order polynomial as
$\mathfrak{B}_{n}f(x)=\textstyle\sum_{k=0}^{n}f_{k}b_{k,n}(x)$, an
approximation of $f(x)$ (if Lipschitz continuous) with maximum error $O(1/n)$
when $f_{k}=f(k/n)$. The coefficients $f_{k}$ are called control points.

The uncertainty considered in this paper comes from net-load, which we
consider a continuous time, continuous space stochastic process $\Xi(t)$. To
make the problem computationally tractable, the $\sigma$-algebra of this
process must be approximated by a discrete and finite state space stochastic
process, whose probability model is truncated over the finite horizon of the
optimization. A scenario tree
$\mathcal{T}=\\{\mathcal{V},\mathcal{E},\mathbb{P};\boldsymbol{\xi}\\}$
represents the basic structure that is used in multi-stage stochastic
optimization to represent the gradual unfolding of information through
filtration $\mathfrak{F}$ of the underlying stochastic process. The tree is a
directed graph from the root node $v_{0}$, (node/vertex in graph theory sense)
of the tree to its leaf nodes, through edges $(v_{-},v)\in\mathcal{E}$, where
$v_{-}$ is the parent of node $v$, and $v,v_{-}\in\mathcal{V}$. The scenario
tree is constructed to meet some minimum approximation error, relative to a
metric that is typically a surrogate for minimizing the optimality loss due to
the inaccuracy of the simplified decision model. In the conventional
construction of a tree (cf. [9]), each node of the tree has an associated
discrete time net-load sample, and the representation of the statistics is
nodal, that is: each node $v\in\mathcal{V}$ of the tree has a certain
probability $\pi_{v}$ that is the marginal probability of that particular
sample for the quantized process at that stage. A better interpretation of the
approximate scenario tree for a continuous time process is that realizations
are mapped into piece-wise constant sample paths.

We extend this notion by applying the Bernstein operator of order $n$ to each
interval mapping the sample net-load paths into a process whose realizations
are random piece-wise polynomial functions of order $n\geq 0$. For the
resulting process, we can construct the corresponding filtration
$\mathfrak{F}$ by representing the sub-$\sigma$ algebras of the set of random
coefficients. More specifically, for all nodes $v\in\mathcal{V}$, we define
the set of its direct successor nodes by $v_{+}\in\mathcal{C}(v)$ and assume
that $t_{v_{-}}\leq t<t_{v}$ represents the continuous time between any two
consecutive nodes $v_{-}\leavevmode\nobreak\ \text{and}\leavevmode\nobreak\
v$. We use $\pi_{v}=\pi_{v_{-},v}$ to represent the joint probability of the
free parameters in the polynomial trajectory and use the notation
$\boldsymbol{\xi}\sim\mathcal{T}$ to refer to the complete construction. We
divide the time between consecutive nodes $[t_{v_{-}},t_{v}]$ into $n+1$
arbitrary subintervals and define a vector of control points:

${\boldsymbol{\xi}}_{v_{-},v}=\left[\xi^{(0)}_{v_{-},v},\ldots,\xi^{(n-1)}_{v_{-},v},\xi^{(n)}_{v_{-},v}\right]^{\mathrm{T}}$
(2)

to represent the coefficients of net-load in the function space of polynomials
of degree $n$. The continuous net-load scenario approximation in each edge’s
time interval is:

$\xi_{v_{-},v}(t)=\mathbf{b}_{n}\left(\frac{t-t_{v_{-}}}{t_{v}-t_{v_{-}}}\right){\boldsymbol{\xi}}_{v_{-},v},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ t_{v_{-}}\leq t<t_{v}$ (3)

where $\left(\frac{t-t_{v_{-}}}{t_{v}-t_{v_{-}}}\right)$ scales the continuous
time index to fit into the time interval corresponding to the edge
$(v_{-},v)$. This approximation has several useful properties, one of which is
that $\xi_{v_{-},v}(t_{v_{-}})=\xi_{v_{-},v}^{(0)}$ and
$\xi_{v_{-},v}(t_{v})=\xi_{v_{-},v}^{(n)}$. In order to maintain $C^{0}$
continuity across edges of the load tree, it is sufficient to enforce that the
control points match at the endpoints. Another useful property of Bernstein
polynomials is that the coefficients of an $n^{\text{th}}$ degree polynomial’s
derivative are the finite differences of the original coefficients:
$\dot{b}_{k,n}(t)=n(b_{k-1,n-1}(t)-b_{k,n-1}(t))$. Therefore, a linear matrix
${\bf M}_{n+1\times n}$, which is rectangular-bidiagonal, can be introduced to
translate the derivatives of ${\bf b}_{n}(t)$ into the same family of
polynomials of degree $n-1$:

$\displaystyle\dot{\bf b}_{n}(t)=n{\bf M}{\bf b}_{n-1}(t),$ (4)

Consequently, the elements of $\dot{\boldsymbol{\xi}}_{v_{-},v}$ can be
expressed as the linear combination of the control points of
${\boldsymbol{\xi}}_{v_{-},v}$ as ${\dot{\boldsymbol{\xi}}_{v_{-},v}=n{\bf
M}{\boldsymbol{\xi}}_{v_{-},v}}$. Specifically, at the end points:

$\displaystyle\dot{\boldsymbol{\xi}}_{v_{-},v}=\left[n\left(\xi^{(1)}_{v_{-},v}-\xi^{(0)}_{v_{-},v}\right),\ldots,n\left(\xi^{(n)}_{v_{-},v}-\xi^{(n-1)}_{v_{-},v}\right)\right]^{\mathrm{T}}$
(5)

Note that neither vector (2) nor (5) are nodal variables, but carry
information across the edge $(v_{-},v)$. Figure 1 shows the relationship
between a continuous (in this case generation) function and it’s coefficients,
along with it’s derivative.

We construct the scenario tree by recusively applying the well known $k$-means
clustering algorithm to a bundle of net-load sample paths. We consider $H$
stages (hours), and have a $\mathbb{R}^{H\times(n+1)}$ dimensional sample
space of polynomial (spline) coefficients. The sample space is partitioned
into Voronoi cells that are constrained to satisfy the structure imposed by
filtration $\mathfrak{F}$ and the piece-wise continuity of polynomial
functions. Considering $L$ net-load sample paths, we approximate each with a
piece-wise continuous polynomial of degree $n$ using the Bernstein basis. The
vector of coefficients for each path and each interval are stored in the
tensor $\boldsymbol{\Xi}^{s}\in\mathbb{R}^{L\times H\times(n+1)}$ spanning the
$H$ stages. The total number of tree nodes is fixed ahead of time as
$|\mathcal{V}|=N+1$, and the number of nodes per stage determined by the
entries of the non decreasing $1\times H$ vector $\boldsymbol{c}$. The edges
${\cal E}$ are allocated optimally while quantizing $\boldsymbol{\Xi}^{s}$.
The structure of the optimum tree is stored in an ancestry (parent) vector
$\boldsymbol{a}$ such that node $v$ is a child of node $a^{1}_{v}$, and
grandchild of node $a^{2}_{v}$. The non-anticipativity property is enforced by
finding an appropriate mapping between each scenario
$\boldsymbol{\xi}^{s}_{\ell}\in\boldsymbol{\Xi}^{s}$ and the tree that ties
the past values to the future samples. To summarize, for the scenario tree
construction, we map our selection of sample discrete-time net-load
trajectories to their continuous time approximations, and then merge paths
until we get our desired tree structure. Figure 3 shows trees constructed for
our numerical comparison.

##  III The Continuous Time Multi-Stage Reserve and Unit Commitment Problem
(CT-MSRUC)

The major differences between a conventional UC formulation and our CT-MSRUC
are (a) we consider continuous load and generation trajectories, approximated
through Bernstein coefficients, and (b) the stochastic UC considers multiple
future scenarios. The stochastic nature of the formulation means only a
handful of decision variables are defined w.r.t. time (i.e. have subscript
$t$), instead, most are defined in _nodal_ form, and thus use subscript $v$ to
identify that they refer to tree node $v$ or edge $(v_{-},v)$.

###  III-A Considerations Regarding Optimizing in Continuous-Time

Ideally, in the continuous-time UC formulation, generation, ramping,
commitment, start-up and shut-down decision variables $\mathrm{x}^{g}(t)$,
$\dot{\mathrm{x}}^{g}(t)$, $\mathrm{y}^{g}(t)$, $\overline{\mathrm{s}}^{g}(t)$
and $\underline{\mathrm{s}}^{g}(t)$, for all generation units
$g\in\mathcal{G}$, may vary continually throughout time $t$. This provides the
system with ultimate flexibility to balance the load in an optimal way. Now
that the continuous load process and its derivative on the tree are
identified, we will extend the idea into the unit commitment formulation for
the generation trajectories. We introduce:

$\displaystyle{\bf x}^{g}_{v_{-},v}$
$\displaystyle=[x^{g(0)}_{v_{-},v},x^{g(1)}_{v_{-},v},x^{g(2)}_{v_{-},v}\ldots,x^{g(n-1)}_{v_{-},v},x^{g(n)}_{v_{-},v}]^{\mathrm{T}}$
(6) $\displaystyle\dot{{\bf x}}^{g}_{v_{-},v}$
$\displaystyle=[\dot{x}^{g(0)}_{v_{-},v},\dot{x}^{g(1)}_{v_{-},v},\dot{x}^{g(2)}_{v_{-},v}\ldots,\dot{x}^{g(n-1)}_{v_{-},v}]^{\mathrm{T}}$
(7)

where ${\bf x}^{g}_{v_{-},v}$ and $\dot{{\bf x}}^{g}_{v_{-},v}$ are the
continuous-time polynomial coefficients for generation and ramping
trajectories respectively. It is clear that the elements of vector $\dot{{\bf
x}}^{g}_{v_{-},v}$ can be expressed as a linear combination of elements of
generation, $\dot{\bf{x}}^{g}_{v_{-},v}=n{\bf M}{\bf{x}}^{g}_{v_{-},v}$, and
the corresponding continuous-time generation and ramp functions are:

$\displaystyle\mathrm{x}^{g}_{v_{-},v}(t)$ $\displaystyle={\bf
b}_{n}\bigg{(}\frac{t-t_{v_{-}}}{t_{v}-t_{v_{-}}}\bigg{)}{\bf
x}^{g}_{v_{-},v}$ (8) $\displaystyle\dot{\mathrm{x}}^{g}_{v_{-},v}(t)$
$\displaystyle=n{\bf M}{\bf
b}_{n-1}\bigg{(}\frac{t-t_{v_{-}}}{t_{v}-t_{v_{-}}}\bigg{)}\dot{\bf
x}^{g}_{v_{-},v}$ (9)

We assume that the commitment and therefore start-up and shut-down variables
are constant within each interval $t_{v_{-}}\leq t<t_{v}$ and equal to the
commitment, start-up and shut-down decisions at the _end_ of the interval. The
reason being that the control point at the end of the interval
$[t_{v_{-}},t_{v}]$, which we conventionally denote by $y^{g}_{v}$,
$\overline{s}^{g}_{v}$ and $\underline{s}^{g}_{v}$ respectively, carries all
the information of the edge $(v_{-},v)$.

Another useful property of Bernstein polynomials states that the entire
generation and ramping trajectories for an edge $(v_{-},v)$ are contained
within the convex hull of their control points ${\bf x}^{g}_{v_{-},v}$ and
$\dot{{\bf x}}^{g}_{v_{-},v}$ respectively. Thus, we can bound generation
dispatch and ramping for
$g\in\mathcal{G},v\in\mathcal{V}^{+}=\mathcal{V}\smallsetminus\\{0\\}$ (all
tree nodes excluding the root), and $t\in[t_{v_{-}},t_{v}]$:

$\displaystyle\min\\{{\bf{x}}^{g}_{v_{-},v}\\}$
$\displaystyle\leq{\mathrm{x}}^{g}_{v_{-},v}(t)\leq\max\\{{\bf{x}}^{g}_{v_{-},v}\\}$
(10) $\displaystyle\min\\{{\bf{\dot{x}}}^{g}_{v_{-},v}\\}$
$\displaystyle\leq\dot{\mathrm{x}}^{g}_{v_{-},v}(t)\leq\max\\{{\bf{\dot{x}}}^{g}_{v_{-},v}\\}$
(11)

###  III-B Design of the combined reserve and power market

Conventional unit commitment problems are solved several hours ahead of the
time horizon they apply to, and should provide generator owners with
production schedules that include expected dispatch, commitment as well as any
reserve power obligations. Solving such a problem while considering several
possible load outcomes (mapped to a tree) naturally yields several distinct
dispatch solutions that may or may not include variations in unit commitment,
depending on how the problem is constrained. The outcome is a set of
decisions, some with respect to time, but others with respect to the different
paths of the tree, which is hard to interpret as a market decision. This is
why we propose to merge the commitment of reserves together with the
traditional procurement of generation to meet the base-load, framing all of
these decisions as part of our stochastic optimization variables.

Beside the continuous generation trajectories, the following summarizes the
difference between the current operational norm and our formulation, as seen
from the perspective of a generator owner:
1) The continuous generation and commitment _schedule_ is the solution that
corresponds to the most likely path on the scenario tree. This is denoted by
the spline $x^{g,\text{schedule}}(t)$ and the vector
$y^{g,\text{schedule}}\in\\{0,1\\}^{H}$ for generation and commitment
respectively, where $H$ is the number of hours.
2) The variable $\overline{y}^{g}\in\\{0,1\\}^{H}$ indicates at which hour $h$
unit $g$ _may_ be committed. As such, $\overline{y}_{h}^{g}=1$ indicates that
for one or more of the paths of the scenario tree that cross hour $h$, the
unit is committed, and the generator owner needs to anticipate this. Note that
a commitment profile will always honor minimum on/off constraints.
3) The continuous reserve capacity of units is denoted by
$\underline{r}^{g}(t)$ and $\overline{r}^{g}(t)$, generator $g$ should
anticipate to be operating in the range
$[x^{g,\text{schedule}}(t)-\underline{r}^{g}(t),x^{g,\text{schedule}}(t)+\overline{r}^{g}(t)]$
at time $t$.
4) We split payments for commitment into (a) payment for the possibility of
being committed during a particular hour $\overline{Y}^{g}$ and (b) the real-
time payment for actually being committed $Y^{g}$.

Note that as we commit reserves it would be natural to also include N-1
contingencies, but for simplicity we omit this in the CT-MSRUC. The extension
is straightforward.

###  III-C Program Formulation

To reduce complexity, the formulation does not consider transmission
constraints, assuming a single-bus system. We assume the problem spans $H$
stages (hours), defining the set $\mathcal{H}\in\\{1,2,\ldots,H\\}$. We define
the set of all scenario tree nodes $\mathcal{V}$, the set of non-root nodes
$\mathcal{V}^{+}=\mathcal{V}\smallsetminus\\{0\\}$, the set of nodes at stage
$h$ as $\mathcal{Z}(h)$, and the preceding stage ${h-1=h_{-}}$. We assume the
root node to be the only member of $\mathcal{Z}(0)$. We use the subscript $h$
to refer to the time instant at the end of stage $h$, which pertains to
decisions made in the interval $t_{h_{-}}\leq t<t_{h}$. The edge between
consecutive nodes $v_{-}$ and $v$ spans the time $t_{v_{-}}\leq t<t_{v}$. As
in the conventional UC formulations, we will consider this interval equal to 1
hour. To make our notation more compact, we consider the case of cubic
Bernstein polynomials ($n=3$). The generalization to an arbitrary $n$ is
straightforward. The following sections write out the nodal unit commitment
formulation.

####  III-C1 Continuity Constraints

We enforce $C^{0}$ and $C^{1}$ continuity across the edge polynomials. For all
$g\in\mathcal{G}$:

$\displaystyle x_{v_{-},v}^{g(n)}=x_{v,v_{+}}^{g(0)}\quad\forall
v\in\mathcal{V}^{+},v_{+}\in\mathcal{C}(v),v\not\in\mathcal{Z}(H)$ (12)
$\displaystyle\dot{x}_{v_{-},v}^{g(n-1)}=\dot{x}_{v,v_{+}}^{g(0)}\quad\forall
v\in\mathcal{V}^{+},v_{+}\in\mathcal{C}(v),v\not\in\mathcal{Z}(H)$ (13)

####  III-C2 Reserve Capacity Constraints

We define the continuous reserve capacity in terms of edge variables as well
as the aforementioned hourly reserve capacity indicators. The edge-wise
continuous reserve coefficients $\hat{\mathrm{r}}_{v_{-},v}^{g(i)}$ and
$\check{\mathrm{r}}_{v_{-},v}^{g(i)}$ define the necessary reserve capacity
required of the units committed on that edge, while
$\overline{\mathrm{r}}_{h_{-},h}^{g(i)}$ and
$\underline{\mathrm{r}}_{h_{-},h}^{g(i)}$ define the continuous hourly reserve
requirements considering all edges $v\in\mathcal{Z}(h)$ spanning
$t_{v_{-}}\leq t<t_{v}$. For all $g\in\mathcal{G}$:

$\displaystyle{\bf\overline{r}}_{h_{-},h}^{g}\geq{\bf
x}_{v_{-},v}^{g}+{\bf{\hat{r}}}_{v_{-},v}^{g}-{\bf
x}^{g,\text{schedule}}_{h_{-},h}\quad\forall
h\in\mathcal{H},v\in\mathcal{Z}(h)$ (14)
$\displaystyle{\bf\underline{r}}_{h_{-},h}^{g}\geq{\bf
x}^{g,\text{schedule}}_{h_{-},h}-{\bf{x}}_{v_{-},v}^{g}+{\bf{\check{r}}}_{v_{-},v}^{g}\quad\forall
h\in\mathcal{H},v\in\mathcal{Z}(h)$ (15)

Here, we assume the continuous functions
$\underline{r}^{g}(t),\overline{r}^{g}(t)$ and $x^{g,\text{schedule}}(t)$ are
modelled using splines, in the same fashion as the edge-wise generation
$x^{g}(t)$, thus (note that $t_{h}-t_{h_{-}}=1$):

$\displaystyle\underline{r}^{g}(t)$ $\displaystyle={\bf
b}_{n}(t-t_{h_{-}}){\bf\underline{r}}_{h_{-},h}^{g}\quad t_{h_{-}}\leq t\leq
t_{h}$ (16) $\displaystyle\overline{r}^{g}(t)$ $\displaystyle={\bf
b}_{n}(t-t_{h_{-}}){\bf\overline{r}}_{h_{-},h}^{g}\quad t_{h_{-}}\leq t\leq
t_{h}$ (17) $\displaystyle x^{g,\text{schedule}}(t)$ $\displaystyle={\bf
b}_{n}(t-t_{h_{-}}){\bf x}_{h_{-},h}^{g,\text{schedule}}\quad t_{h_{-}}\leq
t\leq t_{h}$ (18)

The binary $y_{v}^{g}$ denotes whether a unit is committed along the preceding
edge $(v_{-},v)$ and it’s relationship with the hourly commitment indicator
$\overline{y}_{h}^{g}$ is:

$\displaystyle\overline{y}_{h}^{g}\geq y_{v}^{g}\quad\forall
h\in\mathcal{H},v\in\mathcal{Z}(h),g\in\mathcal{G}$ (19)

####  III-C3 Balance Constraints

We maintain a balance of generation and load by ensuring the coefficients are
equal:

$\displaystyle\textstyle\sum_{g\in\mathcal{G}}{\bf
x}^{g}_{v_{-},v}-{\boldsymbol{\xi}}_{v_{-},v}=0\quad\forall
v\in\mathcal{V}^{+}$ (20)

Furthermore, we denote the estimated root mean square error between
$\boldsymbol{\xi}_{v_{-},v}$ and the bundle of scenarios it approximates as
$\boldsymbol{\epsilon}_{v_{-},v}$. We can require reserves to cover a certain
range around the centroid, proportional to the mean square error:

$\displaystyle\\!\\!\\!\\!\textstyle\sum_{g\in\mathcal{G}}({\bf{x}}^{g}_{v_{-},v}+{\bf{\hat{r}}}_{v_{-},v}^{g})-{\boldsymbol{\xi}}_{v_{-},v}-\rho\,{\boldsymbol{\epsilon}}_{v_{-},v}$
$\displaystyle=0\quad\forall v\in\mathcal{V}^{+}$ (21)
$\displaystyle\\!\\!\\!\\!\textstyle\sum_{g\in\mathcal{G}}({\bf{x}}^{g}_{v_{-},v}-{\bf{\check{r}}}_{v_{-},v}^{g})-{\boldsymbol{\xi}}_{v_{-},v}+\rho\,{\boldsymbol{\epsilon}}_{v_{-},v}$
$\displaystyle=0\quad\forall v\in\mathcal{V}^{+}$ (22)

For an unbiased centroid $\boldsymbol{\xi}$, $\boldsymbol{\epsilon}$
corresponds to the standard deviation of the sample paths. Essentially $\rho$
allows us to put a margin equal to a multiple of the estimated conditional
standard deviation, which can be used as a bound for the probability of
deviating from the $\boldsymbol{\xi}$.

####  III-C4 Generation and Ramping Limits

Generation (MW)$0$$50$$100$$150$$200$$\overline{P}^{g}$$\underline{P}^{g}$Ramp
speed (MW/hr)Time
(hours)$0$$1$$2$$3$$4$$5$$-300$$-200$$-100$$0$$100$$200$$300$$\overline{R}^{g}$$-\overline{R}^{g}$$y_{1}^{g}=1$$\overline{s}_{1}^{g}=0$$\underline{s}_{1}^{g}=0$$y_{2}^{g}=1$$\overline{s}_{2}^{g}=0$$\underline{s}_{2}^{g}=0$$y_{3}^{g}=0$$\overline{s}_{3}^{g}=0$$\underline{s}_{3}^{g}=1$$y_{4}^{g}=0$$\overline{s}_{4}^{g}=0$$\underline{s}_{4}^{g}=0$$y_{5}^{g}=1$$\overline{s}_{5}^{g}=1$$\underline{s}_{5}^{g}=0$$x_{0,1}^{g(2)}$$r_{0,1}^{g(1)}$
Figure 1: The Bernstein coefficients of cubic spline interpolation for
generation and ramping (starred), along with the commitment status and turn
on/off indicators. Note how a turn-on or turn-off indicator relaxes the
ramping constraint (dotted line) for the center ramping coefficient of the
preceding hour.

First off, we define the relationship between generation and ramping:

$\displaystyle\dot{{\bf{x}}}^{g}_{v_{-},v}=n{\bf{M}}{\bf{x}}^{g}_{v_{-},v}\quad\forall
v\in\mathcal{V}^{+},g\in\mathcal{G}$ (23)

Now, for a unit that is offline for a certain edge (hour) $v_{-}$ and online
for the consecutive edge $v$, we can not simply force all the coefficients
${\bf{x}}_{v_{-},v}$ to be zero and ${\bf{x}}_{v,v_{+}}$ to be non-zero, as
this violates the continuity constraints. Instead, we allow units to turn-on
and turn-off as late as possible during the preceding hour. This translates
into the following constraints for all units $g\in\mathcal{G}$, assuming
$n=3$:

$\displaystyle x_{v_{-},v}^{g(i)}+\hat{r}_{v_{-},v}^{g(i)}$
$\displaystyle\leq\overline{P}^{g}y_{v}^{g}\quad\forall
v\in\mathcal{V}^{+},i\in\\{0,1\\}$ (24) $\displaystyle
x_{v_{-},v}^{g(i)}-\check{r}_{v_{-},v}^{g(i)}$
$\displaystyle\leq\underline{P}^{g}y_{v}^{g}\quad\forall
v\in\mathcal{V}^{+},i\in\\{0,1\\}$ (25) $\displaystyle
x_{a^{2}_{v},a^{1}_{v}}^{g(i)}+\hat{r}_{a^{2}_{v},a^{1}_{v}}^{g(i)}$
$\displaystyle\leq\overline{P}^{g}y_{v}^{g}\quad\forall
v\in\mathcal{V},v\not\in\mathcal{Z}(1),i\in\\{2,3\\}$ (26) $\displaystyle
x_{a^{2}_{v},a^{1}_{v}}^{g(i)}-\check{r}_{a^{2}_{v},a^{1}_{v}}^{g(i)}$
$\displaystyle\leq\underline{P}^{g}y_{v}^{g}\quad\forall
v\in\mathcal{V},v\not\in\mathcal{Z}(1),i\in\\{2,3\\}$ (27) $\displaystyle
x_{v_{-},v}^{g(i)}+\hat{r}_{v_{-},v}^{g(i)}$
$\displaystyle\leq\overline{P}^{g}y_{v}^{g}\quad\forall
v\in\mathcal{Z}(H),i\in\\{2,3\\}$ (28) $\displaystyle
x_{v_{-},v}^{g(i)}-\check{r}_{v_{-},v}^{g(i)}$
$\displaystyle\leq\underline{P}^{g}y_{v}^{g}\quad\forall
v\in\mathcal{Z}(H),i\in\\{2,3\\}$ (29)

recalling the ancestry vector $\boldsymbol{a}$ introduced in Section II and
the convention $a^{0}_{v}:=v$. Generalizing this to any $n$ involves replacing
the sets $\\{0,1\\}$ with $\\{0,\ldots,n-C\\}$ and $\\{2,3\\}$ with
$\\{n-C+1,\ldots,n\\}$ where $C$ is the level of continuity ($0$ for $C^{0}$
and so on). The ramping constraints become inter-hourly, but we must ensure
that they are relaxed for on-off transitions. For all $g\in\mathcal{G}$:

$\displaystyle-\overline{R}^{g}$
$\displaystyle\leq\dot{x}_{v_{-},v}^{g(i)}\leq\overline{R}^{g}\quad\forall
v\in\mathcal{V}^{+},i\in\\{0,2\\}$ (30)
$\displaystyle-\overline{R}^{g}\\!\\!-\\!n\overline{P}^{g}\underline{s}_{v}^{g}$
$\displaystyle\leq\dot{x}_{a^{2}_{v},a^{1}_{v}}^{g(1)}\leq\overline{R}^{g}\\!\\!+\\!n\overline{P}^{g}\overline{s}_{v}^{g}\quad\forall
h\in\\{2,\ldots,H\\},v\in\mathcal{Z}(h)$ (31) $\displaystyle\overline{R}^{g}$
$\displaystyle\leq\dot{x}_{v_{-},v}^{g(1)}\leq\overline{R}^{g}\quad\forall
v\in\mathcal{Z}(H)$ (32)

Again here, to generalize this to any $n$ involves applying (31) and (32) to
coefficient $i=n-C$, and leaving the rest
($\\{0,\ldots,n-1\\}\smallsetminus\\{n-C\\}$) to (30). Figure 1 attempts to
visually explain constraints (24)-(32).

####  III-C5 Minimum On and Off Constrains

We formulate conventional minimum on-off constraints for all
$g\in\mathcal{G}$, where $O^{g}_{n}$ and $O^{g}_{f}$ denote the minimum on and
off time respectively:

$\displaystyle\overline{s}_{v}^{g}-\underline{s}_{v}^{g}$
$\displaystyle=y_{v}^{g}-y_{v_{-}}^{g}\quad\forall v\in\mathcal{V}^{+}$ (33)
$\displaystyle y_{v}^{g}$
$\displaystyle\geq\textstyle\sum_{u=0}^{\min(O_{n}^{g},h)-1}\overline{s}^{g}_{a^{u}_{v}}\quad\forall
h\in\mathcal{H},v\in\mathcal{Z}(h)$ (34) $\displaystyle 1-y_{v}^{g}$
$\displaystyle\geq\textstyle\sum_{u=0}^{\min(O_{f}^{g},h)-1}\underline{s}^{g}_{a^{u}_{v}}\quad\forall
h\in\mathcal{H},v\in\mathcal{Z}(h)$ (35)

####  III-C6 Objective Function

Integrating Bernstein polynomials is straight-forward (again assuming
$t_{v}-t_{v_{-}}=1$):

$\int_{t_{v_{-}}}^{t_{v}}C^{g}{\mathrm{x}}(t)\,dt=C^{g}\int_{t_{v_{-}}}^{t_{v}}{\bf{b}}_{n}(t-t_{v_{-}}){\bf{x}}_{v_{-},v}^{g}=\frac{C^{g}}{n+1}\boldsymbol{1}\cdot{\bf{x}}_{v_{-},v}^{g}$

In this case, the objective of the stochastic unit commitment is minimizing
the expected cost of generation and the reserve capacity payment. The
objective becomes:

$\begin{split}\min\sum_{h\in\mathcal{H}}\sum_{g\in\mathcal{G}}\left(\frac{1}{n+1}\sum_{i=0}^{n}\left[\overline{R}^{g}\overline{r}_{h_{-},h}^{g(i)}+\underline{R}^{g}\underline{r}_{h_{-},h}^{g(i)}\right]+\overline{Y}^{g}\overline{y}_{h}^{g}\right)\\\
+\sum_{v\in\mathcal{V}^{+}}\pi_{v}\sum_{g\in\mathcal{G}}\left(Y^{g}y_{v}^{g}+\overline{S}\overline{s}^{g}_{v}+\underline{S}\underline{s}^{g}_{v}+\frac{X^{g}}{n+1}\sum_{i=0}^{n}x_{v_{-},v}^{g(i)}\right)\end{split}$
(36)

where $\pi_{v}$ is the probability of node $v$, $\overline{R}$ and
$\underline{R}$ are up- and down-spinning reserve costs, $\overline{S}$ and
$\underline{S}$ are startup and shutdown costs, $X$ is the linear generation
cost term, $Y$ is the constant (commitment) cost term and $\overline{Y}$ is
the payment for possibly being committed for a particular hour, a reserve
payment of sorts. The time-based segment of (36) denotes up-front payments for
reserves and anticipated operation range, while the nodal based segment is the
expected real-time cost. Without information about real-time prices, we simply
substitute real-time prices with the day-ahead bids, as the best indicator of
what will happen in real-time. However, we wish to remark that having a higher
real-time cost will only make the numerical comparison more favorable to our
framework.

##  IV Numerical Results

###  IV-A Load Scenario Tree

$1700$$1800$$1900$$2000$18:0019:0020:0021:00Power (MW)Time of day (hh:mm)
Figure 2: The spline approximation of three training scenarios for two hours.
The sample data is denoted by diamonds, the zero order piecewise constant
spline is dotted while the $C^{1}$ continuous cubic spline is solid. We see
that the cubic spline is already an extremely good approximation of the sample
path.

To illustrate the proposed framework we use $L=1122$ sample days of _scaled_
and _aggregated_ CAISO net-load measurements spanning the years 2013-2017,
with a 5-minutes resolution[10]. We find the spline approximations of the load
scenarios, where each polynomial spans a single hour of load. We make two
approximations: (a) cubic spline with $C^{1}$ continuity forced across the
nodes and for comparison (b) zero order (discontinuous piece wise constant)
spline. The zero order polynomial reflects standard unit commitment practice,
where the load is averaged over each hour. To solve for the spline
coefficients we use a standard least-squares approach. As an example, Figure 2
captures these approximations for three sample trajectories. We split the
input scenarios (days) into two sets, 70% we use for the tree construction and
the remaining 30% we use as sample days for verifying the solution. The tree
is constructed for a “generic” load day using a simple recursive
$\operatorname{\mathnormal{k}-means}$ algorithm which solves in a matter of
seconds, but could equally be generated from a more targeted dataset, such as
a set of load predictions based on weather forecasts and historical load data.
For both zero and third order spline approximations the training scenarios are
reduced down to the trees shown in Figure 3. The figure shows the root
(arbitrarily chosen to be at hour $-1$), the knots of the polynomials (dotted
vertical lines), the allowed branching points (solid vertical lines) and the
number of nodes per stage indicated by the $c_{h}$ value on top. The different
colored paths show the reduced sample path centroids, with the root mean
square error indicated by the surrounding shaded region.

(a) The zero-order spline scenario tree.

(b) The cubic spline scenario tree.

Figure 3: The zero and third order spline approximation scenario trees (solid
colored lines) along with the root mean square error of the scenario bundle
yielding the corresponding centroid (shaded region). The $c_{h}$ indicator on
top shows how many nodes are defined for each stage.

###  IV-B Solving the Stochastic Unit Commitment

For generation data we use a single area of the IEEE RTS-96 test case [11],
but consider it to be a single-bus system. As mentioned before, the CAISO load
data is scaled by a factor of $16$ in order to fit into the range of
generation available in the RTS system. The cost terms used in the objective
function (36) are summarized in Table I for the different RTS case generation
types. These terms are based on the original case data, however, a piece wise
linear objective is assumed and the quadratic component of the generator cost
curves are ignored. Fixed capacity costs ($\overline{R}$, $\underline{R}$,
$\overline{Y}$) are proportional to the linear and constant cost terms.

Table I: Cost terms of the different RTS case generators. Unit Type | $\overline{S}$ | $\underline{S}$ | $Y$ | $X$ | $\overline{R}$ | $\underline{R}$ | $\overline{Y}$
---|---|---|---|---|---|---|---
U12 | 1500 | 0 | 86.4 | 56.6 | 28.3 | 28.3 | 21.6
U20 | 1500 | 0 | 400.7 | 130.0 | 65.0 | 65.0 | 100.2
U50 | 1500 | 0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0
U76 | 1500 | 0 | 212.3 | 16.1 | 8.0 | 8.0 | 53.1
U100 | 1500 | 0 | 781.5 | 43.7 | 21.8 | 21.8 | 195.4
U155 | 1500 | 0 | 382.2 | 12.4 | 6.2 | 6.2 | 95.6
U197 | 1500 | 0 | 832.8 | 48.6 | 24.3 | 24.3 | 208.2
U350 | 1500 | 0 | 665.1 | 11.8 | 5.9 | 5.9 | 166.3
U400 | 1500 | 0 | 395.4 | 4.4 | 2.2 | 2.2 | 98.8

Using a net-load tree as the input, we solve the unit commitment problem as
described by (12)-(36), for both the cubic spline approximation (CT-MSRUC) and
the discrete time formulation (DT-MSRUC). The latter is simply a special case
of the described formulation, where with respect to the constraints we can
imagine the load and generation to be first-order $C^{0}$ continuous, thus
enforcing the ramping constraints and having the same number of free variables
as a conventional zero order formulation. Modeling was done using Python [12,
13, 14, 15], and for solving we used Gurobi 7.5 [16] with a MIP gap parameter
of $5\%$. The solution time of the DT-MSRUC problem was 40 minutes while for
the CT-MSRUC it solved in 67 hours. The increase in computational time is
expected because the number of continuous variables is roughly quadrupled,
while the number of integer variables remains the same. Reducing the solver
time drastically through decomposition and parameter tuning should be
relatively straightforward, but was not the focus of this work. As an example,
using a $\rho=2.5$ solved in a couple of hours.

(a) DT-MSRUC

(b) CT-MSRUC

Figure 4: Unit commitment solutions for the DT-MSRUC and CT-MSRUC
formulations. Purple (dark) region denotes units that are never committed,
green (gray) denotes units that are committed for one or more paths of the
solution, while yellow (light) denotes units committed for the scheduled
solution.

We solved the two unit commitment problems CT-MSRUC and DT-MSRUC with the
reserve parameter $\rho=3$ and applied the solutions to the set of test
scenarios, considering the real-time cost and potential infeasibility issues.
Figure 4 visualizes the unit commitment of both solutions. The purple shows
uncommitted units for any path during a particular hour and yellow shows
committed units for the _scheduled_ (most likely) path, while green shows what
is not committed for the most likely load, but is committed for one or more of
the alternative load trajectories. We observe that the commitment is fairly
similar between the two cases, the hydro (U50) and nuclear (U400) units are
committed for most paths and a handful of oil and coal plants (U12, U76, U100,
U197) are brought online for the afternoon hours of several paths. A notable
difference lies in the usage of the large coal plants (U155, U350) which are
scheduled differently for the two solutions.

Figure 5 compares several aspects of the two problems. First, the expected
energy, commitment and reserve costs as formulated in (36) are very similar,
with the continuous solution giving a slightly higher expected cost. However,
when we run through the set of test scenarios the opposite seems true, the
mean cost of commitment and energy turns out to be lower (_mean testing_ bar)
and when we add the reserve cost (_total testing_) we see that the CT-MSRUC
solution is actually cheaper, as it offers a more efficient way to service the
realistic test trajectory which is continuous in nature. Moreover, 44% of the
testing load trajectories are simply not feasible for the DT-MSRUC solution,
given the commitment and reserve capacity of scheduled units, compared to a
level of 14% infeasibility for the CT-MSRUC solution.

In real-time, the system operator must maintain generation–consumption
balance. As the load prediction is never perfect, the exact dispatch solution
of the UC is not followed, but the commitment and reserve capacity ”schedule”
should be known ahead of real-time. For the case of stochastic UC, however,
each path of the scenario tree will contain a different commitment schedule.
Observing the load in real-time, we follow the tree to obtain the specific
commitment profile, essentially throwing out any diverging paths until we end
up with what looks like a solution to a deterministic UC. Figure 6 shows a
given test load trajectory (solid green) and the nearest load path on the
scenario tree (dotted blue), compared to the most likely path on the tree
(dashed yellow). While the load will not exactly match the branch of the tree
we follow, we do know the commitment, scheduled dispatch and reserve capacity,
and can transform these into the range of load that our solution can service
(shaded gray region of figure). Here we can observe the CT-MSRUC solution
having two major advantages over the DT-MSRUC solution, leading to higher
likelihood of the solution being feasible for a real load trajectory. First,
we do a much better job of incorporating realistic ramps, so we can better
anticipate load trajectories that change abruptly and have the right fleet of
units available, and secondly, by not modeling in terms of averaged-load of an
hour, we are less likely to ”cut corners” where the test load trajectory falls
outside the feasible region for a brief time. This can be seen in Figure 6
where the discrete-time solution can not deliver desired load during the first
minutes of hour 22.

Figure 5: The various cost components of the solution and the test scenarios.

(a) Applying the DT-MSRUC solution.

(b) Applying the CT-MSRUC solution.

Figure 6: Generator schedule (dashed orange) along with the dispatch (dotted
blue), commitment and reserve solution given along the path of the load
scenario tree that closest resembles a given test load scenario (solid green).
The reserve capacity of committed units is denoted by the shaded gray region.

##  V Conclusions & Future Work

This paper describes a stochastic continuous-time unit commitment formulation,
based on Bernstein polynomials. We briefly discuss how a set of input
scenarios are reduced into a scenario tree, and formulate a multi-stage
stochastic unit commitment, where commitment of units varies between the
different paths of the tree. We show that our formulation outperforms
conventional stochastic discrete-time unit commitment, as we better account
for necessary real-time ramps and we are less likely to underestimate reserves
needed for the hours where load is sharply increasing or decreasing.

Future work will go more in depth into the tree construction, focus on
extending this to multi-bus systems, incorporating the spatial correlation of
load into the scenario tree, as well as looking at the trade-off between tree
sizes, branching structure and higher order polynomials. Further, a natural
progression of this work is incorporating decomposition methods into the
solution algorithm to compensate for the unavoidable increase in complexity.

## References

  * [1] R. Wiebking, “Stochastische modelle zur optimalen lastverteilung in einem kraftwerksverbund,” _Mathematical Methods of Operations Research_ , vol. 21, no. 6, pp. B197–B217, 1977.
  * [2] S. Takriti, J. R. Birge, and E. Long, “A stochastic model for unit commitment problem,” _IEEE Transactions on Power Systems_ , vol. 11, no. 3, pp. 1497–1508, 1996.
  * [3] P. Carpentier, G. Cohen, J.-C. Culioli, and A. Renaud, “Stochastic optimization of unit commitment: a new decomposition framework,” _IEEE Transactions on Power System_ , vol. 11, pp. 1067–1073, 1996.
  * [4] M. Nowak and W. Römisch, “Stochastic lagrangian relaxation applied to power scheduling in a hydro-thermal system under uncertainty,” _Annals of Operations Research_ , vol. 100, no. 251, 2000.
  * [5] T. Shiina and J. R. Birge, “Stochastic unit commitment problem,” _International Transaction on Operations Research_ , vol. 11, pp. 19–32, 2004\.
  * [6] A. Papavasiliou, S. Oren, and R. O’Neill, “Reserve requirements for wind power integration:a scenario-based stochastic programming framework,” _IEEE Transactions on Power Systems_ , vol. 26, no. 4, pp. 2197–2206, 2011.
  * [7] F. Bouffard and F. Galiana, “Stochastic security for operations planning with significant wind power generation,” _IEEE Transaction on Power System_ , vol. 23, no. 2, pp. 306–316, 2008.
  * [8] M. Parvania and A. Scaglione, “Unit commitment with continuous-time generation and ramping trajectory models,” _IEEE Transactions on Power Systems_ , vol. 31, no. 4, pp. 3169–3178, 2016.
  * [9] G. C. Pflug and A. Pichler, “Dynamic generation of scenario trees,” _Computational Optimization and Applications_ , vol. 62, no. 3, pp. 641–668, 2015.
  * [10] C. ISO, “California iso open access same-time information system (oasis),” 2017\. [Online]. Available: http://oasis.caiso.com
  * [11] C. Grigg, P. Wong, P. Albrecht, R. Allan, M. Bhavaraju, R. Billinton, Q. Chen, C. Fong, S. Haddad, S. Kuruganty _et al._ , “The ieee reliability test system-1996. a report prepared by the reliability test system task force of the application of probability methods subcommittee,” _IEEE Transactions on power systems_ , vol. 14, no. 3, pp. 1010–1020, 1999.
  * [12] S. Van Der Walt, S. C. Colbert, and G. Varoquaux, “The numpy array: a structure for efficient numerical computation,” _Computing in Science & Engineering_, vol. 13, no. 2, pp. 22–30, 2011.
  * [13] E. Jones, T. Oliphant, and P. Peterson, “$\\{$SciPy$\\}$: Open source scientific tools for $\\{$Python$\\}$,” 2014.
  * [14] J. D. Hunter, “Matplotlib: A 2d graphics environment,” _Computing in science and engineering_ , vol. 9, no. 3, pp. 90–95, 2007.
  * [15] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay, “Scikit-learn: Machine learning in Python,” _Journal of Machine Learning Research_ , vol. 12, pp. 2825–2830, 2011.
  * [16] I. Gurobi Optimization, “Gurobi optimizer reference manual,” 2017. [Online]. Available: http://www.gurobi.com


BaseMiniExclam[7]#1

#4 $\displaystyle\bodyobj{#2}{#3}{#6}{#5}\BODY$ (1a)

# Linear Programming Bounds for
Almost-Balanced Binary Codes

Venkatesan Guruswami and Andrii Riazanov  Carnegie Mellon University
Computer Science Department
Pittsburgh, PA 15213
Email: {venkatg<EMAIL_ADDRESS>
###### Abstract

We revisit the linear programming bounds for the size vs. distance trade-off
for binary codes, focusing on the bounds for the almost-balanced case, when
all pairwise distances are between $d$ and $n-d$, where $d$ is the code
distance and $n$ is the block length. We give an _optimal_ solution to
Delsarte’s LP for the almost-balanced case with large distance
$d\geq(n-\sqrt{n})/2+1$, which shows that the optimal value of the LP
coincides with the Grey-Rankin bound for self-complementary codes.

We also show that a limitation of the asymptotic LP bound shown by
Samorodnitsky, namely that it is at least the average of the first MRRW upper
bound and Gilbert-Varshamov bound, continues to hold for the almost-balanced
case.

##  I Introduction

This paper concerns the standard size-versus-distance tradeoff for binary
error-correcting codes. A binary code $C$ with block length $n$ and distance
$d$ is a subset of the $n$-dimensional Hamming cube $\\{0,1\\}^{n}$, such that
the Hamming distance between any two elements from $C$ is at least $d$. In
this paper such codes will be referred to as min-distance codes, and we denote
$A(n,d)$ to be the maximal size of such a code. Understanding the behavior of
$A(n,d)$ is one of the most fundamental challenges in combinatorial coding
theory. For the case when $d$ is linear in $n$, one is interested in
asymptotic rates of the code rather than their sizes:

$\displaystyle d,$ (2)

where $\delta\approx d/n\leq 1/2$ is the relative distance of the code.

The best known asymptotic bounds as well as finite-length bounds on $A(n,d)$
for many parameters are given by Delsarte’s _linear programming approach_ [1].
Specifically, McEliece, Rodemich, Rumsey, and Welch [2] obtain the best known
upper bounds on $R(\delta)$ for the entire region $\delta\in(0,1/2)$ by
finding a good feasible solution to Delsarte’s LP. There is, however, a
sizeable gap between the MRRW upper bound on $R(\delta)$ and the Gilbert-
Varshamov lower bound $R(\delta)\geq 1-h(\delta)$ where $h(\cdot)$ is the
binary entropy function, and shrinking this gap has remained a major challenge
for almost 45 years.

In this paper, we are interested in the performance of the linear programming
approach for bounding the size of _almost-balanced_ codes, which obey the
stronger condition that the Hamming distance between any two distinct
codewords lies between $d$ and $(n-d)$. Almost-balanced codes when
$d=(1-\epsilon)n/2$ are closely related to $\epsilon$-biased spaces that are
of fundamental interest in pseudorandomness and derandomization, starting with
the seminal work of Naor and Naor [3]. The recent breakthrough explicit
construction of high distance codes approaching the Gilbert-Varshamov bound by
Ta-Shma crucially proceeds by constructing almost-balanced codes [4].

Almost-balanced codes are closely related to _self-complementary_ codes. The
code $C^{\prime}$ is self-complementary if $u\in C^{\prime}$ implies that the
complementary vector $\overline{u}$ is also in $C^{\prime}$. The connection
between almost-balanced and self-complementary codes is two-way:

  * •

any almost-balanced code $C$ with distance $d$ (i.e. all pairwise distances
lie in the interval $[d,n-d]$) corresponds to a self-complementary code of
size $2|C|$. Indeed, the code $C^{\prime}=\\{u\,|\,u\text{ or
}\overline{u}\text{ is in }C\\}$ is self-complementary and has size $2|C|$.

  * •

any self-complementary code $C^{\prime}$ with distance $d$ corresponds to a
family of almost-balanced codes of size $|C^{\prime}|/2$ and distance $d$.
Such codes can be constructed by taking one codeword from every complementary
pair $u,\overline{u}\in C^{\prime}$.

Thus to translate the bounds for almost-balanced codes to self-complementary
codes, one just needs to multiply the size-versus-distance trade-off (for the
same distance $d$) by $2$.

Let us denote by $B(n,d)$ the maximal size of an almost-balanced code with
block length $n$ and distance $d$. Once again, when $d$ is linear in $n$, we
can study the asymptotic rate of the codes rather than their sizes:

$\displaystyle R^{\text{bal}}(\delta)$
$\displaystyle=\limsup\limits_{n\to\infty}\frac{\log B(n,\lfloor\delta
n\rfloor)}{n}\ .$ (3)

The driving force behind this work is to study upper bounds obtained for
almost-balanced codes via the linear programming method, and compare them to
the ones for the min-distance codes. Our results consists of two parts:

  1. 1.

For almost-balanced codes with large distances ${d\geq\frac{n-\sqrt{n}}{2}+1}$
we find an _optimal solution_ to Delsarte’s linear program. This solution
gives an upper bound on $B(n,d)$ equivalent to the Gray-Rankin bound for self-
complementary codes [5, 6, 7, 8]. Since the solution we obtain is optimal,
this proves a matching _lower bound_ on Delsarte’s LP problem, showing that LP
approach cannot prove better bounds for the almost-balanced case. While this
was already implied for a very limited range of pairs $(n,d)$ for which Grey-
Rankin bound is achievable ([9]), our bound works for arbitrary even $n,d$
satisfying ${d\geq\frac{n-\sqrt{n}}{2}+1}$.

  2. 2.

We show (Theorem 2) that the asymptotic LP bound for almost-balanced codes is
at least the average of the first MRRW bound and Gilbert-Varshamov lower
bound. This is the analog of Samorodnitsky’s result [10] for min-distance
codes, and indicates that the (direct) linear programming approach cannot
attain the best known lower bound on $R^{\text{bal}}(\delta)$. The proof is a
simple adaptation of Samorodnitsky’s argument.

##  II Krawtchouk polynomials and Delsarte’s LP

For fixed $n$, the Krawtchouk polynomials are defined as

$K_{s}(x)=\sum_{j=0}^{s}(-1)^{j}\binom{x}{j}\binom{n-x}{s-j},\qquad
s=0,1,\dots,n.$ (4)

$K_{s}(x)$ is a degree-$s$ polynomial, and $K_{0}$, $K_{1}$, $\dots$, $K_{n}$
form an orthogonal family with respect to a measure
$\mu_{i}=\binom{n}{i}/2^{n}$:

$\sum_{j=0}^{n}\frac{\binom{n}{j}}{2^{n}}K_{s}(j)K_{t}(j)=\delta_{st}\binom{n}{s}.$
(5)

Notice that $K_{0}(x)\equiv 1$ and $K_{i}(0)=\binom{n}{i}$. Below are some of
the properties of Krawtchouk polynomials we will need.

###### Fact 1

$K_{s}(n/2-x)$ is an even function for even $s$, and is an odd function for
odd $s$.

###### Fact 2

$K_{n-s}(i)=K_{s}(i)$ for even $i$, and ${K_{n-s}(i)=-K_{s}(i)}$ for odd $i$.

###### Fact 3

Reciprocity property: $\binom{n}{s}K_{i}(s)=\binom{n}{i}K_{s}(i)$.

The Delsarte’s upper bound on $A(n,d)$ can be formulated as the following
Linear Programming optimization problem: a_i ≥0∑_k=0^n a_k ∑_k=0^n a_k
K_s(k)≥0, s = 0, 1, …, n a_0 = 1 a_k = 01 ≤k ≤(d-1) Denote
$\boldsymbol{\beta}(x)=\sum\limits_{k=0}^{n}\beta_{k}K_{k}(x)$. Then the the
dual of the above LP can be written as β_i ≥0β(0) β(u)≤0, u = d, d+1, …, n β_0
= 1 Denoting by $A_{LP}(n,d)$ the optimum values of the above two problems,
the LP bound claims that $A(n,d)\leq A_{LP}(n,d)$.

The optimal values of the LPs and asymptotic upper bound
$R_{LP}(\delta)=\limsup\limits_{n\to\infty}\frac{1}{n}\log
A_{LP}(n,\lfloor\delta n\rfloor)$ are not known yet, however upper bounds on
$A_{LP}(n,d)$ (and $R_{LP}(\delta)$) can be obtained by finding feasible
solutions to the dual LP (3). In particular, the (first) MRRW upper bound
$R_{MRRW}(\delta)=H\left(\frac{1}{2}-\sqrt{\delta(1-\delta)}\right)$ is
obtained in such a way in [2] (in this paper we do not address the second MRRW
bound).

To formulate the Delsarte’s LP for the almost-balanced case, we point out that
the variables $a_{i}$, $i=0,1,\dots,n$ in the primal form of LP (3) correspond
to the distance distribution of the code $C$, i.e.
$a_{i}=\frac{1}{|C|}\left\lvert\\{(x,y)\in
C\,|\,\Delta(x,y)=i\\}\right\lvert$. The first set of constraints in the LP
(3) then correspond to Delsarte-MacWilliam’s inequalities for the transformed
(dual) distance distribution, while the constraints $a_{k}=0$ for $1\leq k<d$
denote that there are no codewords at distance below $d$.

For an almost-balanced code $C$ there is no two codewords in $C$ at distance
below $d$ or above $(n-d)$, i.e. $a_{k}=0$ for $k<d$ or $k>(n-d)$. Therefore,
the Delsarte’s linear program for almost-balanced binary codes with distance
$d$ is a_i ≥0∑_k=0^n a_k ∑_k=0^n a_k K_s(k)≥0, s = 0, 1, …, n a_0 = 1 a_k = 01
≤k ≤d-1 and k > (n-d). The dual to the above is the following linear program:
β_i ≥0β(0) β(u)≤0, u = d, d+1, …, n - d β_0 = 1 (the difference from the LP
(3) is that $\boldsymbol{\beta}(u)$ for $u>(n-d)$ are no longer required to be
non-positive.)

We denote by $B_{LP}(n,d)$ the optimal value of the above pair of LPs, and
thus $B(n,d)\leq B_{LP}(n,d)$. In analogy with $R_{LP}$, define by
$R^{\text{bal}}_{LP}(\delta)=\limsup\limits_{n\to\infty}\frac{1}{n}\log
B_{LP}(n,\lfloor\delta n\rfloor)$ the LP upper bound on the maximal asymptotic
rate for almost-balanced codes of relative distance $\delta$:
$R^{\text{bal}}(\delta)\leq R^{\text{bal}}_{LP}(\delta)$.

##  III LP bound for almost-balanced codes with large distances

In this section we find the exact optimal solutions to Delsarte’s linear
programs (3) and (3) for almost-balanced codes with large distance.
Specifically, we prove

###### Theorem 1

Let $n,d$ be even such that ${d\geq\frac{n-\sqrt{n}}{2}+1}$. Then the optimal
value of the linear programs 3 and 3 is
$B_{LP}(n,d)=\dfrac{4d(n-d)}{n-(n-2d)^{2}}$.

The upper bound on the size of almost-balanced linear codes $B(n,d)\leq
B_{LP}(n,d)$ for this regime is equivalent to the Grey-Rankin bound [5, 6, 7]
on the size of self-complementary codes. This bound was also obtained by
Delsarte by providing a _feasible_ solution to 3, see, for instance, Problem
(18) in Chapter 17 (p. 544) in [8], or [11] for a more general case. What we
show below is that this solution is in fact  _optimal_ for 3, and therefore
for 3 via duality.

This implies that the bound in Theorem 1 is the best upper bound on $B(n,d)$
one can obtain using (a straightforward application of) Delsarte’s LP, so the
Grey-Rankin bound cannot be improved using this approach. It was already shown
in the literature that the Grey-Rankin bound can sometimes be achieved [9],
however the range of pairs $(n,d)$ for which this happens is limited, and is
related to existence of designs with certain parameters. Theorem 1, on the
other hand, shows that the LP approach cannot improve the Grey-Rankin bound
for any even $n,d$. Besides that, we think that finding optimal (instead of
feasible) solutions to Delsarte’s LP is of independent interest, even for this
narrow regime and the special case of almost-balanced codes.

To prove Theorem 1, we start with the following claim

###### Lemma 1

Let $\theta=(\theta_{0},\theta_{1},\dots,\theta_{n})\in\mathbb{R}^{n+1}$ be a
feasible solution to the LP (3). Define $\overline{\theta}\in\mathbb{R}^{n+1}$
which coincides with $\theta$ on all even coordinates, and has $0$ elsewhere:
$\overline{\theta}=(\theta_{0},0,\theta_{2},0,\theta_{4},\dots)$. Then
$\overline{\theta}$ is a feasible solution to the LP (3) for almost-balanced
codes.

###### Proof:

Denote $\boldsymbol{\theta}(x)=\sum\limits_{k=0}^{n}\theta_{k}K_{k}(x)$, then
$\boldsymbol{\theta}(u)\leq 0$ for any integer $u$ s.t. $d\leq u\leq(n-d)$.
Consider then

$\overline{\boldsymbol{\theta}}(x)=\dfrac{\boldsymbol{\theta}(x)+\boldsymbol{\theta}(n-x)}{2}=\sum\limits_{k=0}^{n}\theta_{k}\dfrac{K_{k}(x)+K_{k}(n-x)}{2}.$

By Fact 1, $K_{k}(n/2-y)$ is an even polynomial for even $k$, and so
$\frac{K_{k}(x)+K_{k}(n-x)}{2}=K_{k}(x)$. On the other hand,
$\frac{K_{k}(x)+K_{k}(n-x)}{2}=0$ for odd $k$, and therefore
${\overline{\boldsymbol{\theta}}(x)=\sum\limits_{k=0}^{n}\overline{\theta}_{k}K_{k}(x)}$
by our definition of $\overline{\theta}$.

For any integer $d\leq u\leq n-d$, we also have $d\leq(n-u)\leq(n-d)$. It
means that $\boldsymbol{\theta}(u)\leq 0$ and $\boldsymbol{\theta}(n-u)\leq
0$, so $\overline{\boldsymbol{\theta}}(u)\leq 0$ for any such $u$. Therefore
$\overline{\theta}$ if feasible for (3). ∎

###### Corollary 1

Optimal solution to (3) has $\beta_{i}=0$ for odd $i$.

###### Proof:

$\theta_{i}\geq 0$ and $K_{k}(0)=\binom{n}{k}\geq 0$ in Lemma 1, thus
$\overline{\boldsymbol{\theta}}(0)=\sum\limits_{k=0}^{n}\overline{\theta}_{k}K_{k}(0)\leq\sum\limits_{k=0}^{n}\theta_{k}K_{k}(0)=\boldsymbol{\theta}(0)$.
So nullifying all the odd-indexed coordinates doesn’t increase the objective
value in 3. ∎

Now, for every even distance $d$ such that $\frac{n-\sqrt{n}}{2}+1\leq d\leq
n/2$ we find an optimal solution to (3) and its dual (3). The condition on $d$
comes from a restriction $|K_{3}(d)|\leq|K_{2}(d)|(n-2d)/3$, which we use in
the proof and prove in Lemma 4.

Observe that $K_{2}(x)=2x^{2}-2nx+\binom{n}{2}$ has roots
$\frac{n\pm\sqrt{n}}{2}$, and so $K_{2}(d)<0$ for distances $d$ of our
interest. We obtain the optimal solution to (3) with a degree-$2$ function
$\boldsymbol{\beta}(x)$ (i.e. $\beta_{i}=0$ for $i>2$) such that $\beta_{1}=0$
(due to Corollary 1). Clearly then, the problem reduces to minimizing
$\beta_{2}$ with the condition that $\boldsymbol{\beta}(u)\leq 0$ for $u$
between $d$ and $(n-d)$. Since the function
$\boldsymbol{\beta}(x)=1+\beta_{2}K_{2}(x)$ is quadratic symmetric around
$n/2$, it is clear that this is equivalent to the condition
$\boldsymbol{\beta}(d)=1+\beta_{2}K_{2}(d)\leq 0$. Thus the best degree-$2$
solution to (3) is exactly $\beta_{0}=1$, $\beta_{2}=-\frac{1}{K_{2}(d)}$, and
$\beta_{i}=0$ for all other $i$. We now prove that this is actually the
overall optimal solution to (3) (for even $d$ and $n$).

###### Lemma 2

$\beta_{0}=1$, $\beta_{2}=-\frac{1}{K_{2}(d)}$, and $\beta_{i}=0$ for all
other $i$ is an optimal solution for (3) when $\frac{n-\sqrt{n}}{2}+1\leq
d\leq n/2$ and $n$, $d$ are even.

###### Proof:

We use duality of linear programming to prove optimally. Namely, we use the
fact that if $\beta$ is feasible for the LP $\eqref{bal_LP}$, some $\alpha$ is
feasible for its dual LP (3), and complementary slackness conditions are
satisfied, then $\beta$ and $\alpha$ are optimal for the LP and its dual.

Complementary slackness conditions for our case are:

$\displaystyle\boldsymbol{\beta}(u)\cdot\alpha_{u}=0$ $\displaystyle
u=d,d+1,\dots,(n-d),$ (5)
$\displaystyle\beta_{s}\cdot\left(\sum_{k=0}^{n}\alpha_{k}K_{s}(k)\right)=0$
$\displaystyle s=1,2,\dots,n.$

Since $\boldsymbol{\beta}(x)=1-K_{2}(x)/K_{2}(d)$ only has roots at $d$ and
$(n-d)$, we immediately see that $\alpha_{u}=0$ for all $u$ other than $d$ and
$(n-d)$. Further, since $\beta_{2}\not=0$, we have that
$\binom{n}{2}+\alpha_{d}\cdot K_{2}(d)+\alpha_{n-d}\cdot K_{2}(n-d)=0$. We
claim that taking

$\alpha_{d}=\alpha_{n-d}=-\frac{\binom{n}{2}}{2K_{2}(d)}$

and $\alpha_{i}=0$ for every other coordinate gives a feasible solution
$\alpha$ to the LP (3). Notice that all the complementary slackness conditions
are satisfied for such $\beta$ and $\alpha$.

We now prove that Delsarte-MacWilliam’s inequalities in (3) are satisfied.
Specifically, we need to show

$\binom{n}{s}+\alpha_{d}K_{s}(d)+\alpha_{n-d}K_{s}(n-d)\geq 0,\qquad
s=0,1,\dots,n.$ (6)

The case $s=0$ is straightforward as $\alpha\geq 0$, and for ${s=2}$ the
equality holds by the choice of $\alpha$. Next, for odd $s$,
$K_{s}(d)=-K_{s}(n-d)$ by Fact 1, and therefore $\binom{n}{s}+\alpha_{d}\cdot
K_{s}(d)+\alpha_{n-d}\cdot K_{s}(n-d)=\binom{n}{s}\geq 0$. Since $n$ is even,
this applies to $s=(n-1)$. Moreover, for $s=n$, $K_{n}(d)=K_{0}(d)$ and
$K_{n}(n-d)=K_{0}(n-d)$ (using Fact 2 and since $n,d$ are even), so the
inequality also holds.

For every even $s$ we have $K_{s}(d)=K_{s}(n-d)$. So it remains to prove for
every even $2<s<n-1$ that

$\binom{n}{s}\geq-2\cdot\alpha_{d}\cdot
K_{s}(d)=\frac{\binom{n}{2}}{K_{2}(d)}\cdot K_{s}(d).$ (7)

Denote for convenience $C=\left\lvert K_{2}(d)/\binom{n}{2}\right\lvert$. We
prove the following statement, from which (7) clearly follows

$|K_{s}(d)|\leq C\cdot\binom{n}{s},\qquad s=2,3,\dots,(n-2).$ (8)

Notice that it is sufficient to show (8) only for $2\leq s\leq n/2$, as the
inequality for all other values of $s$ will follow from Fact 1.

We are going to need the following

###### Lemma 3

Let $|K_{q-1}(d)|\leq\delta\cdot\binom{n}{q-1}$ and
$|K_{q}(d)|\leq\delta\cdot\binom{n}{q}$ for some $\delta>0$, $d<n/2$, and
positive integer $q<n$. Then

$|K_{q+1}(d)|\leq\delta\cdot\binom{n}{q+1}\cdot\frac{n-2d+q}{n-q}.$ (9)

For clarity of exposition, we defer its proof until the end of this section.

We show (8) in two steps. First, we prove the following:

Hypothesis: For every $1\leq t\leq d/2$,

$\displaystyle|K_{2t}(d)|$ $\displaystyle\leq
C\cdot\binom{n}{2t}\cdot\prod_{i=2}^{t}\frac{n-2d+2i}{n-2i},$ (10)
$\displaystyle|K_{2t+1}(d)|$ $\displaystyle\leq
C\cdot\binom{n}{2t+1}\cdot\prod_{i=2}^{t}\frac{n-2d+2i}{n-2i},$

where empty products are treated as $1$.

Base: For $t=1$, we already know $|K_{2}(d)|=C\cdot\binom{n}{2}$. The proof of
above inequality for $|K_{3}(d)|$ is defered to Lemma 4.

Step: Denote $\eta_{t}=\prod\limits_{i=2}^{t}\frac{n-2d+2i}{n-2i}$ for
brevity. So, suppose (10) holds for some $(t-1)$, where $1\leq t\leq d/2$.
First, we derive from Lemma 3 for $q=2t-1$:

$|K_{2t}(d)|\leq
C\cdot\eta_{t-1}\cdot\binom{n}{2t}\cdot\frac{n-2d+(2t-1)}{n-(2t-1)}\leq
C\cdot\eta_{t}\cdot\binom{n}{2t},$ (11)

where we use $\frac{n-2d+(2t-1)}{n-(2t-1)}\leq\frac{n-2d+2t}{n-2t}$. Notice
also that $\frac{n-2d+2t}{n-2t}\leq 1$ for $t\leq d/2$, and so
$\eta_{t}\leq\eta_{t-1}$. Then apply Lemma 3 for $q=2t$ and
$\delta=C\cdot\eta_{t-1}$ again:

$|K_{2t+1}(d)|\leq
C\eta_{t-1}\cdot\binom{n}{2t+1}\frac{n-2d+2t}{n-2t}=C\eta_{t}\binom{n}{2t+1}.$
(12)

Therefore, (10) holds for any $1\leq t\leq d/2$. Denote
$\Phi=\eta_{d/2}=\prod\limits_{i=2}^{d/2}\frac{n-2d+2i}{n-2i}$. We now prove
the following

Hypothesis: For every $s$ such that $d\leq s\leq n/2$,

$|K_{s}(d)|\leq
C\cdot\Phi\cdot\binom{n}{s}\cdot\prod_{k=d+1}^{s-1}\frac{n-2d+k}{n-k}.$ (13)

Base: The cases $s=d,(d+1)$ follow from (10) for $t=d/2$.

Step: Denote $\mu_{s}=\prod\limits_{k=d+1}^{s-1}\frac{n-2d+k}{n-k}$, and
notice that $\mu_{s-1}\leq\mu_{s}$ for any $s$ within the range of interest.
Suppose (13) holds for $(s-2)$ and $(s-1)$, and $(d+2)\leq s\leq n/2$. Then
apply Lemma 3 for $q=(s-1)$ and $\delta=C\cdot\Phi\cdot\mu_{s-1}$ :

$|K_{s}(d)|\leq
C\Phi\cdot\mu_{s-1}\cdot\binom{n}{s}\frac{n-2d+(s-1)}{n-(s-1)}=C\cdot\Phi\cdot\mu_{s}\cdot\binom{n}{s}.$
(14)

We are finally ready to prove (8) for every $s$ between $2$ and $n/2$. For $s$
such that $2\leq s\leq(d+1)$, (10) implies $|K_{s}(d)|\leq C\cdot\eta_{\lfloor
s/2\rfloor}\cdot\binom{n}{s}$. Clearly $\eta_{\lfloor s/2\rfloor}\leq 1$ as
every term in the product is at most $1$, so (8) holds for such $s$.

Next, we have from (13) that $|K_{s}(d)|\leq
C\cdot\Phi\cdot\mu_{s}\cdot\binom{n}{s}\leq
C\cdot\Phi\cdot\mu_{n/2}\cdot\binom{n}{s}$ for every s between $(d+2)$ and
$n/2$, as $\mu_{s}$ is an increasing sequence. So it is sufficient to show
that $\Phi\cdot\mu_{n/2}\leq 1$. Denote $w=n/2-d-1$, so $w<\sqrt{n}/2$, since
$d>(n-\sqrt{n})/2$. Recall that $\Phi=\eta_{d/2}\leq\eta_{w+1}$, and so

$\displaystyle\Phi\cdot\mu_{n/2}$
$\displaystyle\leq\prod_{i=2}^{w+1}\frac{n-2d+2i}{n-2i}\cdot\prod_{v=1}^{w}\frac{n-d+v}{n-d-v}$
(15)
$\displaystyle\leq\left(\frac{n-2d+2w+2}{n-2w-2}\right)^{w}\left(\frac{n-d+w}{n-d-w}\right)^{w}$
(16)
$\displaystyle=\left(\frac{4w+4}{n-2w-2}\cdot\frac{n/2+2w+1}{n/2+1}\right)^{w}\leq
1.$ (17)

The final inequality clearly holds because $w<\sqrt{n}/2$ .

This completes the proof of (8) for the whole range of $2\leq s\leq(n-2)$, and
together with our arguments that (7) holds for $s\in\\{0,1,(n-1),n\\}$, this
means that all Delsarte-McWilliam’s inequalities from (3) are satisfied.
Therefore, we conclude that $\alpha$ is optimal for (3) and $\beta$ is optimal
for (3). ∎

###### Proof:

We use the following recurrence for Krawtchouk polynomials:

$(q+1)K_{q+1}(x)=(n-2x)K_{q}(x)-(n-q+1)K_{q-1}(x),$ (18)

for any positive integer $q<n$, where $K_{0}(x)=1$, and $K_{1}(x)=n-2x$. The
proof now follows from a simple inductive calculation. ∎

###### Lemma 4

Let $\frac{n-\sqrt{n}}{2}+1\leq d\leq n/2$. Then

$|K_{3}(d)|\leq\frac{|K_{2}(d)|}{\binom{n}{2}}\binom{n}{3}=\frac{|K_{2}(d)|\cdot(n-2)}{3}.$
(19)

###### Proof:

We know $K_{1}(d)\geq 0$ and $K_{2}(d)<0$ for such $d$. Using recurrence (18)
and this sign information, obtain

$\displaystyle|K_{3}(d)|=\frac{-(n-2d)K_{2}(d)+(n-1)(n-2d)}{3}$ (20)

We need to show that the above is at most $-\frac{K_{2}(d)\cdot(n-2)}{3}$.
Equivalently, we need to prove

$-2(d-1)K_{2}(d)\overset{?}{\geq}(n-2d)(n-1).$ (21)

Decompose
$K_{2}(d)=2\left(d-\frac{n-\sqrt{n}}{2}\right)\left(d-\frac{n+\sqrt{n}}{2}\right)$,
and using the conditions on $d$, we have $-K_{2}(d)\geq\sqrt{n}$.

Further, using $\sqrt{n}-2\geq(n-2d)$ for (21) we finally get

$\displaystyle-2(d-1)K_{2}(d)\geq(n-\sqrt{n})\sqrt{n}$
$\displaystyle>(\sqrt{n}-2)n$ (22) $\displaystyle\geq(n-2d)(n-1).\hskip
25.00003pt\qed$ (23)

This finally brings us the the proof of our main result.

###### Proof:

Observe that
$K_{2}(d)=2d^{2}-2nd-\binom{n}{2}=-\frac{1}{2}\left(n-(n-2d)^{2}\right)$.
Applying Lemma 2 obtain

$\displaystyle B_{LP}(n,d)=1-\frac{\binom{n}{2}}{K_{2}(d)}$
$\displaystyle=1+\frac{n(n-1)}{n-(n-2d)^{2}}$ (24)
$\displaystyle=\frac{4d(n-d)}{n-(n-2d)^{2}}.\hskip 65.00009pt\qed$ (25)

##  IV Lower bound on LP bound for almost-balanced codes

Consider asymptotic lower and upper bounds on $R(\delta)$

$R_{GV}(\delta)\leq R(\delta)\leq R_{MRRW}(\delta),$

where $R_{GV}(\delta)=1-H(\delta)$ is a Gilbert-Varshamov bound obtained using
a standard packing argument, which is currently the best known lower bound on
$R(\delta)$. The bound
$R_{MRRW}(\delta)=H\left(\frac{1}{2}-\sqrt{\delta(1-\delta)}\right)$ is the
first MRRW [2] bound, which is the best known upper bound for $\delta>0.273$
(the second MRRW bound is the best known for the remaining range of $\delta$).

In [10] Samorodnitsky proved an integrality gap of at most $2$ for the MRRW
bound with respect to the true LP bound:

$\displaystyle\dfrac{R_{GV}(\delta)+R_{MRRW}(\delta)}{2}\leq
R_{LP}(\delta)\leq R_{MRRW}(\delta).$ (26)

Combined with the fact that $R_{MRRW}(\delta)>R_{GV}(\delta)$ for any
$\delta\in(0,\frac{1}{2})$, the above proved that Delsarte’s linear
programming bound cannot attain the currently best known lower bound on
$R(\delta)$.

In this section we prove an analogous result for the linear programming bound
for the almost-balanced codes. Our proof is a slight modification of a proof
from [10].

###### Theorem 2

$R^{\text{bal}}_{LP}(\delta)\geq\dfrac{R_{GV}(\delta)+R_{MRRW}(\delta)}{2}$
for any ${\delta\in(0,\frac{1}{2})}$.

To obtain a lower bound, we derive a feasible solution to (3), closely
following [10].

###### Lemma 5

Let $\varepsilon=\dfrac{1}{4n}\sqrt{\dfrac{\binom{n}{\lfloor
x_{d}\rfloor}}{2^{n}\cdot\binom{n}{d}}}$, where $x_{d}$ is the first
(smallest) root of the polynomial $K_{d}(x)$, and let

  * •

$a_{0}=1$

  * •

$a_{k}=0$ for $0\leq k\leq d-1$ and ${(n-d+1)}\leq k\leq n$

  * •

$a_{d}=a_{n-d}=\varepsilon\cdot(d+1)\cdot\binom{n}{d}$

  * •

$a_{k}=\varepsilon\cdot\binom{n}{k}$ for $d<k<(n-d)$.

Then $a_{0},\dots,a_{n}$ is a feasible solution to the LP (3).

###### Proof:

Clearly, we only need to verify the first set of constraints in (3) (Delsarte-
MacWilliam’s inequalities). The case $s=0$ is immediate, so assume $s\geq 1$.
We have

$\displaystyle\sum_{k=0}^{n}a_{k}K_{s}(k)=$ $\displaystyle
K_{s}(0)+\varepsilon\cdot\sum_{k=d}^{n-d}\binom{n}{k}K_{s}(k)$ (27)
$\displaystyle\hskip 20.0pt+\varepsilon\cdot
d\binom{n}{d}\Big{(}K_{s}(d)+K_{s}(n-d)\Big{)}$ (28)
$\displaystyle=\binom{n}{s}+\varepsilon\sum_{k=d}^{n-d}\binom{n}{s}K_{k}(s)+\varepsilon
d\binom{n}{s}\Big{(}K_{d}(s)+K_{n-d}(s)\Big{)},$ (29)

where we used reciprocity property from Fact 3.

Next we use Facts 1-2, and consider two cases. When $s$ is odd, everything
except the first summand in the RHS of (29) cancels out. Indeed,
$K_{d}(s)=-K_{n-d}(s)$, and the summation in the middle can be written as

$\displaystyle\sum_{k=d}^{n-d}\binom{n}{s}K_{k}(s)=$
$\displaystyle\binom{n}{s}\sum_{k=d}^{\frac{n}{2}-1}\big{(}K_{k}(s)+K_{n-k}(s)\big{)}$
(30) $\displaystyle+\binom{n}{s}\cdot K_{n/2}(s)\cdot\big{[}(n+1)\ \text{mod}\
2\big{]},$

Observe that $K_{k}(s)+K_{n-k}(s)=0$ for any $k$ within the summation range.
Finally, for even $n$ we have
$\binom{n}{s}K_{n/2}(s)=\binom{n}{n/2}K_{s}(n/2)=0$, as $K_{s}(n/2-x)$ is an
odd function. Therefore,
$\sum\limits_{k=0}^{n}a_{k}K_{s}(k)=K_{s}(0)=\binom{n}{s}\geq 0$ for odd $s$.

Now consider the case of even $s$. Using the fact that Krawtchouk polynomials
are orthogonal with respect to the binomial measure
$\mu(k)=\binom{n}{k}/2^{n}$ and that $K_{0}(k)\equiv 1$, obtain
$2^{n}\cdot\sum_{k=0}^{n}\mu(k)K_{s}(k)\cdot
K_{0}(k)=\sum_{k=0}^{n}\binom{n}{k}K_{s}(k)=0$. Then in (27)-(28) for the
summation in the RHS we can write

$\displaystyle\sum_{k=d}^{n-d}$
$\displaystyle\binom{n}{k}K_{s}(k)=-\sum_{k=0}^{d-1}\binom{n}{k}K_{s}(k)-\sum_{k=n-d+1}^{n}\binom{n}{k}K_{s}(k)$
(31)
$\displaystyle=-\sum_{k=0}^{d-1}\binom{n}{s}K_{k}(s)-\sum_{k=n-d+1}^{n}\binom{n}{s}K_{k}(s)$
$\displaystyle=-\binom{n}{s}\sum_{k=0}^{d-1}\Big{(}K_{k}(s)+K_{n-k}(s)\Big{)}=-2\binom{n}{s}\sum_{k=0}^{d-1}K_{k}(s),$

where we used reciprocity and Fact 2 for even $s$. Using these properties
again for last part of RHS in (27)-(28), we obtain

$\sum_{k=0}^{n}a_{k}K_{s}(k)=\binom{n}{s}\left(1-2\varepsilon\sum_{k=0}^{d-1}K_{k}(s)+2\varepsilon\cdot
d\,K_{d}(s)\right).$ (32)

Finally, we notice that the RHS of the above equation is exactly the
expression derived by Samorodnitsky in [10, eq. (22)] (we took $\varepsilon$
exactly two times smaller than in [10] for these expressions to coincide),
where it was proven to be non-negative. Therefore, $a_{0},\dots,a_{n}$ is
feasible for (3). ∎

###### Proof:

Taking the feasible solution $a_{0},\dots,a_{n}$ for (3) from Lemma 5 we
obtain

$B_{LP}(n,d)\geq\sum_{k=0}^{n}a_{k}\geq\varepsilon\sum_{k=d}^{n-d}\binom{n}{k}.$

Consider some fixed $\delta\in(0,1/2)$ and $d=\lfloor\delta n\rfloor$ as $n$
increases. Standard concentration properties of Binomial distribution (e.g.
Chernoff bound) then imply that for large enough $n$, most of the weight will
lie between $\binom{n}{d}/2^{n}$ and $\binom{n}{n-d}/2^{n}$. Then for such
large $n$ we write

$B_{LP}(n,d)\geq\varepsilon\cdot
2^{n-1}\geq\dfrac{1}{8n}\sqrt{\dfrac{\binom{n}{\lfloor x_{d}\rfloor}\cdot
2^{n}}{\binom{n}{d}}}.$

Finally, we use an asymptotic for the first root of Krawtchouk polynomial
$K_{d}$ as $n$ goes to infinity:
${x_{d}=n\left(\frac{1}{2}-\sqrt{\frac{d}{n}\left(1-\frac{d}{n}\right)}\right)+o(n)}$.
Together with an asymptotic
$\lim\limits_{n\to\infty}\frac{1}{n}\log_{2}\binom{n}{\gamma n}=H(\gamma)$, we
derive

$\displaystyle
R^{\text{bal}}_{LP}(\delta)\geq\lim\limits_{n\to\infty}\frac{1}{2n}\Bigg{[}\log_{2}\binom{n}{\lfloor
x_{d}\rfloor}+n-\log_{2}\binom{n}{\delta n}\Bigg{]}$
$\displaystyle=\dfrac{1-H(\delta)+H\left(1/2-\sqrt{\delta\left(1-\delta\right)}\right)}{2}=\dfrac{R_{GV}(\delta)+R_{MRRW}(\delta)}{2}.$

∎

## Acknowledgment

Research supported in part by NSF grants CCF-1563742 and CCF-1814603.

We thank the anonymous reviewers who provided valuable comments about the
paper and brought important references to our attention.

## References

  * [1] P. Delsarte, _An Algebraic Approach to the Association Schemes of Coding Theory_ , ser. Philips journal of research / Supplement. N.V. Philips’ Gloeilampenfabrieken, 1973. [Online]. Available: https://books.google.com/books?id=zna0SgAACAAJ
  * [2] R. McEliece, E. Rodemich, H. Rumsey, and L. Welch, “New upper bounds on the rate of a code via the Delsarte-MacWilliams inequalities,” _IEEE Transactions on Information Theory_ , vol. 23, no. 2, pp. 157–166, 1977\.
  * [3] J. Naor and M. Naor, “Small-bias probability spaces: Efficient constructions and applications,” _SIAM J. Comput._ , vol. 22, no. 4, pp. 838–856, 1993\. [Online]. Available: https://doi.org/10.1137/0222053
  * [4] A. Ta-Shma, “Explicit, almost optimal, epsilon-balanced codes,” in _Proceedings of the 49th Annual ACM SIGACT Symposium on Theory of Computing_ , ser. STOC 2017. New York, NY, USA: Association for Computing Machinery, 2017, p. 238–251. [Online]. Available: https://doi.org/10.1145/3055399.3055408
  * [5] L. D. Grey, “Some bounds for error-correcting codes,” _IRE Trans. Inf. Theory_ , vol. 8, no. 3, pp. 200–202, 1962. [Online]. Available: https://doi.org/10.1109/TIT.1962.1057721
  * [6] R. A. Rankin, “The closest packing of spherical caps in n dimensions,” _Proceedings of the Glasgow Mathematical Association_ , vol. 2, no. 3, p. 139–144, 1955.
  * [7] ——, “On the minimal points of positive definite quadratic forms,” _Mathematika_ , vol. 3, no. 1, p. 15–24, 1956.
  * [8] F. MacWilliams and N. Sloane, _The Theory of Error-Correcting Codes_ , 2nd ed. North-holland Publishing Company, 1978.
  * [9] G. McGuire, “Quasi-symmetric designs and codes meeting the grey–rankin bound,” _Journal of Combinatorial Theory, Series A_ , vol. 78, no. 2, pp. 280–291, 1997. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0097316597927655
  * [10] A. Samorodnitsky, “On the optimum of Delsarte’s linear program,” _Journal of Combinatorial Theory, Series A_ , vol. 96, no. 2, pp. 261 – 287, 2001. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0097316501931760
  * [11] T. Helleseth, T. Kløve, and V. Levenshtein, “A bound for codes with given minimum and maximum distances,” 08 2006, pp. 292 – 296.