https://doi.org/10.5194/ms-15-183-2024
https://doi.org/10.5194/ms-15-183-2024
Research article |  | 19 Mar 2024

# Improved strategies of the Equality Set Projection (ESP) algorithm for computing polytope projection

Binbin Pei, Wenfeng Xu, and Yinghui Li
Abstract

​​​​​​​This paper proposes an optimization method for the Equality Set Projection algorithm to compute the orthogonal projection of polytopes. However, its computational burden significantly increases for the case of dual degeneracy, which limits the application of the algorithm. Two improvements have been proposed to solve this problem for the Equality Set Projection algorithm: first, a new criterion that does not require a discussion of the uniqueness of the solution in linear programming, which simplifies the algorithm process and reduces the computational cost; and second, an improved method that abandons the calculation of a ridge's equality set to reduce the computational burden in the case of high-dimensional dual degeneracy.

1 Introduction

Orthogonal projection is valuable for both theoretical research and engineering applications , particularly as a basic operator for polytopes, which are a type of geometric structure composed of flat boundaries and giving a good representation of linear constraints on variables. The orthogonal projection of polytopes is extensively applied in the field of control and optimization . For example, the calculations of the Minkowski sum , controllable set , reachable set , and parametric linear program can all be transformed into the projection of polytopes.

Traditional projection algorithms can be grouped into three classes, i.e., Fourier elimination , block elimination , and vertex enumeration (Avis1998). The former two algorithms require expensive computation and produce a lot of redundant inequalities. The third algorithm is suitable for polytopes with fewer vertices; however, the number of vertices of a high-dimensional polytope exceeds the number of half-spaces. Thus, it is not an efficient algorithm for the high-dimensional cases.

The Equality Set Projection (ESP) algorithm was proposed in , and it applies the inspiration of the gift-wrapping algorithm to the projection calculation. Under the condition of non-degeneracy, the complexity of ESP is linear in the number of half-spaces, which offers a tremendous advantage compared to the traditional projection algorithms. However, in the case of dual degeneracy, the computation grows exponentially.

Dual degeneracy is a frequently encountered phenomenon in control applications. For instance, when variables have constant maximum and minimum constraints, the occurrence of dual degeneracy is more likely. Unfortunately, this limitation curtails the wider application of ESP algorithms in the control field. Therefore, the primary motivation behind enhancing the ESP method is to improve its performance under dual-degeneracy conditions. We propose a novel criterion based on the non-zero rows of the null-space matrix to address the issue of discussing the uniqueness of solutions under such conditions. Furthermore, we introduce an improved strategy that eliminates the need to search for equation sets under dual-degeneracy conditions. These two advancements render the ESP algorithm more concise, efficient, and easily implementable. Consequently, the algorithm's performance in the situations involving dual degeneracy has been significantly enhanced, which paves the way for further promotion and application of ESP in engineering practices.

The paper is organized as follows. The preliminary knowledge and the workflow of the original ESP algorithm are introduced in Sect. 2. The two proposed improvement strategies are given in Sects. 3 and 4. Numerical simulation comparisons between the original ESP and the improved ESP are reported in Sect. 5, and final conclusions given in Sect. 6.

2 Preliminaries and existing results

## 2.1 Definitions and notation

• 1.

A polytope is a bounded polyhedron defined by the intersection of closed half-spaces. Let P equal $\left\{\left(\mathbit{x},\mathbit{y}\right)\in {R}^{d}×{R}^{k}|\mathbf{C}\mathbit{x}+\mathbf{D}\mathbit{y}\le \mathbit{b}\right\}$; the orthogonal projection of P onto x space is defined as

${\mathit{\pi }}_{d}\left(P\right)=\left\{\mathbit{x}\in {R}^{d}|\exists \mathbit{y}\in {R}^{k},\left(\mathbit{x},\mathbit{y}\right)\in P\right\}.$
• 2.

f is a face of the polytope $P\subset {R}^{d}×{R}^{k}$ if there exists a hyperplane , where cRd, dRk, and bR such that

$f=P\cap \left\{\left(\mathbit{x},\mathbit{y}\right)\in {R}^{d}×{R}^{k}|{\mathbit{c}}^{\mathrm{T}}\mathbit{x}+{\mathbit{d}}^{\mathrm{T}}\mathbit{y}=b\right\}$

and ${\mathbit{c}}^{\mathrm{T}}\mathbit{x}+{\mathbit{d}}^{\mathrm{T}}\mathbit{y}\le b$ for all $\left(\mathbit{x},\mathbit{y}\right)\in P$. The face of the polytope $P\subset {R}^{d}×{R}^{k}$ refers to the face of dimension n−1, whereas the ridge refers to the face of dimension n−2.

• 3.

The affine hull of every face of the polytope is

$\text{aff}f=\left\{\left(\mathbit{x},\mathbit{y}\right)\in {R}^{d}×{R}^{k}|{\mathbit{c}}^{\mathrm{T}}\mathbit{x}+{\mathbit{d}}^{\mathrm{T}}\mathbit{y}=b\right\}.$
• 4.

If A is a matrix and $E\subseteq \left\{\mathrm{1},\mathrm{\dots },q\right\}$, then AE is a matrix formed by rows of A whose indices are in E. If P is a polytope, then the notation PE refers to the set .

• 5.

E is an equality set of P if and only if E=G(E), where

$G\left(E\right)=\left\{i\in \left\{\mathrm{1},\mathrm{\dots },q\right\}|{\mathbf{C}}_{i}\mathbit{x}+{\mathbf{D}}_{i}\mathbit{y}={b}_{i},\forall \left(\mathbit{x},\mathbit{y}\right)\in {P}_{E}\right\}.$
• 6.

N(A) is defined as a matrix whose columns form an orthonormal basis for the null space of A.

## 2.2 Auxiliary results

The following results can be found in , , , and , respectively. Their proofs are omitted for the sake of brevity.

Theorem 2.1 (Jones2005) Given a polytope

$P=\left\{\left(\mathbit{x},\mathbit{y}\right)\in {R}^{d}×{R}^{k}|\mathbf{C}\mathbit{x}+\mathbf{D}\mathbit{y}\le \mathbit{b}\right\}\phantom{\rule{0.125em}{0ex}},$

if E is an equality set of P, then PE is a face of P. Furthermore, if F is a face of P, then there exists a unique equality set such that F=PE.

Proposition 2.1 (Ziegler1995) If $P\subset {R}^{d}×{R}^{k}$ is a polytope, then the projection of P onto Rd is a polytope.

Proposition 2.2 (Webster1994) If ${F}_{\mathrm{1}},\mathrm{\dots },{F}_{t}$ are facets of a polytope $P=\left\{\mathbit{z}\in {R}^{n}|A\mathbit{z}\le \mathbit{b}\right\}$, $\text{aff}\phantom{\rule{0.125em}{0ex}}{F}_{i}=\left\{\mathbit{z}\in {R}^{n}|{\mathbit{a}}_{i}^{\mathrm{T}}\mathbit{z}={b}_{i}\right\}$, where each aiRn, biR, then

Proposition 2.3 (Ziegler1995) If $P\subset {R}^{d}×{R}^{k}$ is a polytope, then for every face F of πd(P), the pre-image ${\mathit{\pi }}_{d}^{-\mathrm{1}}\left(F\right)$ is a face of P.

Proposition 2.4 (Jones2005) If $M=\left\{\left(\mathbit{x},\mathbit{y}\right)\in {R}^{d}×{R}^{k}|{\mathbf{C}}_{E}\mathbit{x}+{\mathbf{D}}_{E}\mathbit{y}={b}_{E}\right\}$ is a affine set, then the projection of M onto Rd is

Proposition 2.5 If E is an equality set of polytope

$P=\left\{\left(\mathbit{x},\mathbit{y}\right)\in {R}^{d}×{R}^{k}|\mathbf{C}\mathbit{x}+\mathbf{D}\mathbit{y}\le \mathbit{b}\right\}\phantom{\rule{0.125em}{0ex}},$

then $\text{dim}\phantom{\rule{0.25em}{0ex}}{\mathit{\pi }}_{d}\left({P}_{E}\right)=\text{dim}\phantom{\rule{0.25em}{0ex}}{P}_{E}-k+\text{rank}\phantom{\rule{0.25em}{0ex}}{\mathbf{D}}_{E}$.

## 2.3 Original ESP algorithm

We introduce the workflow of the original ESP algorithm in this subsection. A brief introduction is given here, and the details can be found in .

The algorithm is initialized by discovering a random facet Fπ in πd(P). Then we compute all ridges contained in the initial facet and search all facets adjacent to it through each ridge. Subsequently, we compute all the ridges of each new facet and let the iterating process continue until all the ridges appear twice. According to the property of the polytope, each ridge exactly connects to two facets. Therefore, all facets of πd(P) can be found by this iteration. From Proposition 2.2, we finally get .

The ESP algorithm consists of three key oracles as given below. Shooting Oracle:

$\left({E}_{{f}_{\mathrm{0}}},{\mathbit{a}}_{f},{b}_{f}\right)=\text{SHOOT}\left(P\right).$

Ridge Oracle:

$L{E}_{r}=\text{RDG}\left({E}_{f},{\mathbit{a}}_{f},{b}_{f}\right).$

$\left({E}_{\mathrm{adj}},{\mathbit{a}}_{\mathrm{adj}},{b}_{\mathrm{adj}}\right)=\text{ADJ}\left(\left({E}_{f},{\mathbit{a}}_{f},{b}_{f}\right)\left({E}_{r},{\mathbit{a}}_{r},{b}_{r}\right)\right).$

