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

. This paper proposes an optimization method for the Equality Set Projection algorithm to compute the orthogonal projection of polytopes. However, its computational burden signiﬁcantly 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: ﬁrst, a new criterion that does not require a discussion of the uniqueness of the solution in linear programming, which simpliﬁes 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.


Introduction
Orthogonal projection is valuable for both theoretical research (Kwan et al., 2022;Xu and Guo, 2021) and engineering applications (Xia et al., 2021;Liu et al., 2023), 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 (Borrelli et al., 2017;Lee and Kouvaritakis, 2000;Sadeghzadeh, 2018).For example, the calculations of the Minkowski sum (Teissandier and Delos, 2011), controllable set (Vincent and Wu, 1990), reachable set (Rakovic et al., 2006), and parametric linear program (Jones et al., 2008) can all be transformed into the projection of polytopes.
Traditional projection algorithms can be grouped into three classes, i.e., Fourier elimination (Keerthi and Sridharan, 1990;Saraswat and Hentenryck, 1995), block elimination (Balas, 1998;Fukuda and Prodon, 1996), and vertex enumeration (Avis, 1998).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 halfspaces.Thus, it is not an efficient algorithm for the highdimensional cases.
The Equality Set Projection (ESP) algorithm was proposed in Jones et al. (2004), and it applies the inspiration of the giftwrapping 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 Published by Copernicus Publications.
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.

Definitions and notation
1.A polytope is a bounded polyhedron defined by the intersection of closed half-spaces.Let P equal (x, y) ∈ R d × R k |Cx + Dy ≤ b ; the orthogonal projection of P onto x space is defined as 3.The affine hull of every face of the polytope is

If
A is a matrix and E ⊆ {1, . .., q}, then A E is a matrix formed by rows of A whose indices are in E. If P is a polytope, then the notation P E refers to the set 5. E is an equality set of P if and only if E = G(E), where 6. N(A) is defined as a matrix whose columns form an orthonormal basis for the null space of A.

Auxiliary results
The following results can be found in Jones (2005), Ziegler (1995), Webster (1994), and Balas and Oosten (1998), respectively.Their proofs are omitted for the sake of brevity.
Theorem 2.1 (Jones, 2005) Given a polytope if E is an equality set of P , then P E is a face of P .Furthermore, if F is a face of P , then there exists a unique equality set such that F = P E .Proposition 2.1 (Ziegler, 1995) If P ⊂ R d × R k is a polytope, then the projection of P onto R d is a polytope.Proposition 2.2 (Webster, 1994) If F 1 , . .., F t are facets of a polytope Proposition 2.5 (Balas and Oosten, 1998) If E is an equality set of polytope

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 Jones (2005).
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 π d The ESP algorithm consists of three key oracles as given below.Shooting Oracle: Ridge Oracle: 15,[183][184][185][186][187][188][189][190][191][192][193]2024 https://doi.org/10.5194/ms-15-183-2024 Adjacency Oracle: The SHOOT oracle randomly obtains an equality set E f 0 of as a facet of π d (P ), the RDG oracle returns all ridges contained in the facet E f , a f , b f of π d (P ), and the ADJ oracle takes two equality sets E r and E f of P and returns the unique equality set E adj , where

Outline of the ESP algorithm
Step 1. Get an initial facet E f 0 , a f , b f of P : E f 0 , a f , b f = SHOOT(P ).Then, all the ridges contained in the facet are calculated: Then, all the ridges contained in the adjacent facet can be calculated: LE r = RDG E adj , a adj , b adj .For each element E r j , a r j , b r j in LE r , if there exists an element in L that already contains E r j , we remove the element from the list.Otherwise, we add element E r j , a r j , b r j , E adj , a adj , b adj to list L. Go to Step 2.

SHOOT oracle
Step 1 (selecting a direction γ and finding a point on the surface of P ) Randomly choose γ ∈ R d and solve the following linear programming.
Step 2 (determining whether a facet of π d (P ) is found) Step 3 (computing the affine hull of the facet and normalizing) Step 4 (handling the dual degeneracy in Step 1) If the linear programming (1) in Step 1 is not dualdegenerate, i.e., (1) has a unique optimizer, then stop.Otherwise, E f 0 needs to be recalculated as follows.
For each row C i x + D i y ≤ b i of P , we compute the following linear programming problem: where

