Exploiting Algebraic Structure in Global Optimization and the Belgian Chocolate Problem
Zachary Charles Nigel Boston


[
    [
    December 30, 2023
=====================




The Belgian chocolate problem involves maximizing a parameter δ over a non-convex region of polynomials. In this paper we detail a global optimization method for this problem that outperforms previous such methods by exploiting underlying algebraic structure. Previous work has focused on iterative methods that, due to the complicated non-convex feasible region, may require many iterations or result in non-optimal δ. By contrast, our method locates the largest known value of δ in a non-iterative manner. We do this by using the algebraic structure to go directly to large limiting values, reducing the problem to a simpler combinatorial optimization problem. While these limiting values are not necessarily feasible, we give an explicit algorithm for arbitrarily approximating them by feasible δ. 
Using this approach, we find the largest known value of δ to date, δ = 0.9808348. 
We also demonstrate that in low degree settings, our method recovers previously known upper bounds on δ and that prior methods converge towards the δ we find.




§ INTRODUCTION


	Global optimization problems of practical interest can often be cast as optimization programs over non-convex feasible regions. Unfortunately, iterative optimization over such regions may require large numbers of iterations and result in non-global maxima. Finding all or even many critical points of such programs is generally an arduous, computationally expensive task. 
	In this paper we show that by exploiting the underlying algebraic structure, we can directly find the largest known values of the Belgian chocolate problem, a famous open problem bridging optimization and control theory. Moreover, this algebraic method does not require any iterative approach. Instead of relying on eventual convergence, our method algebraically identifies points that provide the largest value of the Belgian chocolate problem so far.

	While this approach may seem foreign to the reader, we will show that our algebraic optimization method outperforms prior global optimization methods for solving the Belgian chocolate problem. We will contrast our method with the optimization method of Chang and Sahinidis <cit.> in particular. Their method used iterative branch-and-reduce techniques <cit.> to find what was the largest known value of δ until our new approach. Due to the complicated feasible region, their method may take huge numbers of iterations or converge to suboptimal points. Our method eliminates the need for these expensive iterative computations by locating and jumping directly to the larger values of δ. This approach has two primary benefits over <cit.>. First, it allows us to more efficiently find δ as we can bypass the expensive iterative computations. This also allows us to extend our approach to cases that were not computationally tractable for <cit.>. Second, our approach allows us to produce larger values of δ by finding a finite set of structured limit points. In low-degree cases, this set provably contains the supremum of the problem, while in higher degree cases, the set contains larger values of δ than found in <cit.>.

	The Belgian chocolate problem is a famous open problem in control theory proposed by Blondel in 1994. In the language of control theory, Blondel wanted to determine the largest value of a process parameter for which stabilization of an unstable plant could be achieved by a stable minimum-phase controller <cit.>. Blondel designed the plant to be a low-degree system that was resistant to known stabilization methods, in the hope that a solution would lead to development of new stabilization techniques. Specifically, Blondel wanted to determine the largest value of δ > 0 for which the transfer function P(s) = (s^2-1)/(s^2-2δ s+1) can be stabilized by a proper, bistable controller.

	For readers unfamiliar with control theory, this problem can be stated in simple algebraic terms. To do so, we will require the notion of a stable polynomial. A polynomial is stable if all its roots have negative real part. The Belgian chocolate problem is then as follows.

	

	Belgian chocolate problem: Determine for which δ > 0 there exist real, stable polynomials x(s), y(s), z(s) with (x) ≥ (y) satisfying 
	
    z(s) = (s^2-2δ s+1)x(s)+(s^2-1)y(s).

	We call such δ admissible. In general, stability of x,y,z becomes harder to achieve the larger δ is. Therefore, we are primarily interested in the supremum of all admissible δ. If we fix a maximum degree n for x and y, then this gives us the following global optimization problem for each n.

	

	Belgian chocolate problem (optimization version):
	
    δ, x(s), y(s)maximize      δ
       subject to       x, y, z are stable,
              z(s) = (s^2-2δ s+1)x(s)+(s^2-1)y(s),
             (y) ≤(x) ≤ n.

	Note that we can view a degree n polynomial with real coefficients as a (n+1)-dimensional real vector of its coefficients. Under this viewpoint, the space of polynomials x, y, z that are stable and satisfy (<ref>) is an extremely complicated non-convex space. As a result, it is difficult to employ global optimization methods directly to this problem. The formulation above does suggest an undercurrent of algebra in this problem. This will be exploited to transform the problem into a combinatorial optimization problem by finding points that are essentially local optima.

	Previous work has employed various optimization methods to find even larger admissible δ. 
	Patel et al. <cit.> were the first to show that δ = 0.9 is admissible by x, y of degree at most 11, answering a long-standing question of Blondel.
	They further showed that δ = 0.93720712277 is admissible. In 2005, Burke et al. <cit.> showed that δ = 0.9 is admissible with x,y of degree at most 3. They also improved the record to δ = 0.94375 using gradient sampling techniques. In 2007, Chang and Sahinidis used branch-and-reduce techniques to find admissible δ as large as 0.973974 <cit.>. In 2012, Boston used algebraic techniques to give examples of admissible δ up to 0.97646152 <cit.>. Boston found polynomials that are almost stable and satisfy (<ref>). Boston then used ad hoc methods to perturb these to find stable x,y,z satisfying (<ref>). While effective, no systematic method for perturbing these polynomials to find stable ones was given.

	In this paper, we extend the approach used by Boston in 2012 <cit.> to achieve the largest known value of δ so far. We will refer to this method as the method of algebraic specification. We show that these almost stable polynomials serve as limiting values of the optimization program. Empirically, these almost stable polynomials achieve the supremum over all feasible δ. Furthermore, we give a theoretically rigorous method for perturbing the almost stable polynomials produced by algebraic specification to obtain stable polynomials. Our approach shows that all δ≤ 0.9808348 are admissible. This gives the largest known admissible value of δ to date. We further show that previous global optimization methods are tending towards the limiting values of δ found via our optimization method.

	We do not assume any familiarity on the reader's part with the algebra and control theory and will introduce all relevant notions. While we focus on the Belgian chocolate problem throughout the paper, we emphasize that the general theme of this paper concerns the underlying optimization program. We aim to illustrate that by considering the algebraic structure contained within an optimization problem, we can develop better global optimization methods.



§ MOTIVATION FOR OUR APPROACH


	In order to explain our approach, we will discuss previous approaches to the Belgian chocolate problem in more detail. Such approaches typically perform iterative non-convex optimization in the space of stable controllers in order to maximize δ. In <cit.>, Chang and Sahinidis formulated, for each n, a non-convex optimization program that sought to maximize δ subject to the polynomials x, y, (s^2-2δ s +1)x+(s^2-1)y being stable and such that n ≥(x) ≥(y). For notational convenience, we will always define z = (s^2-2δ s + 1)x+(s^2-1)y. Chang and Sahinidis used branch-and-reduce techniques to attack this problem for n up to 10.

	Examining the roots of the x,y,z they found for (x) = 6,8,10, a pattern emerges. Almost all the roots of these polynomials are close to the imaginary axis and are close to a few other roots. In fact, most of these roots have real part in the interval (-0.01,0). In other words, the x,y,z are approximated by polynomials with many repeated roots on the imaginary axis. It is also worth noting that the only roots of x that were omitted are very close to -δ±√(δ^2-1). This suggests that x should have a factor close to (s^2+2δ s+1).

	This suggests the following approach. Instead of using non-convex optimization to iteratively push x,y,z towards polynomials possessing repeated roots on the imaginary axis, we will algebraically construct polynomials with this property. This will allow us to immediately find large limit points of the optimization problem in (<ref>). While the x,y,z we construct are not stable, they are close to being stable. We will show later that we can perturb x,y,z and thereby push their roots just to the left of the imaginary axis, causing them to be stable. This occurs at the expense of decreasing δ by an arbitrarily small amount.

	Our method only requires examining finitely many such limit points. Moreover, for reasonable degrees of x and y, these limit points can be found relatively efficiently. By simply checking each of these limit points, we reduce to a combinatorial optimization problem. This combinatorial optimization problem provably achieves the supremal values of δ for (x) ≤ 4. For higher degree x, our method finds larger values of δ than any previous optimization method thus far. In the sections below we will further explain and motivate our approach, and show how this leads to the largest admissible δ found up to this point.



§ MAIN RESULTS


	

 §.§ Preliminaries


	Given t ∈, we let (t) denote its real part. We will let [s] denote the set of polynomials in s with real coefficients. For p(s) ∈[s], we call p(s) stable if every root t of p satisfies (t) < 0. We let H denote the set of all stable polynomials in [s]. We call p(s) quasi-stable if every root t of p satisfies (t) ≤ 0. We let H denote the set of quasi-stable polynomials of [s]. We let H^m, H^m denote the sets of stable and quasi-stable polynomials respectively of degree at most m.

	We call δ admissible if there exist x,y ∈ H such that (x) ≥(y) and
	
    (s^2-2δ s+1)x(s) + (s^2-1)y(s) ∈ H.


	We call δ quasi-admissible if there exist x,y ∈H such that (x) ≥(y) and
	
    (s^2-2δ s+1)x(s) + (s^2-1)y(s) ∈H.


	Note that since quasi-stability is weaker than stability, quasi-admissibility is weaker than admissibility. Our main theorem (Theorem <ref> below) will show that if δ is quasi-admissible, then all smaller δ are admissible. Note that this implies that the Belgian chocolate problem is equivalent to finding the supremum of all admissible δ. We will then find quasi-admissible δ in order to establish which δ are admissible. This is the core of our approach. These quasi-admissible δ are easily identified and are limit points of admissible δ.

	In practice, one verifies stability by using the Routh-Hurwitz criteria. Suppose we have a polynomial p(s) = a_0s^n + a_1s^n-1 + … + a_n-1s + a_n ∈[s] such that a_0 > 0. Then we define the n× n Hurwitz matrix A(p) as

	
    A(p) = [   a_1   a_3   a_5     …     …     0     0;   a_0   a_2   a_6     …     …     0     0;     0   a_1   a_3     …     …     0     0;     0   a_0   a_2     …     …     0     0;     ⋮     ⋮     ⋮     ⋱     ⋱     ⋮     ⋮;     0     0     0     …     … a_n-2   a_n ].


	Adolf Hurwitz showed that a real polynomial p with positive leading coefficient is stable if and only if all leading principal minors of A(p) are positive.
	While it may seem natural to conjecture that p is quasi-stable if and only if all leading principal minors are nonnegative, this only works in one direction.

	Suppose p is a real polynomial with positive leading coefficient. If p is quasi-stable then all the leading principal minors of A(p) are nonnegative.
	If p(s) is quasi-stable, then for all ϵ > 0, p(s+ϵ) is stable. Therefore, for all ϵ > 0, the leading minors of A(p(s+ϵ)) are all positive. Note that
	
    lim_ϵ→ 0 A(p(s+ϵ)) = A(p).


	Since the minors of a matrix are expressible as polynomial functions of the entries of the matrix, the leading principal minors of A are limits of positive real numbers. They are therefore nonnegative.

	To see that the converse doesn't hold, consider p(s) = s^4 + 198s^2 + 101^2. Its Hurwitz matrix has nonnegative leading principal minors, but p is not quasi-stable. This example, as well as a more complete characterization of quasi-stability given below, can be found in <cit.>. In particular, it is shown in <cit.> that a real polynomial p with positive leading coefficient is quasi-stable if and only if for all ϵ > 0, A(p(s+ϵ)) has positive leading principal minors.

	

 §.§ Quasi-admissible and admissible δ


	We first present the following theorem concerning which δ are admissible. We will defer the proof until later as it is a simple corollary to a stronger theorem about approximating polynomials in H by polynomials in H.

	If δ is admissible then all δ̂ < δ are also admissible.

	For δ = 1, note that the Belgian chocolate problem reduces to whether there are x,y ∈ H with (x) ≥(y) such that (s-1)^2x + (s^2-1)y ∈ H. This cannot occur for non-zero x,y since (s-1)^2x + (s^2-1)y has a root at s = 1. Theorem <ref> then implies that any δ≥ 1 is not admissible. In 2012, Bergweiler and Eremenko showed that any admissible δ must satisfy δ < 0.999579 <cit.>.

	On the other hand, if we fix x,y then there is no single largest admissible δ associated to x,y. Standard results from control theory show that if δ is admissible by x, y then for ϵ small enough, δ+ϵ is admissible by the same polynomials.

	Therefore, supremum δ^* over all admissible δ will not be associated to stable x,y. From an optimization point of view, the associated optimization program in (<ref>) has an open feasible region. In particular, the set of admissible δ for (<ref>) is of the form (0,δ_n^*) for some δ_n^* that is not admissible by x,y of degree at most n. However, as we will later demonstrate, quasi-admissible δ lie on the boundary of this feasible region. Moreover, quasi-admissible δ naturally serve as analogues of local maxima. We will therefore find quasi-admissible δ and use these to find admissible δ. In Section <ref> we will prove the following theorem relating admissible and quasi-admissible δ. The following is the main theorem of our work and demonstrates the utility of searching for quasi-admissible δ.

	If δ is quasi-admissible, then all δ̂ < δ are admissible. Moreover, if δ is quasi-admissible by quasi-stable x,y of degree at most n, then any δ̂ < δ is admissible by stable x̂, ŷ of degree at most n. 

	This theorem shows that to find admissible δ, we need only to find quasi-admissible δ. In fact our theorem will show that if δ is quasi-admissible via x,y of degree at most n, then all δ̂ < δ are admissible via x,y of degree at most n as well. In short, quasi-admissible δ serve as upper limit points of admissible δ. Also note that since admissible implies quasi-admissible, Theorem <ref> implies Theorem <ref>.

	The proof of Theorem <ref> will be deferred until Section <ref>. In fact, we will do more than just prove the theorem. We will given an explicit algoritm for approximating quasi-stable δ̂ by stable δ within any desired tolerance. We will also be able to use the techniques in Section <ref> to prove the following theorem showing that admissible δ are always smaller than some quasi-admissible δ.

	If δ is admissible by x,y of degree at most n then there is some δ̂ > δ that is quasi-admissible by x̂, ŷ of degree at most n. Moreover, this δ̂ is not admissible by these polynomials.

	In other words, for any admissible δ, there is a larger δ̂ that is quasi-admissible but not necessarily admissible. Therefore, we can restrict to looking at polynomials x,y,z with at least one root on the imaginary axis.



§ LOW DEGREE EXAMPLES


	In this section we demonstrate that in low-degree settings, the supremum of all admissible δ in (<ref>) is actually a quasi-admissible δ. By looking at quasi-stable polynomials that are not stable, we can greatly reduce our search space and directly find the supremum of the optimization program in (<ref>). 
	For small degrees of x, y, we will algebraically design quasi-stable polynomials that achieve previously known bounds on the Belgian chocolate problem in these degrees.

	Burke et al. <cit.> showed that for x∈ H^3, y ∈ H^0, any admissible δ must satisfy δ < √(2+√(2))/2 and for x ∈ H^4, y ∈ H^0, δ must satisfy δ < √(10+2√(5))/4. He et al. <cit.> later found x ∈ H^4, y ∈ H^0 admitting δ close to this bound.

	In fact, these upper bounds on admissible δ are actually quasi-admissible δ that can be obtained in a straightforward manner. For example, suppose we restrict to x of degree 3, y of degree 0. Then for some A, B, C, k ∈, we have
	
    x(s) = s^3+ As^2 + Bs + C
    
    	y(s) = k


	Instead of trying to find admissible δ using this x and y, we will try to find quasi-admissible δ. That is, we want δ such that
	
    z(s) = (s^2-2δ s + 1)x(s) + (s^2-1)y(s) ∈H.


	In other words, this z(s) can be quasi-stable instead of just stable. Note that z(s) must be of degree 5. We will specify a form for z(s) that ensures it is quasi-stable. Consider the case z(s) = s^5. This is clearly quasi-stable as its only roots are at s = 0. To ensure that z(s) = s^5 and equation (<ref>) holds, we require
	
    (s^2-2δ s + 1)(s^3+ As^2 + Bs + C)+(s^2-1)k = s^5


	Equating coefficients gives us the following 5 equations in 5 unknowns.
	
    A-2δ=0
    
    	-2Aδ + B + 1=0
    
    	A - 2Bδ + C + k=0
    
    	B-2Cδ=0
    
    	C-k=0


	In fact, ensuring that we have as many equations as unknowns was part of the motivation for letting z(s) = s^5. Solving for A,B,C,k,δ, we find
	
    8δ^4-8δ^2+1=0
    
    	A = 2δ
    
    	B = 4δ^2-1
    
    	C = 4δ^3-2δ
    
    	k = 4δ^3-2δ


	Taking the largest real root of 8δ^4-8δ^2+1 gives δ = √(2+√(2))/2. Taking A,B,C,k as above yields polynomials x, y, z with real coefficients. One can verify that x is stable (via the Routh-Hurwitz test, for example), while y is degree 0 and therefore stable. Note that since z(s) = s^5, z is only quasi-stable. Therefore, there is x ∈ H^3, y ∈ H^0 for which √(2+√(2))/2 is quasi-admissible. This immediately gives the limiting value for x ∈ H^3, y ∈ H^0 discovered by Burke et al <cit.>. Combining this with Theorem <ref>, we have shown the following theorem.
	For (x) ≤ 3, δ = √(2+√(2))/2 is quasi-admissible and all δ < √(2+√(2))/2 are admissible.
	Next, suppose that x has degree 4 and y has degree 0. For A, k, δ∈, define
	
    x(s) = (s^2+2δ s + 1)(s^2+A)
    
    	y(s) = k


	Note that as long as A ≥ 0, x will be quasi-stable and y will be stable for any k. As above, we want quasi-admissible δ. We let z(s) = s^6, so that z(s) is quasi-stable. Finding A, δ, k amounts to solving
	
    (s^2-2δ s + 1)x(s) + (s^2-1)y(s) = z(s)
    ⇔ (s^2-2δ s + 1)(s^2+2δ s+1)(s^2+A)+(s^2-1)k = s^6
    ⇔ s^6 + (A - 4δ^2 + 2)s^4 + (-4Aδ^2 + 2A + k + 1)s^2 + (A-k) = s^6


	Note that the (s^2+2δ s + 1) term in x is used to ensure that the left-hand side will have zero coefficients in its odd degree terms. Since (s^2+2δ s + 1) is stable, it does not affect stability of x. Equating coefficients and manipulating, we get the following equations.
	
    16δ^4 -20δ^2+5=0
    
    	A -4δ^2+2=0
    
    	k -A=0


	Taking the largest real root of 16δ^4 -20δ^2+5 gives δ = √(10+2√(5))/4. For this δ one can easily see that A = 4δ^2 - 2 ≥ 0, so x is quasi-stable, as are y and z by design. Once again, we were able to easily achieve the limiting value discovered by Burke et al. <cit.> discussed in Section <ref> by searching for quasi-admissible δ. Combining this with Theorem <ref>, we obtain the following theorem.

	For (x) ≤ 4, δ  = √(10+2√(5))/4 is quasi-admissible and all δ < √(10+2√(5))/4 are admissible.

	The examples above demonstrate how, by considering quasi-stable x,y and z, we can find quasi-admissible δ that are limiting values of admissible δ. Moreover, the quasi-stable δ above were found by solving relatively simple algebraic equations instead of having to perform optimization over the space of stable x and y.



§ ALGEBRAIC SPECIFICATION


	The observations in Section <ref> and Section <ref> and the examples in Section <ref> suggest the following approach which we refer to as algebraic specification. This method will be used to find the largest known values of δ found for any given degree. We wish to construct quasi-stable x(s), y(s), z(s) with repeated roots on the imaginary line satisfying (<ref>). For example, we may wish to find polynomials of the following form:
	
    x(s) = (s^2+2δ s+1)(s^2+A_1)^4(s^2+A_2)^2(s^2+A_3)^2(s^2+A_4)
    
    	y(s) = k(s^2+B_1)^3(s^2+B_2)^2
    
    	z(s) = s^14(s^2+C_1)^2(s^2+C_2)(s^2+C_3)


	We refer to such an arrangement of x,y,z as an algebraic configuration. As long as δ > 0, the parameters {A_i}_i=1^4, {B_i}_i=1^2, and {C_i}_i=1^3 are all nonnegative, and k is real, x(s), y(s), z(s) will be real, quasi-stable polynomials. We then wish to solve
	
    (s^2-2δ s+1)x(s)+(s^2-1)y(s)=z(s)


	Recall that the (s^2+2δ s+1) factor in x(s) is present to ensure that the left-hand side has only even degree terms, as the right-hand side clearly only has even degree terms. Expanding (<ref>) and equating coefficients, we get 11 equations in 11 unknowns. Using PHCPack <cit.> to solve these equations and selecting the solution with the largest δ such that the A_i, B_i, C_i ≥ 0, we get the following solution, rounded to seven decimal places:
	
    δ = 0.9808348
    
    	A_1 = 1.1856917
    
    	A_2 = 6.6228807
    
    	A_3 = 0.3090555
    
    	A_4 = 0.2292503
    
    	B_1 = 0.5430391
    
    	B_2 = 0.2458118
    
    	C_1 = 4.4038385
    
    	C_2 = 0.7163490
    
    	C_3 = 7.4637156
    
    	k = 196.1845537


	The actual solution has δ = 0.980834821202…. This is the largest δ we have found to date using this method. By Theorem <ref>, we conclude the following theorem.

	All δ≤ 0.9808348 are admissible.

	In general, we can form an algebraic configuration for x(s), y(s), z(s) as
	
    x(s) = (s^2+2δ s +1) ∏_i=1^m_1 (s^2+A_i)^j_i.

	
    y(s) = k∏_i=1^m_2(s^2+B_i)^k_i.

	
    z(s) = s^c∏_i=1^m_3(s^2+C_i)^ℓ_i.


	For fixed degrees of x, y, note there are only finitely many such configurations. Instead of performing optimization over the non-convex feasible region of the Belgian chocolate problem, we instead tackle the combinatorial optimization problem of maximizing δ among the possible configurations. 

	Note that c in (<ref>) is whatever exponent is needed to make (z) = (x)+2. We want x,y,z to satisfy (<ref>). Expanding and equating coefficients, we get equations in the undetermined variables above. As long as the number of unknown variables equals the number of equations, we can solve and look for real solutions with δ and all A_i, B_i, C_i nonnegative.

	Not all quasi-stable polynomials can be formed via algebraic specification. In particular, algebraic specification forces all the roots of y,z and all but two of the roots of x to lie on the imaginary axis. However, more general quasi-stable x,y,z could have some roots with negative real part and some with zero real part. This makes the possible search space infinite and, as discussed in Section <ref>, empirically does not result in larger δ. Further evidence for this statement will be given in Section <ref>.

	While the method of algebraic specification has demonstrable effectiveness, it becomes computationally infeasible to solve these general equations for very large n. In particular, the space of possible algebraic configurations of x,y,z grows almost exponentially with the degree of the polynomials. For large n, an exhaustive search over the space of possible configurations becomes infeasible, especially as the equations become more difficult to solve.

	We will describe an algebraic configuration via the shorthand 
	
    [j_1,…,j_m_1],[k_1,…, k_m_2],[ℓ_1,…, ℓ_m_3].

	This represents the configuration described in (<ref>),(<ref>),(<ref>) above. In particular, if the second term of (<ref>) is empty then y = k, while if the third term of (<ref>) is empty then z is a power of s. For example, the following configuration is given by [3,1],[2],[1]:
	
    x(s) = (s^2+2δ s + 1)(s^2+A_1)^3(s^2+A_2)
    
    	y(s) = k(s^2+B_1)^2
    
    	z(s) = s^10(s^2+C_1)


	A table containing the largest quasi-admissible δ we have found and their associated algebraic configuration for given degrees of x is given below. Note that for each entry of the table, given (x) = n and quasi-admissible δ, Theorem <ref> implies that all δ̂ < δ are admissible with x,y of degree at most n.

	



§ APPROXIMATING QUASI-ADMISSIBLE Δ BY ADMISSIBLE Δ


	In this section we will prove Theorem <ref>. Our proof will be algorithmic in nature. We will describe an algorithm that, given δ that is quasi-admissible by quasi-stable polynomials x, y, will produce for any δ̂ < δ stable polynomials x̂, ŷ admitting δ̂. Moreover, given (x) = n, we will ensure that (x̂) ≤ n.

	Suppose that for a given δ there are x,y,z ∈H with (x) ≥(y) satisfying (<ref>). Let n = (x). Define
	
    R(s) := (s^2-1)y(s)z(s).

	
	Note that for any s ∈, R(s) = 0 iff (s^2-1)y(s) = 0, R(s) = 1 iff (s^2-2δ s+1)x(s) = 0, and R(s) is infinite iff z(s) = 0. Since x,y,z are quasi-stable, we know that for (s) > 0, R(s) = 1 iff s = δ± i√(1-δ^2) and R(s) = 0 iff s = 1. All other points where R(s) is 0, 1, or infinite satisfy (s) ≤ 0.
	Precomposing R(s) with the fractional linear transformation f(s) = (1+s)/(1-s), we get the complex function
	
    D(s) := R(1+s1-s).


	Note that this fractional linear transformation maps the unit disk {s | |s| = 1} to the imaginary axis { s | (s) = 0}. Also note that f^-1(1) = 0, f^-1(δ± i√(1-δ^2)) = ± it where t = √(1-δ)/√(1+δ). Therefore, D(s) satisfies the following properties:
	
		
  * For |s| < 1, D(s) = 0 iff s = 0.
		
  * For |s| < 1, D(s) = 1 iff s = ± it.
		
  * |D(s)| < ∞ for |s| < 1.
	

	Note that the last holds by the quasi-stability of z(s). Since z(s) = 0 implies (s) ≤ 0, D(s) = ∞ implies |s| ≥ 1. In particular, the roots of x, y, z that have 0 real part now correspond to points |s| = 1 such that D(s) = 1, 0, ∞ respectively. For any ϵ > 0, let
	
    D_ϵ(s) := D(s/1+ϵ).

	D_ϵ(s) then satisfies
	
		
  * For |s| ≤ 1, D_ϵ(s) = 0 iff s = 0.
		
  * For |s| ≤ 1, D_ϵ(s) = 1 iff s = ± i(1+ϵ)t.
		
  * |D(s)| < ∞ for |s| ≤ 1.
	
	Precomposing with the inverse fractional linear transformation f^-1(s) = (s-1)/(s+1), we get
	
    R_ϵ(s) := D_ϵ(s-1s+1).

	By the properties of D_ϵ(s) above, we find that R_ϵ(s) satisfies
	
		
  * For (s) ≥ 0, R_ϵ(s) = 0 iff s = 1.
		
  * For (s) ≥ 0, R_ϵ(s) = 1 iff s = δ_ϵ± i√(1-δ_ϵ^2) where
		
    δ_ϵ = 1-(1+ϵ)^2t^21+(1+ϵ^2)t^2.

		
  * For (s) ≥ 0, |R_ϵ(s)| < ∞.
	

	Moreover, R_ϵ(s) ≠ 0, 1, ∞ for any s such that (s) < 0. We can rewrite R_ϵ(s) as R_ϵ(s) = p(s)/q(s).
	Note that by the first property of R_ϵ, the only root of p(s) in {s | (s) ≥ 0} is at s = 1. By properties of f(s), f^-1(s), one can show that p(-1) = 0. This follows from the fact that R(-1) = 0, which implies that lim_s→∞ D(s) = lim_s→∞D_ϵ(s) = 0, and therefore R_ϵ(-1)  = 0. Therefore, p(s) = (s^2-1)y_ϵ(s) where y_ϵ(s) has no roots in {s | (s) ≥ 0}.
	By the second property of R_ϵ, the only roots of q-p in {s | (s) ≥ 0} are at ±δ_ϵ + i√(1-δ_ϵ^2). Therefore, q-p = (s^2-2δ_ϵ s+1)x_ϵ(s) where x_ϵ(s) has no roots in {s | (s) ≥ 0}.
	Finally, by the third property of R_ϵ we find that z_ϵ(s) = (s^2-2δ_ϵ s+1)x_ϵ(s)+(s^2-1)y_ϵ(s) is stable. Moreover, basic properties of fractional linear transformations show that if (x) = n ≥(y) = m, then x_ϵ, y_ϵ are both of degree n. Therefore, x_ϵ, y_ϵ, z_ϵ are stable polynomials satisfying (<ref>) for δ_ϵ. For any δ̂ < δ, we can take ϵ such that δ_ϵ = δ̂, proving the desired result.

	Note that if we start with δ admissible by stable x,y,z of degree at most n, then we can do the reverse of this procedure to perturb x,y,z to quasi-stable x̂, ŷ, ẑ. By the reverse of the arguments above, x̂, ŷ, ẑ will be quasi-stable but at least one of these polynomials will not be stable. These polynomials will be associated to some quasi-admissible δ̂ > δ. This gives the proof of Theorem <ref>.

	The proof above describes the following algorithm for perturbing quasi-stable x,y,z satisfying (<ref>) to obtain stable x̂, ŷ,ẑ satisfying (<ref>).

	

	Input: Real numbers δ, ϵ > 0 and real polynomials x,y,z ∈H satisfying (<ref>).

	Output: δ̂ and real polynomials x̂,ŷ,ẑ∈ H satisfying (<ref>).
	
		
  * Let R(s) = (s^2-1)y(s)/z(s). For ϵ > 0, compute
		
    R_ϵ(s) = R((2+ϵ)s + ϵϵ s + (2+ϵ)).

		
  * Reduce R_ϵ(s) to lowest terms. Suppose that in lowest terms R_ϵ(s) = p(s)/q(s).
		
  * Factor p(s) as (s^2-1)ŷ(s) and factor q(s)-p(s) as (s^2-2δ̂s+1)x̂(s). Let ẑ(s) = q(s).
	

	

	To further illustrate the method of algebraic specification and this algorithm for perturbing to get quasi-stable polynomials, we give the following detailed example.

	Say we are interested in x of degree 4. We may then give the following algebraic specification of x, y, z discussed in Section <ref>. In the shorthand of (<ref>), this is the configuration [1],[],[].
	
    x(s) = (s^2+2δ s+1)(s^2+A)
    
    	y(s) = k
    
    	z(s) = s^6

	As in Section <ref>, we solve (s^2-2δ s+1)x(s) + (s^2-1)y(s) = z(s). This implies that δ, A, k satisfy 16δ^4-20δ^2+5 = 0, A = 4δ^2-2, k = 4δ^2 - 2. Taking the largest root of 16δ^4-20δ^2+5 gives δ = √(10+2√(5))/4, A = k = (√(5)+1)/2. Given numerically to six decimal places, δ = 0.951057. Computing R(s) using exact arithmetic, we get
	
    R(s) = (s^2-1)y(s)z(s) = (s^2-1)(√(5)+1)2s^6

	We then use a fractional linear transformation s ↦ (1+s)/(1-s) to get:
	
    D(s)    = R((1+s)/(1-s))
       = 2s(√(5)+1)(s-1)^4s^6+6s^5+15s^4+20s^3+15s^2+6s+1


	One can verify that D(s) can equal 1 on the boundary of the unit circle, so we push these away from the boundary (with ϵ = 0.01) by defining
	
    D_ϵ(s)    =D(s/1+0.01)
       =6.40805(0.99010s-1)^4s0.942045s^6+…+5.94054s


	While we gave an approximate decimal form above for brevity, this computation can and should be done with exact arithmetic. We let R_ϵ(s) = f_ϵ((s-1)/(s+1)). Writing R_ϵ(s) as p(s)/q(s) in lowest terms, we get:
	
    p(s) = 64080.55401(0.990990s+199.00990)^4(s^2-1)
    
    	q(s) = 0.62122× 10^14s^6 + … +0.94204


	As proved above, p(s) will equal (s^2-1)ŷ(s). Dividing p(s) by the s^2-1 factor, we get a polynomial ŷ(s) such that its only root is at s = -201. Therefore ŷ(s) is stable. The denominator, ẑ(s) is easily verified to only have roots with negative part. Finally, the polynomial q(s) - p(s) will equal (s^2-2δ̂s+1)x̂(s). Finding its roots, one can show that q(s)-p(s) only has roots with negative real part, except for roots at s = 0.950097 ± 0.311954i. These roots are of the form δ̂±√(δ̂^2-1) for δ̂ = 0.950097. Therefore δ̂ = 0.950097 is admissible via the stable polynomials x̂,ŷ,ẑ. While we have decreased δ slightly, we have achieved stability in the process. By decreasing ϵ, we can get arbitrarily close to our original δ.
	



§ OPTIMALITY OF ALGEBRAIC SPECIFICATION


	Not only does our method of algebraic specification find larger δ than have been found before, one can view previous approaches to the Belgian chocolate problem as approximating algebraic specification. In particular, previously discovered admissible δ can be seen as approximating some quasi-admissible δ' that can be found via algebraic specification.

	For example, in <cit.>, Chang and Sahinidis found that δ = 0.9739744 is admissible by
	
    x(s)    =s^10 + 1.97351109136261s^9
       +5.49402092964662s^8 + 8.78344232801755s^7
       + 11.67256448604672s^6 + 13.95449016040116s^5
       +11.89912895529042s^4 + 9.19112429409894s^3
       +5.75248874640322s^2+2.03055901420484s
       +1.03326203778346,
    
    	y(s)   =0.00066128189295s^5+3.611364710425s^4
       +0.03394722108511s^3+3.86358782861648s^2
       +0.0178174691792s+1.03326203778319.


	The roots of x,y,z were discussed in Section <ref>. As previously noted, x,y,z are close to polynomials with repeated roots on the imaginary axis. Examining the roots of x,y,z, one can see that x,y,z are tending towards quasi-stable polynomials x', y', z' that have the same root structure as the algebraic configuration [3,1],[2],[1]. In other words, we will consider the following quasi-stable polynomials:

	
    x'(s) = (s^2+2δ' s + 1)(s^2+A_1)^3(s^2+A_2)
    
    	y'(s) = k(s^2+B)^2
    
    	z'(s) = s^10(s^2+C)


	Solving for the free parameters and finding the largest real δ' such that A_1, A_2, B, C ≥ 0, we obtain the following values, given to seven decimal places.
	
    δ' = 0.9744993
    
    	A_1 =  1.3010813
    
    	A_2 =  0.4475424
    
    	B =    0.5345301
    
    	C =    2.5521908
    
    	k =    3.4498736.


	One can easily verify that taking these values of the parameters, the roots of x, y, z are close to the roots of x',y',z'. These algebraically designed x', y', z' possess the root structure that x,y,z are tending towards. Moreover, the x', y', z' show that δ' is quasi-stable and their associated δ' gives an upper bound for the δ found by Chang and Sahinidis. This demonstrates that the stable polynomials found by Chang and Sahinidis are tending towards the quasi-stable ones listed above. Moreover, by Theorem <ref> all δ < 0.9744993 are admissible.

	In fact, many examples of admissible δ given in previous work are approximating quasi-admissible δ found via algebraic specification. This includes the previously mentioned examples in <cit.> and all admissible values of δ given by Chang and Sahinidis in <cit.>. We further conjecture that for all admissible δ, there is a quasi-admissible δ' > δ that can be achieved by algebraically specified x,y,z.

	More formally, if we fix x, y to be of degree at most n, let δ_n^* denote the supremum of the optimization problem in (<ref>). Note that as discussed in Section <ref>, δ_n^* is not admissible by x,y of degree at most n. The empirical evidence given in this section and in Sections <ref> and <ref> suggests that this δ_n^* is quasi-admissible and can be obtained through algebraic specification. This leads to the following conjecture.

	For all n, δ_n^* is quasi-admissible by some x,y,z that are formed via algebraic specification.



§ CONCLUSION


	The Belgian chocolate problem has remained resilient to direct global optimization techniques for over a decade. Most prior work attempts to maximize δ subject to the stability constraints by applying iterative methods to complicated non-convex regions. By contrast, we find the largest known value of δ in a more direct fashion. We do this by reducing our problem to combinatorial optimization over a finite set of algebraically constructed limit points. Our key algebraic insight is that quasi-admissible δ are limiting values of the admissible δ. In fact, previous methods actually find admissible δ that approach quasi-admissible δ. We give the method of algebraic specification to design quasi-stable polynomials and directly find these quasi-admissible δ by solving a system of equations. We then show that we can perturb these quasi-stable polynomials to obtain stable polynomials with admissible δ that are arbitrarily close to the quasi-admissible δ. We show that this method recovers the largest admissible δ known to date and gives a much better understanding of the underlying landscape of admissible and quasi-admissible δ. We conjecture that for all n, the supremum of all δ admissible by x,y of degree at most n is a quasi-admissible δ that can be found through our method of algebraic specification.



§ ACKNOWLEDGMENTS


	The authors would like to thank Bob Barmish for his valuable feedback, discussions, and advice. The first author was partially supported by the National Science Foundation grant DMS-1502553. The second author was partially supported by the Simons Foundation grant MSN179747.

spmpsci