The SHOOT oracle randomly obtains an equality set ${E}_{{f}_{\mathrm{0}}}$ of P satisfying ${\mathit{\pi }}_{d}\left({P}_{{E}_{{f}_{\mathrm{0}}}}\right)=\left\{\mathbit{x}|{\mathbit{a}}_{f}^{\mathrm{T}}\mathbit{x}={b}_{f}\right\}\cap {\mathit{\pi }}_{d}\left(P\right)$ as a facet of πd(P), the RDG oracle returns all ridges contained in the facet $\left({E}_{f},{\mathbit{a}}_{f},{b}_{f}\right)$ of πd(P), and the ADJ oracle takes two equality sets Er and Ef of P and returns the unique equality set Eadj, where ${\mathit{\pi }}_{d}\left({P}_{{E}_{f}}\right)=\left\{\mathbit{x}|{\mathbit{a}}_{f}^{\mathrm{T}}\mathbit{x}={b}_{f}\right\}\cap {\mathit{\pi }}_{d}\left(P\right)$ is a facet of πd(P) and ${\mathit{\pi }}_{d}\left({P}_{{E}_{r}}\right)=\left\{\mathbit{x}|{\mathbit{a}}_{r}^{\mathrm{T}}\mathbit{x}={b}_{r}\right\}\cap {\mathit{\pi }}_{d}\left({P}_{{E}_{f}}\right)$ is a ridge of πd(P) with ${\mathit{\pi }}_{d}\left({P}_{{E}_{r}}\right)\subset {\mathit{\pi }}_{d}\left({P}_{{E}_{f}}\right)$.

### 2.3.1 Outline of the ESP algorithm

Step 1. Get an initial facet $\left({E}_{{f}_{\mathrm{0}}},{\mathbit{a}}_{f},{b}_{f}\right)$ of P: $\left({E}_{{f}_{\mathrm{0}}},{\mathbit{a}}_{f},{b}_{f}\right)=\text{SHOOT}\left(P\right)$. Then, all the ridges contained in the facet are calculated: $L{E}_{r}=\text{RDG}\left({E}_{{f}_{\mathrm{0}}},{\mathbit{a}}_{f},{b}_{f}\right)$. For each element $\left({E}_{{r}_{j}},{a}_{{r}_{j}},{b}_{{r}_{j}}\right)$ in LEr, record $\left(\left({E}_{f\mathrm{0}},{\mathbit{a}}_{f\mathrm{0}},{b}_{f\mathrm{0}}\right),\left({E}_{{r}_{j}},{a}_{{r}_{j}},{b}_{{r}_{j}}\right)\right)$ in list L and let $\mathbf{G}={\mathbit{a}}_{f}^{\mathrm{T}}$, g=bf.

Step 2. If Lϕ, select an element $\left(\left({E}_{{f}_{i}},{a}_{{f}_{i}},{b}_{{f}_{i}}\right),\left({E}_{{r}_{i}},{a}_{{r}_{i}},{b}_{{r}_{i}}\right)\right)$ randomly from list L and let $\left(\left({E}_{\mathrm{fn}},{a}_{\mathrm{fn}},{b}_{\mathrm{fn}}\right),\left({E}_{\mathrm{rn}},{a}_{\mathrm{rn}},{b}_{\mathrm{rn}}\right)\right)=\left(\left({E}_{{f}_{i}},{a}_{{f}_{i}},{b}_{{f}_{i}}\right),\left({E}_{{r}_{i}},{a}_{{r}_{i}},{b}_{{r}_{i}}\right)\right)$. If L=ϕ, stop and report G and g.

Step 3. Search the adjacent facet such that $\left({E}_{\mathrm{adj}},{\mathbit{a}}_{\mathrm{adj}},{b}_{\mathrm{adj}}\right)=\text{ADJ}\left(\left({E}_{\mathrm{fn}},{a}_{\mathrm{fn}},{b}_{\mathrm{fn}}\right),\left({E}_{\mathrm{rn}},{a}_{\mathrm{rn}},{b}_{\mathrm{rn}}\right)\right)$ and let

$\mathbf{G}=\left[\begin{array}{c}\mathbf{G}\\ {\mathbit{a}}_{\mathrm{adj}}^{\mathrm{T}}\end{array}\right],\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\mathbit{a}=\left[\begin{array}{c}\mathbit{a}\\ {b}_{\mathrm{adj}}\end{array}\right].$

Then, all the ridges contained in the adjacent facet can be calculated: $L{E}_{r}=\text{RDG}\left({E}_{\mathrm{adj}},{\mathbit{a}}_{\mathrm{adj}},{b}_{\mathrm{adj}}\right)$. For each element $\left({E}_{{r}_{j}},{a}_{{r}_{j}},{b}_{{r}_{j}}\right)$ in LEr, if there exists an element in L that already contains ${E}_{{r}_{j}}$, we remove the element from the list. Otherwise, we add element $\left(\left({E}_{{r}_{j}},{a}_{{r}_{j}},{b}_{{r}_{j}}\right),\left({E}_{\mathrm{adj}},{\mathbit{a}}_{\mathrm{adj}},{b}_{\mathrm{adj}}\right)\right)$ to list L. Go to Step 2.

### 2.3.2 SHOOT oracle

Step 1 (selecting a direction γ and finding a point on the surface of P)

Randomly choose γRd and solve the following linear programming.

$\begin{array}{}\text{(1)}& \begin{array}{l}\left({r}^{*},{\mathbit{y}}^{*}\right)=\underset{\left(r,\mathbit{y}\right)\in R×{R}^{k}}{\mathrm{arg}\phantom{\rule{0.125em}{0ex}}max}r,\\ \text{s.t.}\phantom{\rule{0.25em}{0ex}}\mathbf{C}\mathbit{\gamma }r+\mathbf{D}\mathbit{y}\le \mathbit{b}\end{array}\end{array}$

Step 2 (determining whether a facet of πd(P) is found)

Get ${E}_{{f}_{\mathrm{0}}}$ as ${E}_{{f}_{\mathrm{0}}}=\left\{i|{\mathbf{C}}_{i}\mathbit{\gamma }r+{\mathbf{D}}_{i}\mathbit{y}={b}_{i}\right\}$. If $\text{Rank}\left(\mathbf{N}{\left({\mathbf{D}}_{{E}_{{f}_{\mathrm{0}}}}^{\mathrm{T}}\right)}^{\mathrm{T}}{\mathbf{C}}_{{E}_{{f}_{\mathrm{0}}}}\right)=\mathrm{1}$, go to Step 3. Otherwise, go to Step 1.

Step 3 (computing the affine hull of the facet and normalizing)

$\begin{array}{l}\left[{\mathbit{a}}_{f}^{\mathrm{T}},{b}_{f}\right]=\mathbf{N}{\left({\mathbf{D}}_{{E}_{{f}_{\mathrm{0}}}}^{\mathrm{T}}\right)}^{\mathrm{T}}\left[{\mathbf{C}}_{{E}_{{f}_{\mathrm{0}}}},{b}_{{E}_{{f}_{\mathrm{0}}}}\right]\\ \left[{\mathbit{a}}_{f}^{\mathrm{T}},{b}_{f}\right]={\frac{\text{sign}\phantom{\rule{0.25em}{0ex}}{b}_{f}}{∥{\mathbit{a}}_{f}∥}}_{\mathrm{2}}\left[{\mathbit{a}}_{f}^{\mathrm{T}},{b}_{f}\right]\end{array}$

Step 4 (handling the dual degeneracy in Step 1)

If the linear programming (1) in Step 1 is not dual-degenerate, i.e., (1) has a unique optimizer, then stop. Otherwise, ${E}_{{f}_{\mathrm{0}}}$ needs to be recalculated as follows.

For each row ${\mathbf{C}}_{i}\mathbit{x}+{\mathbf{D}}_{i}\mathbit{y}\le {b}_{i}$ of P, we compute the following linear programming problem:

${J}^{*}\left(i\right)=\underset{\mathbit{x},\mathbit{y}}{min}{\mathbf{C}}_{i}\mathbit{x}+{\mathbf{D}}_{i}\mathbit{y}-{b}_{i},$

$\begin{array}{rl}\text{s.t.}& \phantom{\rule{0.25em}{0ex}}{\mathbf{C}}_{{E}_{c}^{i}}\mathbit{x}+{\mathbf{D}}_{{E}_{c}^{i}}\mathbit{y}\le b,\\ & \mathbit{x}={\mathbit{x}}^{*}\phantom{\rule{0.125em}{0ex}},\end{array}$

where

${E}_{c}^{i}=\left\{\mathrm{1},\mathrm{2},\mathrm{\dots },q\right\}\i,{E}_{{f}_{\mathrm{0}}}=\left\{i\in \left\{\mathrm{1},\mathrm{2},\mathrm{\dots },q\right\}|{J}^{*}\left(i\right)=\mathrm{0}\right\}.$

### 2.3.3 RDG oracle

Step 1 (computing the dimension of ${P}_{{E}_{f}}$)

If $k+d-\text{rank}\left[{\mathbf{C}}_{{E}_{f}}\phantom{\rule{0.25em}{0ex}}{\mathbf{D}}_{{E}_{f}}\right]=d-\mathrm{1}$, continue. Otherwise, go to step 3.

Step 2 ($\text{dim}{P}_{{E}_{f}}=d-\mathrm{1}$) For all $i\in {E}_{f}^{c}=\left\{\mathrm{1},\mathrm{2},\mathrm{\dots },q\right\}\{E}_{f}$,

and, solving the following linear programming,

where c is an arbitrary positive constant, ${Q}_{c}\left(i\right)={E}_{c}\Q\left(i\right)$. If τ(i)<0, then ${\mathbit{a}}_{r}={\mathbf{S}}_{\left\{i\right\}},{b}_{r}={t}_{\left\{i\right\}}$ and add $\left(Q\left(i\right),{a}_{r},{b}_{r}\right)$ to LEr.

Step 3 ($\text{dim}{P}_{{E}_{f}}>d-\mathrm{1}$)

Conduct singular value decomposition ${\mathbit{a}}_{f}^{\mathrm{T}}=\left[\begin{array}{cccc}\mathit{\sigma }& \mathrm{0}& \mathrm{\dots }& \mathrm{0}\end{array}\right]{\mathbit{V}}^{\mathrm{T}}$, $\mathbit{V}=\left[\begin{array}{cc}\stackrel{\mathrm{^}}{\mathbit{V}}& \stackrel{\mathrm{̃}}{\mathbit{V}}\end{array}\right]$. Then, we get

$\stackrel{\mathrm{̃}}{P}=\left\{\left(\stackrel{\mathrm{̃}}{\mathbit{x}},\stackrel{\mathrm{̃}}{\mathbit{y}}\right)\in {R}^{d-\mathrm{1}}×{R}^{n}|{\mathbf{S}}_{f}\stackrel{\mathrm{̃}}{\mathbit{V}}\stackrel{\mathrm{̃}}{\mathbit{x}}+{\mathbit{L}}_{f}\stackrel{\mathrm{̃}}{\mathbit{y}}\le {t}_{f}-\mathbf{S}\stackrel{\mathrm{^}}{\mathbit{V}}{b}_{f}/\mathit{\sigma }\right\},$