RDG oracle
Step 1 (computing the dimension of and, solving the following linear programming, min (τ,x) where Step 4 (normalizing and orthogonalizing) .

ADJ oracle
Step 1 (finding a point on the adjacent facet π d P E adj and its corresponding pre-image) Solve the following linear programming: where ε is a sufficiently small positive constant.
Step 2 (the case without dual degeneracy) compute the following linear programming problem: where Step 4 (computing the affine hull of the facet and normalizing)

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.

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 (x * , y * ) is on π d P E adj (Jones, 2005).If the optimizer (x * , y * ) of linear programming is unique, E adj should be taken as the equality set corresponding to the face P E adj where (x * , y * ) is located.If the optimizer is not unique, E adj should be taken as the equality set corresponding to the face P E 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 , the orthogonal projection π d (P ) = {Gx ≤ g} onto x space R d is shown in Fig. 1, where Figure 2 shows that the ADJ oracle takes the facet a T f x = b f ∩ π (P ) of P (the blue line segment in this example) and ridge a T r x = b r ∩ π d P E f (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 a T f x = b f (1 − ε) and the yellow line segment.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 adj , and its corresponding equality set is E adj = {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 (x * , y * ) on the intersection of the blue facet and yellow facet is taken as the unique optimal solution, the P E adj that (x * , y * ) belongs to is the intersection line of the blue and yellow facets of P , E adj = {1, 4}.According to Fig. 3, P E adj ⊂ P E adj and π d P E adj ⊂ π d P E adj .
Owing to affπ d P E adj = affπ d P E adj , 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.

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 pro-gramming.Although several achievements have been made on the unique conditions of LP solutions (Mangasarian, 1979), 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 (Borrelli et al., 2017), it is conservative to determine the uniqueness of the solution by testing whether the dual degeneracy occurs in Jones et al. (2004).
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 adj does not vary.In this case, an extra computation is added to handle the non-uniqueness of the optimal solution of the LP.

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 E f satisfying this condition can be denoted as where The discussion on the uniqueness of optimizers of LP is to ensure that the obtained E f is the element in F π with the maximal dimP E f .The primary result can be stated as given below.Theorem 3.1 Consider the F π defined in Expression (4), and E f is the subset of F π that maximizes the value of dimP E f if and only if there are no zero rows in N D T E f .Proof (sufficient part): the sufficiency of Theorem 3.1 has been proven by the following inverse negative proposition.
If there exist zero rows in N D , there is a one-to-one correspondence between N row i and D T ei ; i.e., there is a one-to-one correspondence between the rows of N D T E f and the elements of E f .Let E f = {e 1 , e 2 , . .., e n−m }.
. .N n−r forms a basis for the null space of D T E f , implying that N 1 N 2 . . .N n−r forms a basis for the null space of D T E f .Next, we will prove that E f is an equality set, i.e., Since E f is an equation set, we only need to prove that has the same number of bases, and thus rank (5) Equation ( 5) implies that rank Then, combining Proposition 2.5 with Expression ( 6), we get 7) Combining Expressions ( 7) and ( 8) gives us and thus E f is an equality set.We can verify that From Proposition 2.4, affP E f = affP E f , and thus and the sufficient part is proven.
Necessary part: the necessary part is proven by contradiction.
Suppose that each row in * and E f have the same affine hull.This can be obtained from Proposition 2.5: There are no zero rows in N D T E f , and thus N 1 , N 2 , . .., N n−r cannot be expressed linearly by N * 1 , N * 2 , . .., N * n−r , which is a contradiction.Therefore, the necessary part is proven.
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 E f that maximizes the value of dimP E f in F π .

Improved SHOOT oracle
Step 1 (selecting a direction γ and finding a point on the surface of P ) By randomly choosing γ ∈ R d , we solve the following linear programming.
Step 2 (whether a facet of π d (P ) is found) Step 1.
Step 3 (ensuring the E f 0 is the element with the maximal dimP Step 4 (computing the affine hull of the facet and normalizing) Step 1 (finding a point on the adjacent facet of π d (P ) and its corresponding pre-image) Solve the following linear programming.
Step 2 (searching for the E f of F π that maximizes the value of dimP E f ) If there exist zero rows in N D T E adj , let E adj ← E adj as Corollary 3.1.Otherwise, continue.
Step 3 (computing the affine hull of the adjacent facet and normalizing) In this section, we first explain the obstacles to searching for E r in the case of dimP E f > d − 1, an improved strategy that does not require E r is proposed, and then a new criterion is given to judge whether the two ridges are the same.

The barrier to searching for E r
In the RDG oracle, if dimP E f > d − 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 a r and b r can easily be obtained through this strategy, but it is often expensive to search for E r .The only way to do this is to list all of the equality sets To improve the computational efficiency of the algorithm when the optimizer of LP is not unique, the search for E r has been abandoned.However, in the original ESP algorithm, E r 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 E r .Lemma 4.1 Let P be a polytope.
Proof Only the proof of sufficiency is given here, and the necessity is obvious.Since Combining Expression (11) with Expression (12) gives us Note that and thus R 1 = R 2 .Since the two facets of the polytope can only have one common ridge, 4.3 Some other patches caused by the abandonment of searching for E r Since E r 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.We can verify that the modified algorithm still gets the correct result, though it also adds certain extra computation cost in the case of dimP E f = d − 1.We note that dim P E f > d − 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.

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.

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 halfspaces onto a two-dimensional plane iii.The projection of an n-dimensional cube onto a threedimensional 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.
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.

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.
x 1 (k + 1) T and u are the state and input vectors, respectively, subject to the constants as follows.
Then, we define the target set as follows.
The N -step controllable set C(N) is defined as the collection of all initial states x 0 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. .
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 Nstep 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 .Moreover, the running times of the original ESP and the improved ESP are shown in Fig. 6.When N = 1, 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 = 1, 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.

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 https://doi.org/10.5194/ms-15-183-2024 Mech.Sci., 15, 183-193, 2024  equality set that maximizes the value of dimP 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 E r in the case of dimP E f > d − 1, which significantly reduces the burden of computation.Finally, the effectiveness of the proposed improved ESP algorithm has been validated through numerical simulations.
= b and c T x + d T y ≤ b for all (x, y) ∈ P .The face of the polytope P ⊂ R d × R k refers to the face of dimension n − 1, whereas the ridge refers to the face of dimension n − 2.
stop and report G and g.Step 3. Search the adjacent facet such that E adj , a adj , b adj = ADJ ((E fn , a fn , b fn ) , (E rn , a rn , b rn )) and let π d−1 ( P ) can be computed using a recursive call to the ESP algorithm.Each facet of π d−1 ( P ) can be converted into the ridge of π d (P ) as x| Ṽ ãT f x = bf ∩ π d−1 ( P ).https://doi.org/10.5194/ms-15-183-2024Mech.Sci., 15, 183-193, 2024 Then, we get a r = ãf , b r = bf , and (Q(i), a r , b r ) to LE r .

Figure 1 .
Figure 1.Polytope and its projection onto R d .

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

4
An improved strategy without searching for E r and E i ⊂ E c , and searching for the π d E i p is equal to x ãT f x = bf ∩ π d−1 ( P ).The computational effort increases in an exponential manner in the number of dimensions.https://doi.org/10.5194/ms-15-183-2024Mech.Sci., 15, 183-193, 2024 4.2 New criteria without E r

Figure 4 .
Figure 4. Comparative running times for typical cases.

Figure 5 .
Figure 5. Results of the N -step controllable set.

Figure 6 .
Figure 6.Comparative running times for the calculation of the Nstep controllable set.
T E f , then E f is not the element in F π that maximizes the value of dimP E f .Let E f = {e 1 , e 2 , . .., e n , }, D T E = D T e 1 D T e 2 . . .D T e n ; without loss of generality, we https://doi.org/10.5194/ms-15-183-2024Mech.Sci., 15, 183-193, 2024 B. Pei et al.: Improved strategies of the ESP algorithm can assume that the last m rows of N D T E f are zero rows.