where ${\mathbf{S}}_{f}={\mathbf{C}}_{{E}_{f}^{c}}-{\mathbf{D}}_{{E}_{f}^{c}}{\mathbf{D}}_{{E}_{f}}^{\mathit{†}}{\mathbf{C}}_{{E}_{f}}$, ${\mathbit{L}}_{f}={\mathbf{D}}_{{E}_{f}^{c}}\mathbf{N}\left({\mathbf{D}}_{E}\right)$, and ${t}_{f}={b}_{{E}_{f}^{c}}-{\mathbf{D}}_{{E}_{f}^{c}}{\mathbf{D}}_{{E}_{f}}^{\mathit{†}}{b}_{{E}_{f}}$. The ${\mathit{\pi }}_{d-\mathrm{1}}\left(\stackrel{\mathrm{̃}}{P}\right)$ can be computed using a recursive call to the ESP algorithm. Each facet of ${\mathit{\pi }}_{d-\mathrm{1}}\left(\stackrel{\mathrm{̃}}{P}\right)$ can be converted into the ridge of πd(P) as

$\left\{\mathbit{x}|\stackrel{\mathrm{̃}}{\mathbit{V}}{\stackrel{\mathrm{̃}}{\mathbit{a}}}_{f}^{\mathrm{T}}\mathbit{x}={\stackrel{\mathrm{̃}}{b}}_{f}\right\}\cap {\mathit{\pi }}_{d-\mathrm{1}}\left(\stackrel{\mathrm{̃}}{P}\right).$

Then, we get ${\mathbit{a}}_{r}={\stackrel{\mathrm{̃}}{\mathbit{a}}}_{f}$, ${b}_{r}={\stackrel{\mathrm{̃}}{b}}_{f}$, and $\left(Q\left(i\right),{\mathbit{a}}_{r},{b}_{r}\right)$ to LEr.

Step 4 (normalizing and orthogonalizing)

$\begin{array}{l}\left[\begin{array}{cc}{\mathbit{a}}_{r}^{\mathrm{T}}& {b}_{r}\end{array}\right]←\left[\begin{array}{cc}{\mathbit{a}}_{r}^{\mathrm{T}}& {b}_{r}\end{array}\right]-{\mathbit{a}}_{f}^{\mathrm{T}}{\mathbit{a}}_{r}\left[\begin{array}{cc}{\mathbit{a}}_{f}^{\mathrm{T}}& {b}_{f}\end{array}\right]\\ \left[\begin{array}{cc}{\mathbit{a}}_{r}^{\mathrm{T}}& {b}_{r}\end{array}\right]←\frac{\text{sign}\phantom{\rule{0.125em}{0ex}}{b}_{r}}{\mathit{ł}|{\mathbit{a}}_{r}{‖}_{\mathrm{2}}}\left[\begin{array}{cc}{\mathbit{a}}_{r}^{\mathrm{T}}& {b}_{r}\end{array}\right]\end{array}\phantom{\rule{0.125em}{0ex}}.$

Step 1 (finding a point on the adjacent facet ${\mathit{\pi }}_{d}\left({P}_{{E}_{\mathrm{adj}}}\right)$ and its corresponding pre-image)

Solve the following linear programming:

$\begin{array}{}\text{(2)}& \begin{array}{l}\left({\mathbit{x}}^{*},{\mathbit{y}}^{*}\right)=\underset{\mathbit{x},\mathbit{y}}{\mathrm{arg}\phantom{\rule{0.25em}{0ex}}max}\phantom{\rule{0.25em}{0ex}}{\mathbit{a}}_{r}^{\mathrm{T}}\mathbit{x}\phantom{\rule{0.125em}{0ex}},\\ \begin{array}{rl}\text{s.t.}& \phantom{\rule{0.25em}{0ex}}{\mathbf{C}}_{{E}_{r}}\mathbit{x}+{\mathbf{D}}_{{E}_{r}}\mathbit{y}\le {b}_{{E}_{r}}\phantom{\rule{0.125em}{0ex}},\\ & {\mathbit{a}}_{f}^{\mathrm{T}}\mathbit{x}={b}_{f}\left(\mathrm{1}-\mathit{\epsilon }\right),\end{array}\end{array}\end{array}$

where ε is a sufficiently small positive constant.

If LP (2) is not dual-degenerate, i.e., LP (2) has a unique optimizer $\left({\mathbit{x}}^{*},{\mathbit{y}}^{*}\right)$, then continue. Otherwise, go to Step 3.
Step 2 (the case without dual degeneracy)

Go to Step 4.
Step 3 (the case of dual degeneracy)

$\mathit{\gamma }=-\frac{{\mathbit{a}}_{r}^{\mathrm{T}}{\mathbit{x}}^{*}-{b}_{r}}{{\mathbit{a}}_{f}^{\mathrm{T}}{\mathbit{x}}^{*}-{b}_{f}}$

For each row ${\mathbf{C}}_{i}\mathbit{x}+{\mathbf{D}}_{i}\mathbit{y}\le {b}_{i}$ of P, compute the following linear programming problem:

$\begin{array}{l}{J}^{*}\left(i\right)=\underset{\mathbit{x},\mathbit{y}}{min}{\mathbf{C}}_{i}\mathbit{x}+{\mathbf{D}}_{i}\mathbit{y}-{b}_{i}\phantom{\rule{0.125em}{0ex}},\\ \begin{array}{rl}\text{s.t.}& \phantom{\rule{0.25em}{0ex}}{\mathbf{C}}_{{E}_{c}^{i}}\mathbit{x}+{\mathbf{D}}_{{E}_{c}^{i}}\mathbit{y}\le {b}_{{E}_{c}^{i}},\\ & {\left({\mathbit{a}}_{r}+\mathit{\gamma }{\mathbit{a}}_{f}\right)}^{\mathrm{T}}\mathbit{x}={b}_{r}+\mathit{\gamma }{b}_{f},\end{array}\end{array}$

where ${E}_{c}^{i}=\left\{\mathrm{1},\mathrm{2},\mathrm{\dots },q\right\}\\left\{i\right\}$, ${E}_{\mathrm{adj}}=\left\{i|{J}^{*}\left(i\right)=\mathrm{0}\right\}$.
Step 4 (computing the affine hull of the facet and normalizing)

$\begin{array}{l}\left[{\mathbit{a}}_{f}^{\mathrm{T}},{b}_{f}\right]=\mathbf{N}{\left({\mathbf{D}}_{{E}_{\mathrm{0}}}^{\mathrm{T}}\right)}^{\mathrm{T}}\left[{\mathbf{C}}_{{E}_{\mathrm{0}}},{b}_{{E}_{\mathrm{0}}}\right]\\ \left[{\mathbit{a}}_{f}^{\mathrm{T}},{b}_{f}\right]=\frac{\text{sign}\phantom{\rule{0.25em}{0ex}}{b}_{f}}{{∥{\mathbit{a}}_{f}∥}_{\mathrm{2}}}\left[{\mathbit{a}}_{f}^{\mathrm{T}},{b}_{f}\right]\end{array}$
3 Improvement of the uniqueness discussion

In the original ESP algorithm, expressions (1) and (2) in the SHOOT and ADJ oracles need to judge the uniqueness of the optimizer of linear programming (LP), and different steps are taken according to whether the optimizer is unique or not.

In this section, an illustrative example is given at first, the barriers of the original ESP algorithm in the case of dual degeneracy are listed, and the alternative criteria are proposed.

## 3.1 An illustrative example

The linear programming (2) in Step 1 of the ADJ oracle is given for illustration, and the solution of (2) satisfies

The projection of $\left({\mathbit{x}}^{*},{\mathbit{y}}^{*}\right)$ is on ${\mathit{\pi }}_{d}\left({P}_{{E}_{\mathrm{adj}}}\right)$ (Jones2005). If the optimizer $\left({\mathbit{x}}^{*},{\mathbit{y}}^{*}\right)$ of linear programming is unique, Eadj should be taken as the equality set corresponding to the face ${P}_{{E}_{\mathrm{adj}}}$ where $\left({\mathbit{x}}^{*},{\mathbit{y}}^{*}\right)$ is located. If the optimizer is not unique, Eadj should be taken as the equality set corresponding to the face ${P}_{{E}_{\mathrm{adj}}},$ where the entire optimal solution space is located. Next, a simple example is provided to illustrate the reasons for the discussion of the uniqueness of the optimizer and leads to its alternative strategy.
Example 1

Given a polytope $P=\left\{\left(\mathbit{x},\mathbit{y}\right)\in {R}^{d}×{R}^{k}|\mathbf{C}\mathbit{x}+\mathbf{D}\mathbit{y}\le \mathbit{b}\right\}$,

$\mathbf{C}=\left[\begin{array}{cc}\mathrm{0}& -\mathrm{2}\\ \mathrm{0}& \mathrm{0.667}\\ \mathrm{0.4}& \mathrm{0}\\ -\mathrm{2}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}\end{array}\right],\phantom{\rule{0.25em}{0ex}}\mathbf{D}=\left[\begin{array}{c}\mathrm{2}\\ -\mathrm{0.667}\\ \mathrm{0.4}\\ \mathrm{0}\\ -\mathrm{1}\end{array}\right],\phantom{\rule{0.25em}{0ex}}\mathbit{b}=\left[\begin{array}{c}\mathrm{1}\\ \mathrm{1}\\ \mathrm{1}\\ -\mathrm{1}\\ -\mathrm{1}\end{array}\right],$

the orthogonal projection ${\mathit{\pi }}_{d}\left(P\right)=\left\{\mathbf{G}\mathbit{x}\le \mathbit{g}\right\}$ onto x space Rd is shown in Fig. 1, where

$\mathbf{G}=\left[\begin{array}{cc}-\mathrm{1}& \mathrm{0}\\ \mathrm{0.707}& \mathrm{0.707}\\ \mathrm{1}& \mathrm{0}\\ \mathrm{0}& -\mathrm{1}\end{array}\right],\phantom{\rule{0.25em}{0ex}}\mathbit{g}=\left[\begin{array}{c}-\mathrm{0.5}\\ \mathrm{2.83}\\ \mathrm{1.5}\\ -\mathrm{0.5}\end{array}\right].$

Figure 1Polytope and its projection onto Rd.

Figure 2 shows that the ADJ oracle takes the facet $\left\{{\mathbit{a}}_{f}^{\mathrm{T}}\mathbit{x}={b}_{f}\right\}\cap \mathit{\pi }\left(P\right)$ of P (the blue line segment in this example) and ridge $\left\{{\mathbit{a}}_{r}^{\mathrm{T}}\mathbit{x}={b}_{r}\right\}\cap {\mathit{\pi }}_{d}\left({P}_{{E}_{f}}\right)$ (the intersection of the blue line segment and the yellow line segment in this example) to solve the adjacent facet (the yellow line segment in this example is the facet to be solved), and x* is the intersection of plane ${\mathbit{a}}_{f}^{\mathrm{T}}\mathbit{x}={b}_{f}\left(\mathrm{1}-\mathit{\epsilon }\right)$ and the yellow line segment.

Figure 2The ADJ oracle in Example 1.

All the points on the red line segment of P can be projected onto x*, which indicates that all points on the red line segment are optimal solutions of the LP problem. The optimal solution domain belongs to the yellow facet of ${P}_{{E}_{\mathrm{adj}}}$, and its corresponding equality set is Eadj={4}.

However, most of the LP solvers will only give one optimal solution for the LP problem with multiple optimal solutions. If this solution is mistakenly considered the only optimal solution, this may lead to the iteration cycle and increase the computational burden.

In this case, if $\left({x}^{*},{y}^{*}\right)$ on the intersection of the blue facet and yellow facet is taken as the unique optimal solution, the ${P}_{{{E}^{\prime }}_{\mathrm{adj}}}$ that $\left({x}^{*},{y}^{*}\right)$ belongs to is the intersection line of the blue and yellow facets of P, ${E}_{\mathrm{adj}}^{\prime }=\left\{\mathrm{1},\mathrm{4}\right\}$. According to Fig. 3, ${P}_{{{E}^{\prime }}_{\mathrm{adj}}}\subset {P}_{{E}_{\mathrm{adj}}}$ and ${\mathit{\pi }}_{d}\left({P}_{{{E}^{\prime }}_{\mathrm{adj}}}\right)\subset {\mathit{\pi }}_{d}\left({P}_{{E}_{\mathrm{adj}}}\right)$.

Figure 3The non-uniqueness of the optimal solution is ignored.

Owing to $\text{aff}{\mathit{\pi }}_{d}\left({P}_{{{E}^{\prime }}_{\mathrm{adj}}}\right)=\text{aff}{\mathit{\pi }}_{d}\left({P}_{{E}_{\mathrm{adj}}}\right)$, some facets may appear multiple times in the following iterative process, which may increase the number of iterations, and even the algorithm cannot be terminated.

## 3.2 Shortcomings of the original algorithm

The barriers of the original ESP algorithm in the case of dual degeneracy can be listed as follows.

• i.

Currently, most of the LP solvers do not have the ability to judge the uniqueness of the solutions of linear programming. Although several achievements have been made on the unique conditions of LP solutions , these methods are difficult to embed into the existing LP solvers, since they induce additional computation.

• ii.

Since the dual degeneracy is only a necessary but insufficient condition for linear programming to have multiple optimizers , it is conservative to determine the uniqueness of the solution by testing whether the dual degeneracy occurs in .

• iii.

Expression (3) implies that, in the case of d≥3, the value of x* is always not unique; i.e., the LP problem always has multiple optimal solutions. However, if only the value of x* is not unique and the value of y* is unique, irrespective of the point in the optimal solution space, it can be regarded as the unique optimal solution, and the obtained ${P}_{{E}_{\mathrm{adj}}}$ does not vary. In this case, an extra computation is added to handle the non-uniqueness of the optimal solution of the LP.

## 3.3 Alternative criteria

According to Example (1), there may be multiple faces ${P}_{{E}_{f}}$ in P, and the projections of their affine hulls are all aff Fπ. The collection of all Ef satisfying this condition can be denoted as

where

$\begin{array}{rl}& G\left({E}_{f}\right)=\mathit{\left\{}i\in \left\{\mathrm{1},\mathrm{2},\mathrm{\dots },q\right\}|{\mathbf{C}}_{i}\mathbit{x}+{\mathbf{D}}_{i}\mathbit{y}={b}_{i},\\ & \forall \left(\mathbit{x},\mathbit{y}\right)\in {P}_{{E}_{f}}\mathit{\right\}}.\end{array}$

The discussion on the uniqueness of optimizers of LP is to ensure that the obtained Ef is the element in ${\mathrm{\Xi }}_{{F}_{\mathit{\pi }}}$ with the maximal $\text{dim}{P}_{{E}_{f}}$. The primary result can be stated as given below.
Theorem 3.1 Consider the ${\mathrm{\Xi }}_{{F}_{\mathit{\pi }}}$ defined in Expression (4), and Ef is the subset of ${\mathrm{\Xi }}_{{F}_{\mathit{\pi }}}$ that maximizes the value of $\text{dim}{P}_{{E}_{f}}$ if and only if there are no zero rows in $\mathbf{N}\left({\mathbf{D}}_{{E}_{f}}^{\mathrm{T}}\right)$.
Proof (sufficient part): the sufficiency of Theorem 3.1 has been proven by the following inverse negative proposition.

If there exist zero rows in $\mathbf{N}\left({\mathbf{D}}_{{E}_{f}}^{\mathrm{T}}\right)$, then Ef is not the element in ${\mathrm{\Xi }}_{{F}_{\mathit{\pi }}}$ that maximizes the value of $\text{dim}{P}_{{E}_{f}}$.

Let ${E}_{f}=\left\{{e}_{\mathrm{1}},{e}_{\mathrm{2}},\mathrm{\dots },{e}_{n},\right\}$, ${\mathbf{D}}_{E}^{\mathrm{T}}=\left[\begin{array}{cccc}{\mathbf{D}}_{{e}_{\mathrm{1}}}^{\mathrm{T}}& {\mathbf{D}}_{{e}_{\mathrm{2}}}^{\mathrm{T}}& \mathrm{\dots }& {\mathbf{D}}_{{e}_{n}}^{\mathrm{T}}\end{array}\right]$; without loss of generality, we can assume that the last m rows of $\mathbf{N}\left({\mathbf{D}}_{{E}_{f}}^{\mathrm{T}}\right)$ are zero rows.

$\begin{array}{rl}& \mathbf{N}\left({\mathbf{D}}_{{E}_{f}}^{\mathrm{T}}\right)=\left[\begin{array}{cccc}{N}_{\mathrm{1}}& {N}_{\mathrm{2}}& \mathrm{\dots }& {N}_{n-r}\end{array}\right]\\ & =\left[\begin{array}{cccc}{N}_{\mathrm{11}}& {N}_{\mathrm{12}}& \mathrm{\cdots }& {N}_{\mathrm{1}\left(n-r\right)}\\ {N}_{\mathrm{21}}& {N}_{\mathrm{22}}& \mathrm{\cdots }& {N}_{\mathrm{2}\left(n-r\right)}\\ \mathrm{⋮}& \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ {N}_{\left(n-m\right)\mathrm{1}}& {N}_{\left(n-m\right)\mathrm{2}}& \mathrm{\cdots }& {N}_{\left(n-m\right)\left(n-r\right)}\\ \mathrm{0}& \mathrm{0}& \mathrm{\cdots }& \mathrm{0}\\ \mathrm{⋮}& \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ \mathrm{0}& \mathrm{0}& \mathrm{\cdots }& \mathrm{0}\end{array}\right]\\ & =\left[\begin{array}{c}{N}_{\mathrm{1}}^{\mathrm{row}}\\ {N}_{\mathrm{2}}^{\mathrm{row}}\\ \mathrm{⋮}\\ {N}_{n-m}^{\mathrm{row}}\\ {N}_{n-m+\mathrm{1}}^{\mathrm{row}}\\ \mathrm{⋮}\\ {N}_{n}^{\mathrm{row}}\end{array}\right]\end{array}.$

According to ${\mathbf{D}}_{{E}_{f}}^{\mathrm{T}}\mathbf{N}\left({\mathbf{D}}_{{E}_{f}}^{\mathrm{T}}\right)=\sum _{i=\mathrm{1}}^{n}{\mathbf{D}}_{{e}_{i}}^{\mathrm{T}}{N}_{i}^{\mathrm{row}}=\mathrm{0}$, there is a one-to-one correspondence between ${N}_{i}^{\mathrm{row}}$ and ${\mathbf{D}}_{ei}^{\mathrm{T}}$; i.e., there is a one-to-one correspondence between the rows of $\mathbf{N}\left({\mathbf{D}}_{{E}_{f}}^{\mathrm{T}}\right)$ and the elements of Ef.

Let ${E}_{f}^{\prime }=\left\{{e}_{\mathrm{1}},{e}_{\mathrm{2}},\mathrm{\dots },{e}_{n-m}\right\}$.

$\begin{array}{rl}& \mathbf{N}\left({\mathbf{D}}_{{E}_{f}}^{\mathrm{T}}\right)=\left[\begin{array}{cccc}{N}_{\mathrm{1}}& {N}_{\mathrm{2}}& \mathrm{\cdots }& {N}_{n-r}\end{array}\right]\\ & =\left[\begin{array}{cccc}{N}_{\mathrm{11}}& {N}_{\mathrm{12}}& \mathrm{\cdots }& {N}_{\mathrm{1}\left(n-r\right)}\\ {N}_{\mathrm{21}}& {N}_{\mathrm{22}}& \mathrm{\cdots }& {N}_{\mathrm{2}\left(n-r\right)}\\ \mathrm{⋮}& \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ {N}_{\left(n-m\right)\mathrm{1}}& {N}_{\left(n-m\right)\mathrm{2}}& \mathrm{\cdots }& {N}_{\left(n-m\right)\left(n-r\right)}\\ \mathrm{0}& \mathrm{0}& \mathrm{\cdots }& \mathrm{0}\\ \mathrm{⋮}& \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ \mathrm{0}& \mathrm{0}& \mathrm{\cdots }& \mathrm{0}\end{array}\right]\\ & =\left[\begin{array}{c}{N}_{\mathrm{1}}^{\mathrm{row}}\\ {N}_{\mathrm{2}}^{\mathrm{row}}\\ \mathrm{⋮}\\ {N}_{n-m}^{\mathrm{row}}\\ {N}_{n-m+\mathrm{1}}^{\mathrm{row}}\\ \mathrm{⋮}\\ {N}_{n}^{\mathrm{row}}\end{array}\right]\end{array}$

$\left[\begin{array}{cccc}{N}_{\mathrm{1}}& {N}_{\mathrm{2}}& \mathrm{\dots }& {N}_{n-r}\end{array}\right]$ forms a basis for the null space of ${\mathbf{D}}_{{E}_{f}}^{\mathrm{T}}$, implying that $\left[\begin{array}{cccc}{N}_{\mathrm{1}}^{\prime }& {N}_{\mathrm{2}}^{\prime }& \mathrm{\dots }& {N}_{n-r}^{\prime }\end{array}\right]$ forms a basis for the null space of ${\mathbf{D}}_{{E}_{f}^{\prime }}^{\mathrm{T}}$.

Next, we will prove that ${E}_{f}^{\prime }$ is an equality set, i.e.,

${E}_{f}^{\prime }=\left\{i\in \left\{\mathrm{1},\mathrm{\cdots },q\right\}|{\mathbf{C}}_{i}\mathbit{x}+{\mathbf{D}}_{i}\mathbit{y}={b}_{i},\forall \left(\mathbit{x},\mathbit{y}\right)\in {P}_{{E}_{f}^{\prime }}\right\}.$

Since Ef is an equation set, we only need to prove that

$\forall {e}_{i}\in {E}_{f}^{\prime \prime },\exists \left(\mathbit{x},\mathbit{y}\right)\in {P}_{{E}_{f}^{\prime }},\text{s.t.}\phantom{\rule{0.25em}{0ex}}{\mathbf{C}}_{i}\mathbit{x}+{\mathbf{D}}_{i}\mathbit{y}\ne {b}_{i},$

where

${E}_{f}^{\prime \prime }=\left\{{e}_{n-m+\mathrm{1}},{e}_{n-m+\mathrm{2}},\mathrm{\cdots },{e}_{n}\right\}.$

Since ${\mathbf{D}}_{{E}_{f}\\left\{{e}_{i}\right\}}^{\mathrm{T}}\mathbf{N}\left({\mathbf{D}}_{{E}_{f}\\left\{\left\{{e}_{i}\right\}\right\}}^{\mathrm{T}}\right)=\mathrm{0}$, the null space of ${\mathbf{D}}_{{E}_{f}}^{\mathrm{T}}$ and ${\mathbf{D}}_{{E}_{f}\\left\{{e}_{i}\right\}}^{\mathrm{T}}$ has the same number of bases, and thus

$\begin{array}{}\text{(5)}& \text{rank}\phantom{\rule{0.25em}{0ex}}{\mathbf{D}}_{{E}_{f}}^{\mathrm{T}}-\text{rank}\phantom{\rule{0.25em}{0ex}}{\mathbf{D}}_{{E}_{f}\\left\{{e}_{i}\right\}}^{\mathrm{T}}=\mathrm{1}.\end{array}$

Equation (5) implies that

$\begin{array}{}\text{(6)}& \text{rank}\left[{\mathbf{C}}_{{E}_{f}}\phantom{\rule{0.25em}{0ex}}{\mathbf{D}}_{{E}_{f}}\right]-\text{rank}\left[{\mathbf{C}}_{{E}_{f}\\left\{{e}_{i}\right\}}\phantom{\rule{0.25em}{0ex}}{\mathbf{D}}_{{E}_{f}\\left\{{e}_{i}\right\}}\right]=\mathrm{1}.\end{array}$

Then, combining Proposition 2.5 with Expression (6), we get

$\begin{array}{}\text{(7)}& \text{dim}{P}_{{E}_{f}}=k+d-\text{rank}\left[{\mathbf{C}}_{{E}_{f}}\phantom{\rule{0.25em}{0ex}}{\mathbf{D}}_{{E}_{f}}\right]=\text{dim}{P}_{{E}_{f}\\left\{{e}_{i}\right\}}-\mathrm{1}.\end{array}$

${E}_{{E}_{f}\\left\{{e}_{i}\right\}}\subset {E}_{f}$ implies that

$\begin{array}{}\text{(8)}& {P}_{{E}_{f}\\left\{{e}_{i}\right\}}\supseteq {P}_{{E}_{f}}.\end{array}$

Combining Expressions (7) and (8) gives us

${P}_{{E}_{f}\\left\{{e}_{i}\right\}}\supset {P}_{{E}_{f}}.$

Thus, ${P}_{{E}_{f}\\left\{{e}_{i}\right\}}\{P}_{{E}_{f}}\ne \mathit{\varphi }$. Let $\mathbit{x}\in {P}_{{E}_{f}\\left\{{e}_{i}\right\}}\{P}_{{E}_{f}}$, ${\mathbf{C}}_{i}\mathbit{x}+{\mathbf{D}}_{i}\mathbit{y}\ne {b}_{i}$, and thus ${E}_{f}^{\prime }$ is an equality set. We can verify that

$\begin{array}{c}\mathbf{N}{\left({\mathbf{D}}_{{E}_{f}}^{\mathrm{T}}\right)}^{\mathrm{T}}{\mathbf{C}}_{{E}_{f}}=\mathbf{N}{\left({\mathbf{D}}_{{{E}^{\prime }}_{f}}^{\mathrm{T}}\right)}^{\mathrm{T}}{\mathbf{C}}_{{{E}^{\prime }}_{f}}\phantom{\rule{0.125em}{0ex}},\\ \mathbf{N}{\left({\mathbf{D}}_{{E}_{f}}^{\mathrm{T}}\right)}^{\mathrm{T}}{b}_{{E}_{f}}=\mathbf{N}{\left({\mathbf{D}}_{{{E}^{\prime }}_{f}}^{\mathrm{T}}\right)}^{\mathrm{T}}{b}_{{{E}^{\prime }}_{f}}.\end{array}$

From Proposition 2.4, $\text{aff}{P}_{{E}_{f}^{\prime }}=\text{aff}{P}_{{E}_{f}}$, and thus ${E}_{f}^{\prime }\in {\mathrm{\Xi }}_{{F}_{\mathit{\pi }}}$. Therefore, $\text{dim}{P}_{{E}_{f}^{\prime }}\ge \text{dim}{P}_{{E}_{f}\\left\{{e}_{i}\right\}}>\text{dim}{P}_{{E}_{f}}$, and the sufficient part is proven.

Necessary part: the necessary part is proven by contradiction.

Suppose that each row in $\mathbf{N}\left({\mathbf{D}}_{{E}_{f}}^{\mathrm{T}}\right)$ is not a zero row. If $\exists {{E}_{f}}^{*}\in {\mathrm{\Xi }}_{{F}_{\mathit{\pi }}}$, then $\text{dim}{P}_{{{E}_{f}}^{*}}>\text{dim}{P}_{{E}_{f}}$. Let $\text{dim}{P}_{{{E}_{f}}^{*}}-\text{dim}{P}_{{E}_{f}}={n}_{c}$, since the projections of ${{E}_{f}}^{*}$ and Ef have the same affine hull. This can be obtained from Proposition 2.5:

$\text{rank}\phantom{\rule{0.25em}{0ex}}{\mathbf{D}}_{{E}_{f}}-\text{rank}\phantom{\rule{0.25em}{0ex}}{\mathbf{D}}_{{{E}_{f}}^{*}}={n}_{c},$

since ${\mathit{\pi }}_{d}\left(\text{aff}{P}_{{E}_{f}}\right)={\mathit{\pi }}_{d}\left(\text{aff}{P}_{{E}_{f}}^{*}\right)$. Proposition 2.4 implies that $\mathbf{N}\left({\mathbf{D}}_{{E}_{f}}^{\mathrm{T}}\right)$ and $\mathbf{N}\left({\mathbf{D}}_{{{E}_{f}}^{*}}^{\mathrm{T}}\right)$ have the same number of columns, i.e., ${n}^{*}-\text{rank}\phantom{\rule{0.25em}{0ex}}{\mathbf{D}}_{{{E}_{f}}^{*}}=n-\text{rank}\phantom{\rule{0.25em}{0ex}}{\mathbf{D}}_{{E}_{f}}$, $n-{n}^{*}={n}_{c}$.

We can denote $\mathbf{N}\left({\mathbf{D}}_{{E}_{f}^{*}}^{\mathrm{T}}\right)$ as follows.

$\begin{array}{rl}& \mathbf{N}\left({\mathbf{D}}_{{E}_{f}^{*}}^{\mathrm{T}}\right)=\left[\begin{array}{cccc}{N}_{\mathrm{1}}^{*}& {N}_{\mathrm{2}}^{*}& \mathrm{\dots }& {N}_{n-r}^{*}\end{array}\right]\\ & =\left[\begin{array}{cccc}{N}_{\mathrm{11}}^{*}& {N}_{\mathrm{12}}^{*}& \mathrm{\cdots }& {N}_{\mathrm{1}\left(n-r\right)}^{*}\\ {N}_{\mathrm{21}}^{*}& {N}_{\mathrm{22}}& \mathrm{\cdots }& {N}_{\mathrm{2}\left(n-r\right)}^{*}\\ \mathrm{⋮}& \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ {N}_{\left({n}^{*}\right)\mathrm{1}}& {N}_{\left({n}^{*}\right)\mathrm{2}}& \mathrm{\cdots }& {N}_{\left({n}^{*}\right)\left(n-r\right)}\end{array}\right]\end{array}$

Then, we add nc zero elements to each column of $\mathbf{N}\left({\mathbf{D}}_{{E}_{f}^{*}}^{\mathrm{T}}\right)$ as follows.

$\begin{array}{rl}& \mathbf{N}\left({\mathbf{D}}_{{E}_{f}^{*}}^{\prime \mathrm{T}}\right)=\left[\begin{array}{cccc}{N}_{\mathrm{1}}^{\prime *}& {N}_{\mathrm{2}}^{\prime *}& \mathrm{\dots }& {N}_{n-r}^{\prime *}\end{array}\right]\\ & =\left[\begin{array}{cccc}{N}_{\mathrm{11}}^{*}& {N}_{\mathrm{12}}^{*}& \mathrm{\cdots }& {N}_{\mathrm{1}\left(n-r\right)}^{*}\\ {N}_{\mathrm{21}}^{*}& {N}_{\mathrm{22}}& \mathrm{\cdots }& {N}_{\mathrm{2}\left(n-r\right)}^{*}\\ \mathrm{⋮}& \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ {N}_{\left({n}^{*}\right)\mathrm{1}}& {N}_{\left({n}^{*}\right)\mathrm{2}}& \mathrm{\cdots }& {N}_{\left({n}^{*}\right)\left(n-r\right)}\\ \mathrm{0}& \mathrm{0}& \mathrm{\cdots }& \mathrm{0}\\ \mathrm{⋮}& \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ \mathrm{0}& \mathrm{0}& \mathrm{\cdots }& \mathrm{0}\end{array}\right]\end{array}$

Since ${N}_{\mathrm{1}}^{\prime *}$, ${N}_{\mathrm{2}}^{\prime *},\mathrm{\dots },{N}_{n-r}^{\prime *}$ are linearly independent, ${N}_{\mathrm{1}}^{\prime *}$, ${N}_{\mathrm{2}}^{\prime *},\mathrm{\dots },{N}_{n-r}^{\prime *}$ are bases of ${\mathbf{D}}_{{E}_{f}}^{\mathrm{T}}$. There are no zero rows in $\mathbf{N}\left({\mathbf{D}}_{{E}_{f}}^{\mathrm{T}}\right)$, and thus N1, ${N}_{\mathrm{2}},\mathrm{\dots },{N}_{n-r}$ cannot be expressed linearly by ${N}_{\mathrm{1}}^{\prime *}$, ${N}_{\mathrm{2}}^{\prime *},\mathrm{\dots },{N}_{n-r}^{\prime *}$, which is a contradiction. Therefore, the necessary part is proven.

Corollary 3.1 Let ${E}_{f}=\left\{{e}_{\mathrm{1}},{e}_{\mathrm{2}},\mathrm{\dots },{e}_{n},\right\}\in {\mathrm{\Xi }}_{{F}_{\mathit{\pi }}}$. If there exist zero rows in $\mathbf{N}\left({\mathbf{D}}_{{E}_{f}}^{\mathrm{T}}\right)$, then ${E}_{f}^{\prime }$, which consists of elements of Ef corresponding to non-zero rows in $\mathbf{N}\left({\mathbf{D}}_{{E}_{f}}^{\mathrm{T}}\right)$, is the element in ${\mathrm{\Xi }}_{{F}_{\mathit{\pi }}}$ that maximizes the value of $\text{dim}{P}_{{E}_{f}}$.

Corollary 3.1 is a direct consequence of Theorem 3.1, and it implies that it is no longer necessary to discuss the uniqueness of the solutions in the SHOOT and ADJ oracles. However, it only takes an arbitrary solution as the unique solution and then gets the Ef that maximizes the value of $\text{dim}{P}_{{E}_{f}}$ in ${\mathrm{\Xi }}_{{F}_{\mathit{\pi }}}$.

### 3.3.1 Improved SHOOT oracle

Step 1 (selecting a direction γ and finding a point on the surface of P)

By randomly choosing γRd, we solve the following linear programming.

$\begin{array}{}\text{(9)}& \begin{array}{l}\left({r}^{*},{\mathbit{y}}^{*}\right)=\underset{\left(r,\mathbit{y}\right)\in R×{R}^{k}}{\mathrm{arg}\phantom{\rule{0.25em}{0ex}}max}r\\ \text{s.t.}\phantom{\rule{0.25em}{0ex}}\mathbf{C}\mathbit{\gamma }r+\mathbf{D}\mathbit{y}\le \mathbit{b}\end{array}\end{array}$

Step 2 (whether a facet of πd(P) is found)

$\begin{array}{}\text{(10)}& {E}_{{f}_{\mathrm{0}}}=\left\{i|{\mathbf{C}}_{i}\mathit{\gamma }r+{\mathbf{D}}_{i}\mathbit{y}={b}_{i}\right\}\end{array}$

If $\text{Rank}\left(\mathbf{N}{\left({\mathbf{D}}_{{E}_{{f}_{\mathrm{0}}}}^{\mathrm{T}}\right)}^{\mathrm{T}}{\mathbf{C}}_{{E}_{{f}_{\mathrm{0}}}}\right)=\mathrm{1}$, continue. Otherwise, go to Step 1.

Step 3 (ensuring the ${E}_{{f}_{\mathrm{0}}}$ is the element with the maximal $\text{dim}{P}_{{{E}_{f}}_{\mathrm{0}}}$ in ${\mathrm{\Xi }}_{{F}_{\mathit{\pi }}}$)

If there exist zero rows in $\mathbf{N}\left({\mathbf{D}}_{{E}_{{f}_{\mathrm{0}}}}^{\mathrm{T}}\right)$, let ${E}_{{f}_{\mathrm{0}}}←{E}_{{f}_{\mathrm{0}}}^{\prime }$ as Corollary 3.1. Otherwise, continue.

Step 4 (computing the affine hull of the facet and normalizing)

$\begin{array}{l}\left[{\mathbit{a}}_{f}^{\mathrm{T}},{b}_{f}\right]=\mathbf{N}{\left({\mathbf{D}}_{{E}_{{f}_{\mathrm{0}}}}^{\mathrm{T}}\right)}^{\mathrm{T}}\left[{\mathbf{C}}_{{E}_{{f}_{\mathrm{0}}}},{b}_{{E}_{{f}_{\mathrm{0}}}}\right]\\ \left[{\mathbit{a}}_{f}^{\mathrm{T}},{b}_{f}\right]=\frac{\text{sign}{b}_{f}}{{∥{\mathbit{a}}_{f}∥}_{\mathrm{2}}}\left[{\mathbit{a}}_{f}^{\mathrm{T}},{b}_{f}\right]\end{array}$

Step 1 (finding a point on the adjacent facet of πd(P) and its corresponding pre-image)

Solve the following linear programming.

$\begin{array}{l}\left({\mathbit{x}}^{*},{\mathbit{y}}^{*}\right)=\underset{\mathbit{x},\mathbit{y}}{\mathrm{arg}\phantom{\rule{0.25em}{0ex}}max}\phantom{\rule{0.25em}{0ex}}{\mathbit{a}}_{r}^{\mathrm{T}}\mathbit{x}\\ \begin{array}{rl}\text{s.t.}& \phantom{\rule{0.25em}{0ex}}{\mathbf{C}}_{{E}_{r}}\mathbit{x}+{\mathbf{D}}_{{E}_{r}}\mathbit{y}\le {b}_{{E}_{r}},\\ & {\mathbit{a}}_{f}^{\mathrm{T}}\mathbit{x}={b}_{f}\left(\mathrm{1}-\mathit{\epsilon }\right)\end{array}\end{array}$

Step 2 (searching for the Ef of ${\mathrm{\Xi }}_{{F}_{\mathit{\pi }}}$ that maximizes the value of $\text{dim}{P}_{{E}_{f}}$)

If there exist zero rows in $\mathbf{N}\left({\mathbf{D}}_{{E}_{\mathrm{adj}}}^{\mathrm{T}}\right)$, let ${E}_{\mathrm{adj}}←{E}_{\mathrm{adj}}^{\prime }$ as Corollary 3.1. Otherwise, continue.

Step 3 (computing the affine hull of the adjacent facet and normalizing)

$\begin{array}{l}\left[{\mathbit{a}}_{\mathrm{adj}}^{\mathrm{T}},{b}_{\mathrm{adj}}\right]=\mathbf{N}{\left({\mathbf{D}}_{{E}_{\mathrm{adj}}}^{\mathrm{T}}\right)}^{\mathrm{T}}\left[{\mathbf{C}}_{{E}_{\mathrm{adj}}},{b}_{{E}_{\mathrm{adj}}}\right]\\ \left[{\mathbit{a}}_{\mathrm{adj}}^{\mathrm{T}},{b}_{\mathrm{adj}}\right]=\frac{\text{sign}{b}_{\mathrm{adj}}}{{∥{\mathbit{a}}_{\mathrm{adj}}∥}_{\mathrm{2}}}\left[{\mathbit{a}}_{\mathrm{adj}}^{\mathrm{T}},{b}_{\mathrm{adj}}\right]\end{array}$
4 An improved strategy without searching for Er

In this section, we first explain the obstacles to searching for Er in the case of $\text{dim}{P}_{{E}_{f}}>d-\mathrm{1}$, an improved strategy that does not require Er is proposed, and then a new criterion is given to judge whether the two ridges are the same.

## 4.1 The barrier to searching for Er

In the RDG oracle, if $\text{dim}{P}_{{E}_{f}}>d-\mathrm{1}$, this requires a recursive call to the ESP algorithm that reduces the dimension onto which we are projecting at each step.

We note that ar and br can easily be obtained through this strategy, but it is often expensive to search for Er. The only way to do this is to list all of the equality sets ${E}_{c}=\left\{\mathrm{1},\mathrm{2},\mathrm{\dots },q\right\}\{E}_{f}$, $\left\{{E}_{p}^{i}={E}_{f}\cup {E}_{i}\right\}$, and EiEc, and searching for the ${\mathit{\pi }}_{d}\left({E}_{p}^{i}\right)$ is equal to . The computational effort increases in an exponential manner in the number of dimensions.

## 4.2 New criteria without Er

To improve the computational efficiency of the algorithm when the optimizer of LP is not unique, the search for Er has been abandoned. However, in the original ESP algorithm, Er is needed to judge the occurrence times of ridges, and hence we have proposed a new criterion to judge whether the two ridges in πd(P) are the same without Er.

Lemma 4.1 Let P be a polytope. ${F}_{\mathrm{1}}=\left\{\mathbit{x}|{\mathbit{a}}_{{f}_{\mathrm{1}}}^{\mathrm{T}}\mathbit{x}={b}_{{f}_{\mathrm{1}}}\right\}\cap {\mathit{\pi }}_{d}\left(P\right)$ and ${F}_{\mathrm{2}}=\left\{\mathbit{x}|{\mathbit{a}}_{{f}_{\mathrm{2}}}^{\mathrm{T}}\mathbit{x}={b}_{{f}_{\mathrm{2}}}\right\}\cap {\mathit{\pi }}_{d}\left(P\right)$ are the facets of P, ${R}_{\mathrm{1}}=\left\{\mathbit{x}|{\mathbit{a}}_{{r}_{\mathrm{1}}}^{\mathrm{T}}\mathbit{x}={b}_{{r}_{\mathrm{2}}}\right\}\cap {F}_{\mathrm{1}}$ and ${R}_{\mathrm{2}}=\left\{\mathbit{x}|{\mathbit{a}}_{{r}_{\mathrm{2}}}^{\mathrm{T}}\mathbit{x}={b}_{{r}_{\mathrm{2}}}\right\}\cap {F}_{\mathrm{2}}$ are the ridges of P, and R1F1 and R2F2. Then, ${F}_{\mathrm{1}}\cap {F}_{\mathrm{2}}={R}_{\mathrm{1}}={R}_{\mathrm{2}}$ if and only if

$\begin{array}{}\text{(11)}& \text{rank}\left[\begin{array}{c}{\mathbit{a}}_{{f}_{\mathrm{1}}}^{\mathrm{T}}\\ {\mathbit{a}}_{{f}_{\mathrm{2}}}^{\mathrm{T}}\\ {\mathbit{a}}_{{r}_{\mathrm{1}}}^{\mathrm{T}}\\ {\mathbit{a}}_{{r}_{\mathrm{2}}}^{\mathrm{T}}\end{array}\right]=\text{rank}\left[\begin{array}{cc}{\mathbit{a}}_{{f}_{\mathrm{1}}}^{\mathrm{T}}& {b}_{{f}_{\mathrm{1}}}\\ {\mathbit{a}}_{{f}_{\mathrm{2}}}^{\mathrm{T}}& {b}_{{f}_{\mathrm{2}}}\\ {\mathbit{a}}_{{r}_{\mathrm{1}}}^{\mathrm{T}}& {b}_{{r}_{\mathrm{1}}}\\ {\mathbit{a}}_{{r}_{\mathrm{2}}}^{\mathrm{T}}& {b}_{{r}_{\mathrm{2}}}\end{array}\right]=\mathrm{2}.\end{array}$

Proof Only the proof of sufficiency is given here, and the necessity is obvious. Since ${\mathbit{a}}_{{f}_{\mathrm{1}}}\perp {\mathbit{a}}_{{r}_{\mathrm{1}}},{\mathbit{a}}_{{f}_{\mathrm{2}}}\perp {\mathbit{a}}_{{r}_{\mathrm{2}}}$,

$\begin{array}{}\text{(12)}& \text{rank}\left[\begin{array}{cc}{\mathbit{a}}_{{f}_{\mathrm{1}}}^{\mathrm{T}}& {b}_{{f}_{\mathrm{1}}}\\ {\mathbit{a}}_{{r}_{\mathrm{1}}}^{\mathrm{T}}& {b}_{{r}_{\mathrm{1}}}\end{array}\right]=\text{rank}\left[\begin{array}{cc}{\mathbit{a}}_{{f}_{\mathrm{2}}}^{\mathrm{T}}& {b}_{{f}_{\mathrm{2}}}\\ {\mathbit{a}}_{{r}_{\mathrm{2}}}^{\mathrm{T}}& {b}_{{r}_{\mathrm{2}}}\end{array}\right]=\mathrm{2}.\end{array}$

Combining Expression (11) with Expression (12) gives us

$\begin{array}{rl}& \left\{\mathbit{x}|{\mathbit{a}}_{{r}_{\mathrm{1}}}^{\mathrm{T}}\mathbit{x}={b}_{{r}_{\mathrm{1}}}\right\}\cap \left\{\mathbit{x}|{\mathbit{a}}_{{f}_{\mathrm{2}}}^{\mathrm{T}}\mathbit{x}={b}_{{f}_{\mathrm{2}}}\right\}\\ & =\left\{\mathbit{x}|{\mathbit{a}}_{{r}_{\mathrm{2}}}^{\mathrm{T}}\mathbit{x}={b}_{{r}_{\mathrm{2}}}\right\}\cap \left\{\mathbit{x}|{\mathbit{a}}_{{f}_{\mathrm{2}}}^{\mathrm{T}}\mathbit{x}={b}_{{f}_{\mathrm{2}}}\right\}.\end{array}$

Note that

$\begin{array}{l}{R}_{\mathrm{2}}=\left\{\mathbit{x}|{\mathbit{a}}_{{r}_{\mathrm{2}}}^{\mathrm{T}}\mathbit{x}={b}_{{r}_{\mathrm{2}}}\right\}\cap \left\{\mathbit{x}|{\mathbit{a}}_{{f}_{\mathrm{2}}}^{\mathrm{T}}\mathbit{x}={b}_{{f}_{\mathrm{2}}}\right\}\cap {\mathit{\pi }}_{d}\left(P\right),\\ {R}_{\mathrm{1}}=\left\{\mathbit{x}|{\mathbit{a}}_{{r}_{\mathrm{1}}}^{\mathrm{T}}\mathbit{x}={b}_{{r}_{\mathrm{1}}}\right\}\cap \left\{\mathbit{x}|{\mathbit{a}}_{{f}_{\mathrm{1}}}^{\mathrm{T}}\mathbit{x}={b}_{{f}_{\mathrm{1}}}\right\}\cap {\mathit{\pi }}_{d}\left(P\right),\end{array}$

and thus R1=R2. Since the two facets of the polytope can only have one common ridge, ${F}_{\mathrm{1}}\cap {F}_{\mathrm{2}}={R}_{\mathrm{1}}={R}_{\mathrm{2}}$.

## 4.3 Some other patches caused by the abandonment of searching for Er

Since Er is required in other steps in the original ESP algorithm, certain other patches need to be made.

In the ADJ oracle, LP (9) and Eq. (10) can be changed to LP (13) and Eq. (14), respectively.

$\begin{array}{}\text{(13)}& \begin{array}{l}\left({\mathbit{x}}^{*},{\mathbit{y}}^{*}\right)=\underset{\mathbit{x},\mathbit{y}}{\mathrm{arg}\phantom{\rule{0.25em}{0ex}}max}{\mathbit{a}}_{r}^{\mathrm{T}}\mathbit{x}\\ \begin{array}{rl}\text{s.t.}& \phantom{\rule{0.25em}{0ex}}\mathbf{C}\mathbit{x}+\mathbf{D}\mathbit{y}\le \mathbit{b},\\ & {\mathbit{a}}_{f}^{\mathrm{T}}\mathbit{x}={b}_{f}\left(\mathrm{1}-\mathit{\epsilon }\right)\end{array}\end{array}\end{array}$

We can verify that the modified algorithm still gets the correct result, though it also adds certain extra computation cost in the case of $\text{dim}{P}_{{E}_{f}}=d-\mathrm{1}$. We note that $\text{dim}\left({P}_{{E}_{f}}\right)>d-\mathrm{1}$ if and only if the optimizer of LP is non-unique. A switching strategy can be adopted to take the original method until the optimizer of LP is non-unique.

5 Numerical simulations

In this section, the effectiveness of the proposed method is verified based on certain typical orthogonal projection cases and applications to N-step controllable set calculations in the control field.

## 5.1 Numerical simulations for the typical orthogonal projection cases

We select three cases for numerical simulation, and two of them are dual-degenerate, though one is not dual-degenerate.

• i.

The projection of a three-dimensional n prism with n+2 half-spaces onto a two-dimensional plane

• ii.

The projection of an n-dimensional cube with 2n half-spaces onto a two-dimensional plane

• iii.

The projection of an n-dimensional cube onto a three-dimensional space

In cases (i) and (ii), the dual degeneracy occurs, whereas there is no dual degeneracy in case (iii). The simulation operation has been implemented in MATLAB 2019, and the running times of the original ESP and the improved ESP are shown in Fig. 4.

Figure 4Comparative running times for typical cases.

The running time of the original ESP algorithm increases exponentially with dimensions and the number of half-planes of the input polytopes under the dual-degeneracy cases, whereas the improved ESP algorithm significantly reduces the computational burden in the case of dual degeneracy. The running time is also optimized in the case of non-dual degeneracy, since it is no longer necessary to judge the uniqueness of linear programming solutions in the improved ESP algorithm.

## 5.2 Applications to N-step controllable set calculations

We apply the proposed method for the calculation of N-step controllable sets in the control field, and the problem setup is described as given below.

Consider the following linear discrete system.

$\left[\begin{array}{c}{x}_{\mathrm{1}}\left(k+\mathrm{1}\right)\\ {x}_{\mathrm{2}}\left(k+\mathrm{1}\right)\end{array}\right]=\underset{\mathbf{A}}{\underbrace{\left[\begin{array}{cc}\mathrm{1.1}& \mathrm{2}\\ \mathrm{0}& \mathrm{0.95}\end{array}\right]}}\left[\begin{array}{c}{x}_{\mathrm{1}}\left(k\right)\\ {x}_{\mathrm{2}}\left(k\right)\end{array}\right]+\underset{\mathbf{B}}{\underbrace{\left[\begin{array}{c}\mathrm{0}\\ \mathrm{0.0787}\end{array}\right]}}\mathbit{u}$

$\mathbit{x}={\left[{x}_{\mathrm{1}},{x}_{\mathrm{2}}\right]}^{\mathrm{T}}$ and u are the state and input vectors, respectively, subject to the constants as follows.

$\underset{{\mathbf{A}}_{x}}{\underbrace{\left[\begin{array}{cc}\mathrm{0.125}& \mathrm{0}\\ -\mathrm{0.125}& \mathrm{0}\\ \mathrm{0}& \mathrm{0.125}\\ \mathrm{0}& -\mathrm{0.125}\end{array}\right]}}\left[\begin{array}{c}{x}_{\mathrm{1}}\\ {x}_{\mathrm{2}}\end{array}\right]\le \left[\begin{array}{c}\mathrm{1}\\ \mathrm{1}\\ \mathrm{1}\\ \mathrm{1}\end{array}\right]\underset{{\mathbf{A}}_{u}}{\underbrace{\left[\begin{array}{c}-\mathrm{1}\\ \mathrm{1}\end{array}\right]}}\mathbit{u}\le \left[\begin{array}{c}\mathrm{1}\\ \mathrm{1}\end{array}\right]$

Then, we define the target set as follows.

$\underset{{\mathbf{A}}_{f}}{\underbrace{\left[\begin{array}{cc}\mathrm{0.25}& \mathrm{0}\\ -\mathrm{0.25}& \mathrm{0}\\ \mathrm{0}& \mathrm{0.25}\\ \mathrm{0}& -\mathrm{0.25}\end{array}\right]}}\left[\begin{array}{c}{x}_{\mathrm{1}}\\ {x}_{\mathrm{2}}\end{array}\right]\le \left[\begin{array}{c}\mathrm{1}\\ \mathrm{1}\\ \mathrm{1}\\ \mathrm{1}\end{array}\right]$

The N-step controllable set C(N) is defined as the collection of all initial states x0 that can evolve into the target set in N steps while satisfying the state and input constraints, and it can be transformed into the orthogonal projection problem as given below.

$C\left(N\right)={\mathit{\pi }}_{d}\left(P=\left\{\left(\mathbit{x},\mathbit{u}\right)\in {R}^{d}×{R}^{k}|\mathbf{C}\mathbit{x}+\mathbf{D}\mathbit{u}\le \mathbit{b}\right\}\right),$

where d=2, k=1,

$\begin{array}{rl}& \mathbf{D}=\left[\begin{array}{cccc}{\mathbf{A}}_{u}& \mathrm{0}& \mathrm{\cdots }& \mathrm{0}\\ \mathrm{0}& {\mathbf{A}}_{u}& \mathrm{\cdots }& \mathrm{0}\\ \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}\\ \mathrm{0}& \mathrm{0}& \mathrm{\cdots }& {\mathbf{A}}_{u}\\ \mathrm{0}& \mathrm{0}& \mathrm{\cdots }& \mathrm{0}\\ {\mathbf{A}}_{x}\mathbf{B}& \mathrm{0}& \mathrm{\cdots }& \mathrm{0}\\ {\mathbf{A}}_{x}\mathbf{AB}& {\mathbf{A}}_{x}\mathbf{B}& \mathrm{\cdots }& \mathrm{0}\\ \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}\\ {\mathbf{A}}_{f}{\mathbf{A}}^{N-\mathrm{1}}\mathbf{B}& {\mathbf{A}}_{f}{\mathbf{A}}^{N-\mathrm{2}}\mathbf{B}& \mathrm{\cdots }& {\mathbf{A}}_{f}\mathbf{B}\end{array}\right],\\ & \mathbf{C}=\left[\begin{array}{c}\mathrm{0}\\ \mathrm{0}\\ \mathrm{⋮}\\ \mathrm{0}\\ {\mathbf{A}}_{x}\\ {\mathbf{A}}_{x}\mathbf{A}\\ {\mathbf{A}}_{x}{\mathbf{A}}^{\mathrm{2}}\\ \mathrm{⋮}\\ {\mathbf{A}}_{f}{\mathbf{A}}^{N}\end{array}\right],\end{array}$

and

$\mathbit{b}=\left[\begin{array}{c}\mathrm{1}\\ \mathrm{1}\\ \mathrm{⋮}\\ \mathrm{1}\end{array}\right].$

The N-step controllable set can be derived, as shown in Fig. 5, using either the improved ESP algorithm or the original ESP algorithm. When N is small, the region of the N-step controllable set changes significantly with the increasing N, though the trend of change gradually weakens and when N≥12. The N-step controllable set no longer changes with the increase in N.

Figure 5Results of the N-step controllable set.

Moreover, the running times of the original ESP and the improved ESP are shown in Fig. 6. When $N=\mathrm{1},\mathrm{2}$, the two algorithms have similar computational efficiencies. When N>2, the computational efficiency of the improved ESP method is higher than that of the original ESP method. This is due to the fact that, when $N=\mathrm{1},\mathrm{2}$, dual degeneracy does not occur or occurs very rarely, while, when N>3, dual degeneracy occurs frequently. Clearly, the improved ESP algorithm is better than the original ESP method in the application of N-step controllable set computation.

Figure 6Comparative running times for the calculation of the N-step controllable set.

6 Conclusions

In this paper, the improved ESP algorithm becomes simpler, faster, and easier to implement in the case of dual degeneracy. This is achieved by the following methods. First, an alternative criterion has been proposed to directly obtain the equality set that maximizes the value of $\text{dim}{P}_{{E}_{f}}$. It avoids the uniqueness discussion of the optimizer in linear programming, which is not easy to implement. Further, it also reduces the operation times and simplifies the algorithm process. Second, an improved strategy has been presented to solve the difficulty in searching for the equality set Er in the case of $\text{dim}{P}_{{E}_{f}}>d-\mathrm{1}$, which significantly reduces the burden of computation. Finally, the effectiveness of the proposed improved ESP algorithm has been validated through numerical simulations.

Code availability

Software code in this study can be requested from the corresponding author.

Data availability

Research data in this study can be requested from the corresponding author.

Author contributions

BP and WX proposed improved strategies for the ESP algorithm. WX finished the simulation part. YL participated in the correction of the paper.

Competing interests

The contact author has declared that none of the authors has any competing interests.

Disclaimer

Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

Acknowledgements

The authors acknowledge the financial support of the National Natural Science Foundation of China. The authors would like to thank the editors and the anonymous reviewers for their constructive comments and suggestions that have improved the quality of the paper.

Financial support

This research has been supported by the National Natural Science Foundation of China (grant no. 62103440).

Review statement

This paper was edited by Daniel Condurache and reviewed by two anonymous referees.

References

Avis, D.: Computational experience with the reverse search vertex enumeration algorithm, Optim. Method. Softw., 10, 107–124, https://doi.org/10.1080/10556789808805706, 1998. a

Balas, E.: Projection with a minimal system of inequalities, Comput. Optim. Appl., 10, 189–193, https://doi.org/10.1023/a:1018368920203, 1998. a

Balas, E. and Oosten, M.: On the dimension of projected polyhedra, Discrete Appl. Math., 87, 1–9, https://doi.org/10.1016/s0166-218x(98)00096-1, 1998. a, b

Borrelli, F., Bemporad, A., and Morari, M.: Predictive Control for Linear and Hybrid Systems, Cambridge University Press, New York, https://doi.org/10.1017/9781139061759, 2017. a, b

Fukuda, K. and Prodon, A.: Double description method revisited, Springer, Berlin, https://doi.org/10.1007/3-540-61576-8_77, 1996. a

Jones, C. N.: Polyhedral Tools for Control, thesis, Department of Engineering University of Cambridge, Corpus ID: 117208061, 2005. a, b, c, d, e

Jones, C. N., Kerrigan, E. C., and Maciejowski, J. M.: Equality set projection: a new algorithm for the projection of polytopes in halfspace representation, Report CUED/F-INFENG/TR. 463, University of Cambridge, Corpus ID: 14790429, 2004. a, b

Jones, C. N., Kerrigan, E. C., and Maciejowski, J. M.: On polyhedral projection and parametric programming, J. Optimiz. Theory App., 138, 207–220, https://doi.org/10.1007/s10957-008-9384-4, 2008. a

Keerthi, S. S. and Sridharan, K.: Solution of parametrized linear inequalities by Fourier elimination and its applications, J. Optimiz. Theory App., 65, 161–169, https://doi.org/10.1007/bf00941167, 1990.  a

Kwan, M., Sauermann, L., and Zhao, Y. F.: Extension complexity of low-dimensional polytopes, T. Am. Math. Soc., 375, 4209–4250, https://doi.org/10.1090/tran/8614, 2022. a

Lee, Y. I. and Kouvaritakis, B.: Robust receding horizon predictive control for systems with uncertain dynamics and input saturation, Automatica, 36, 1497–1504, https://doi.org/10.1016/s0005-1098(00)00064-9, 2000. a

Liu, H. Y., He, P. Q., Jiang, Y., Zhang, Q. L., and Xia, Q. Q.: A Simple Check Polytope Projection Penalized Algorithm for ADMM Decoding of LDPC Codes, IEEE Access, 11, 2524–2530, https://doi.org/10.1109/access.2023.3234182, 2023. a

Mangasarian, O. L.: Uniqueness of solution in linear programming, Linear Algebra Appl., 25, 151–162, https://doi.org/10.1016/0024-3795(79)90014-4, 1979. a

Rakovic, S. V., Kerrigan, E. C., Mayne, D. Q., and Lygeros, J.: Reachability analysis of discrete-time systems with disturbances, IEEE T. Automat. Contr., 51, 546–561, https://doi.org/10.1109/tac.2006.872835, 2006. a

Sadeghzadeh, A.: Gain-scheduled continuous-time control using polytope-bounded inexact scheduling parameters, Int. J. Robust Nonlin., 28, 5557–5574, https://doi.org/10.1002/rnc.4333, 2018. a

Saraswat, V. and Hentenryck, V. P.: Principles and Practice of Constraint Programming, MIT Press, Boston, ISBN: 978-0-262-19361-0, 1995. a

Teissandier, D. and Delos, V.: Algorithm to calculate the Minkowski sums of 3-polytopes based on normal fans, Comput. Aided Design, 43, 1567–1576, https://doi.org/10.1016/j.cad.2011.06.016, 2011. a

Vincent, T. L. and Wu, Z. Y.: Estimating projections of the controllable set, J. Guid. Control Dynam., 13, 572–575, https://doi.org/10.2514/3.25372, 1990. a

Webster, R. J.: Convexity, Oxford University Press, Oxford, https://doi.org/10.1093/oso/9780198531470.001.0001, 1994. a, b

Xia, Q. Q., Wang, X., Liu, H. J., and Zhang, Q. L.: A Hybrid Check Polytope Projection Algorithm for ADMM Decoding of LDPC Codes, IEEE Commun. Lett., 25, 108–112, https://doi.org/10.1109/lcomm.2020.3023600, 2021. a

Xu, Y. and Guo, Q.: On Monotone Translation-Projection Covariant Minkowski Valuations, Discrete Comput. Geom., 65, 713–729, https://doi.org/10.1007/s00454-019-00069-y, 2021. a

Ziegler, G. M.: Lectures on Polytopes, Springer, New York, https://doi.org/10.1007/978-1-4613-8431-1, 1995. a, b, c