CONSTRAINED GLOBAL OPTIMIZATION OF MULTIVARIATE POLYNOMIALS USING POLYNOMIAL B-SPLINE FORM AND B-SPLINE CONSISTENCY PRUNE APPROACH

. In this paper, we propose basic and improved algorithms based on polynomial B-spline form for constrained global optimization of multivariate polynomial functions. The proposed algorithms are based on a branch-and-bound framework. In improved algorithm we introduce several new ingredients, such as B-spline box consistency and B-spline hull consistency algorithm to prune the search regions and make the search more efficient. The performance of the basic and improved algorithm is tested and compared on set of test problems. The results of the tests show the superiority of the improved algorithm over the basic algorithm in terms of the chosen performance metrics for 7 out-off 11 test problems. We compare optimal value of global minimum obtained using the proposed algorithms with CENSO, GloptiPoly and several state-of-the-art NLP solvers, on set of 11 test problems. The results of the tests show the superiority of the proposed algorithm and CENSO solver (open source solver for global optimization of B-spline constrained problem) in that it always captures the global minimum to the user-specified accuracy.


Introduction
Generally constrained global optimization of nonlinear programming problems (NLP) is the study of how to find the best (optimum) solution to a problem. The constrained global optimization of NLPs is stated as follows. Branch-and-bound framework is commonly used for solving constrained global optimization problems [13,17]. For instance, several interval methods [14,18,19,35] use this framework to find the global minimum of NLPs. Since interval analysis methods require function evaluation, which leads to a computationally slow algorithm. Compared with these methods the global optimization algorithm of multivariate polynomial using Bernstein form, e.g. [26,27] has the advantage that it avoids function evaluations which might be costly if the degree of the polynomial is high. Global optimization of polynomials using the Bernstein approach needs transformation of the given multivariate polynomial from its power form into its Bernstein form. The minimum and maximum values of the Bernstein coefficients provide lower and upper bounds for the range of polynomial. Generally, this range enclosure (i.e. bounds) obtained is overestimated in nature and can be improved by degree elevation. Unfortunately, this process implies the increase of computation time.
In this paper we propose polynomial B-spline as an inclusion function [7,9,11,25]. The minimum and maximum B-spline coefficients provide lower and upper bounds for the range of polynomial. The range enclosure obtained using the polynomial B-spline form can be sharpened by increasing the number of B-spline segments, i.e. without degree elevation as shown in Figure 1, which motivates us to use polynomial B-spline form as an inclusion function. In the B-spline approach for unconstrained global optimization [29] a B-spline is used to approximate the objective function with randomly scattered data using the least-square and pseudo-inverse methods. The use of B-spline approach for constrained global optimization is given in [11,13] and references therein. Strength of the B-spline form, and thus solvers operating on the B-spline form, is the possibility for exact representation of multivariate (piecewise) polynomials and approximate representation of any function by sampling. Whereas we follow the procedure given in Section 2 by [21,22] to obtain the B-spline representation of a multivariate polynomial. This procedure do not require sample points and corresponding function evaluations. To the best of our knowledge, there are few papers in the literature on B-spline constrained global optimization [11,13].
In this work, we propose B-spline based algorithms for solving non-convex nonlinear multivariate polynomial programming problems, where the objective function and constraints ( & ℎ ) are limited to be polynomial functions. The proposed work extends the Bernstein method in [26,27] and B-spline method in [11][12][13] for constrained NLPs. The extensions are based on tools such as B-spline hull consistency (BsHC) and B-spline box consistency (BsBC) to contract the variable domains. The merits of the proposed approach are: (i) it avoids evaluation of the objective function and constraints; (ii) an initial guess to start optimization is not required, only an initial search box bounding the region of interest; (iii) it guarantees that the global minimum is found to a user-specified accuracy, and (iv) prior knowledge of stationary points is not required. Numerical performance of the proposed basic and improved algorithms are tested on 11 standard benchmark problems taken from [1,4,5,10,16] with dimensions varying from 2 to 7 and the number of constraints varying from 1 to 11. The optimal value of global minimum obtained with the proposed algorithms for 3 test problems are also compared with CENSO, GloptiPoly and some of the well-known NLP solvers. The rest of the paper is organized as follows. In Section 2, we give the notations and definitions of the B-spline form. In Section 3, we present the Bspline constrained global optimization and outline range enclosure property, subdivision procedure, the cut-off test, and the basic algorithms. In Section 4, we present the B-spline hull, B-spline box consistency techniques, and the B-spline constrained global optimization. In Section 5, we first compare the performances of the proposed basic and improved algorithms, and then we compare the optimal value of global minimum obtained using proposed algorithm with GloptiPoly [15] and several state-of-the-art NLP solvers. We give the conclusion of the work in Section 6.

Univariate case
Firstly, we consider a univariate polynomial to be expressed in terms of the B-spline basis of the space of polynomial splines of degree ≥ (i.e. order + 1). By substituting (2.1) into (2.7), we get

Multi segment B-splines
Let us consider a polynomial ( ) = 3.371 3 − 10.10 2 + 8.233 + 2 and ∈ [0, 2]. Its polynomial B-spline plot with number of segments equal to 1, 2 and 3 are shown in Figure 1. The B-spline with a single segment has four control points, while the one with two segments has five control points, and the one with three segments has six control points. The advantage of B-spline with more number of segments is that we have more control points. This gives a tight range enclosure without having to increase the degree of the B-spline. The drawback of having more segments is the increase in computation time with the number of B-spline coefficients. In our application to global minimization, we find that the B-spline having a number of segment equal to the order of B-spline plus one is a good option.
(1) equal to order of B-spline ( = + 1): Let us continue considering same polynomial ( ) as above. Its polynomial B-spline form of degree = 3 and order 4 with the number of segments taken equal to order of B-spline, will consist of seven B-spline coefficients, i.e. seven B-spline control points and seven B-spline basis functions. The plot of these seven basis functions is shown in Figure 2. As seen from this figure, one of the B-spline basis function that is 3 3 lies on the entire domain of .  (2) taken equal to order of B-spline plus one ( = + 2): We continue with the same polynomial ( ) as above and obtain the polynomial B-spline form with the number of segments equal to order plus one. Now, the B-spline has eight B-spline coefficients and eight B-spline basis functions. As shown in Figure 3, half the B-spline basis functions are having the support of lower domain value of , whereas the other half has the support of upper domain value of . Because of this symmetry, a B-spline with the number of segments equal to order plus one is a good option for our application of global minimization.

Multivariate case
Now, we derive the B-spline representation of a given multivariate polynomial . . . ∑︁  where := ( 1 , 2 , . . . , ), and := ( 1 , 2 , . . . , ). By substituting (2.1) for each , equation (2.10) can be written as  The B-spline form of a multivariate polynomial is defined by (2.11). The partial derivative of a polynomial in a particular direction can be found from the B-spline coefficients of the original polynomial on a box b ⊆ x, the first partial derivative with respect to of a polynomial ( ) in B-spline form in [34] where represents knot vector of . Now, ′ (b) contains an enclosure of the range of the partial derivative of on b.

Example
We consider following example to explain the above ideas.
Example 2.1. Let ( , ) = 2 + 2 − − + 0.34 and , ∈ [0.5, 1.5]. We want to obtain the polynomial B-spline form having two B-spline segments for given power form polynomial. The matrix form representation of ( , ) is The degree of variable in the given polynomial is 1 = 2 and that of is 2 = 2. The degree of B-spline in direction 1 , can be greater or equal to the degree of . Similarly, the B-spline in direction will have a degree 2 equal to or greater than degree of variable . In practice, we can therefore take 1 = 1 = 2 and 2 = 2 = 2. Therefore, the order, of B-spline will be = + 1 = 3 in both the directions. As = 2, the B-spline will have two segments in each direction. As , ∈ [0.5, 1.5], the knot vector for both the variables will be the same.
From (2.11), B-spline representation of ( , ) can be expressed in matrix form as From (2.13), we calculate the values of B-spline coefficients as A simple computation leads to the B-spline coefficient matrix

B-spline constrained global optimization
Let ∈ N be the number of variables and = ( 1 , 2 , . . . , ) ∈ R . A multi-index is defined as = where IR denotes the set of compact intervals. Let wid x denotes the width of x , that is wid Global optimization of polynomials using the polynomial B-spline approach needs transformation of the given multivariate polynomial from its power form into its polynomial B-spline form. Then B-spline coefficients are collected in an array (x) = ( (x)) ∈ , where = { : ≤ }. This array is called a patch. We denote 0 as a special subset of the index set comprising indices of the vertices of this array, that is (3.1)

Range enclosure property
The following lemma describes the range enclosure property of the B-spline coefficients.
Obtaining the B-spline coefficients of multivariate polynomials by transforming the polynomial from power form to B-spline form provides an enclosure of the range of the multivariate polynomial on x. Then by Lemma 3.1, the minimum and the maximum values of B-spline coefficient provide lower and upper bounds for the range of polynomial . This range enclosure will be sharp if and only if min( (x)) ∈ (respectively max( (x)) ∈ ) is attained at the indices of the vertices of the array (x), as described in Lemma 3.2. This condition is known as the vertex condition.
Based on Bernstein coefficients range enclosure proofs [32], the proof of Lemma 3.1 is given next.
Proof. Let = + 1 represent the set of integers between and ( < ). Also let be the B-spline representation for polynomial and let¯be its range, that is min Proof. Let be the B-spline representation for polynomial and let¯be its range, We first note that (0) = 1 , We have that (0) = 1 = = + = (1) and the property is valid. and assume that min for simplicity. Then the bound is sharp, i.e.
the same argument is valid if min Definition 3.3. The vertex condition is said to be met within a given tolerance , if As said earlier set 0 comprises of indices of the vertices of the B-spline coefficient array (x), where min gives the minimum value of B-spline coefficient at any vertex point. If the difference between the minimum value of B-spline coefficient at any vertex and the minimum value in B-spline coefficient array is less than , then vertex condition is said to be met within a given tolerance .

B-spline subdivision procedure
The proposed algorithm is based on branch and bound framework of global optimization. Therefore we need to use domain subdivision. Generally, the range enclosure obtained from Lemma 3.1 is over-estimated and can be improved either by subdividing the domain, degree elevation of the B-spline or by increasing the number of B-spline segments. Subdivision is generally more efficient than degree elevation strategy [6] or increasing the number of B-spline segments. Therefore subdivision strategy is preferred over the latter two. A subdivision in the th direction (1 ≤ ≤ ) is a bisection perpendicular to this direction. Let be any subbox. Further suppose that x is bisected along the th component direction then, two subboxes x and x are generated as

The cut-off test
As mentioned earlier, the minimum and maximum B-spline coefficients provide the range enclosure of the function. Let˜be the current minimum estimate, and {b, (b)} be the current item for processing. We denote the minimum over the second entry of item {b, (b)} as . If is greater than˜, then this item cannot contain global minimum and it can be discarded. If the maximum over the second entry of item {b, (b)} is lesser thañ , then the current minimum estimate can be updated, and˜takes this maximum value as the new value. Next we present the cut-off test algorithm.

A basic B-spline constrained global optimization algorithm
In this subsection, we present the basic B-spline algorithm for constrained global optimization of multivariate nonlinear polynomials. The algorithm is inspired by the one described in [31,33].
This basic algorithm uses the polynomial coefficients of the objective function, the inequality constraints, and the equality constraints. The inputs to the algorithm are the polynomial degrees and the initial search box, while the outputs are the global minimum and global minimizers. The polynomial degree is used to compute the B-spline segment number, as the B-spline is constructed with number of segments equal to order of the B-spline plus one. As equality constraints ℎ ( ) = 0 are difficult to verify on computers with finite precision, the equality constraints ℎ ( ) = 0 in (1.3) are replaced by relaxed constraints ℎ ( ) ∈ [− zero , zero ], = 1, 2, . . . , , where zero > 0 is a very small number.
The basic algorithm works as follows. We start the algorithm by computing the B-spline segment vectors , and ℎ , where = [ 1 , . . . , ], for each variable occurring in the objective, inequality and equality polynomials. We keep it as + 1 for each variable, giving = [ 1 + 2, . . . , + 2]. Then, we compute the B-spline coefficients of objective, inequality and equality constraint on the initial search box. We store them in arrays (x), (x) and ℎ (x) respectively. We initialize the current minimum estimate˜to the maximum B-spline coefficient of the objective function on x. Next, we initialize a flag vector with each component to zero, a working list ℒ with the item {x, (x), (x), ℎ (x), }, and a solution list ℒ sol to the empty list. We then pick the last item from the list ℒ and delete its entry from ℒ. For this item, we subdivide the box x along the longest width direction creating two subboxes b 1 and b 2 . We compute the B-spline coefficients arrays {b , (b ), (b ), ℎ (b )}, = 1, 2 for b 1 ,b 2 and the B-spline range enclosures D (b ), D (b ) and D ℎ (b ) of objective, inequality and equality constraint polynomials respectively. We check the feasibility of the inequality and equality constraints for b 1 , b 2 using the B-spline coefficients of the constraint polynomials functions by doing the following tests: zero ] for all = 1, 2, . . . , and = 1, 2, . . . , , then b is a feasible box.
-If D (b ) > 0 for some , then b is a infeasible box and can be deleted.      {Set local current minimum estimate}    We then have the following theorem, The following theorem discusses the convergence properties of basic constrained global optimization algorithm. Theorem 3.9. Let basic constrained global optimization algorithm applied to the box x, the inclusion functions , , and of , , and ℎ, respectively. Let the contraction assumption (3.5) and (3.6) and condition ( ) be satisfied. Then the sequence ( ) is nested and ∩ ∞ =1 = D which means that → D. Furthermore˜→ * with˜≤ * and ↘ * as → ∞.
The similar proof for convergence of algorithm for global optimization of B-spline constrained problems is given in [13].

Improved algorithm with B-spline consistency techniques
We can apply the concept of consistency to each constraints of the problem to eliminate subboxes of the given box that cannot contains the solution. Let ( ) = ( 1 , . . . , ) = 0, ∈ x be a constraint that must be satisfied in an optimization problem we use B-spline box and hull consistency to eliminate subboxes of x that cannot contain a point satisfying ( ) = ( 1 , . . . , ) = 0. The B-spline box and hull consistency can be applied for inequalities. Suppose that in place of the equality we have an inequality ( ) ≤ 0. We can replace this inequality by ( ) = [−∞, 0] and obtain the equation ( ) + [0, +∞] = 0. B-spline box consistency requires application of interval Newton method, generally it does not perform well when the variable bound is very wide. Next we present the B-spline box consistency (BsBC) and B-spline hull consistency (BsHC) techniques with an examples which show that B-spline hull consistency gives 60% pruning and B-spline box consistency gives 38% pruning in variable bound. These can be viewed as extensions of the ideas in [27] in the context of the polynomial B-spline based approach to global optimization.

B-spline hull consistency
In this subsection, we present the B-spline hull consistency to reduce the search region. B-spline hull consistency can be viewed as extension of interval hull consistency in the context of the B-spline form. We present next the procedure of interval hull consistency described in [14,30].
Let us start by considering to apply hull consistency for domain reduction of a variable 1 using a two variable equality constraint ℎ( 1 , 2 ) = 0, ℎ( 1 , 2 ) = 0,0 + 1,0 1 + 0,1 2 + 1,1 1 2 + · · · + 0,2 2 2 + · · · + ,0 1 + · · · + , 2 = 0. (4.1) The implementation of the hull consistency involves the constraint inversion. To obtain constraint inversion we select term having only one variable with highest degree. Lets consider ,0 1 then constraint inversion is given as, Then we use B-spline expansion to obtain the range h ′ (x 1 , x 2 ) of constraint inversion function ℎ ′ ( 1 , 2 ). The new value for interval x 1 is estimated in [14] as This procedure is repeated to contract the domain of 2 variable using the same constraint function ℎ( 1 , 2 ) = 0. This domain box is further contracted in a similar way using the remaining equality constraints. The hull consistency method can also be applied to inequality constraints, we need only replace an inequality of the form ( ) ≤ 0 by the equation ( ) = [−∞, 0] in [14]. Then B-spline hull consistency can be applied as before to this equality constraint. We illustrate the B-spline hull consistency method for equality constraint via an example.   To get a new interval for 2 1 , the interval methods require the evaluation of ℎ 1 (x), so computing the B-spline coefficients for ℎ 1 (x) we get Since (x ′ 1 ) 2 must be non-negative The updated value of x 1 is therefore Next we present B-spline hull consistency (BsHC) algorithm.
Algorithm 3 : BsHC ( , , , x, , , ). Input : Here is a cell structure containing the coefficients array of all the constraints, is a cell structure, containing degree vector for all constraints. Where elements of degree vector defines the degree of each variable occurring in all constraints polynomial, is a cell structure containing vectors corresponding to all constraints, i.e. , ℎ . Where elements of this vector define the number of B-spline segments in each variable direction, the number of constraint variables , number of inequality constraints and number of equality constraints . Output: A box x that is contracted using B-spline box consistency technique for the given constraints. Begin Algorithm {Formulate the constraint inverse polynomial} From the coefficient matrix ( ), choose the monomial term in i.e. having the highest degree and substitute zero for its coefficient i.e = 0, then multiply the coefficient matrix ( ) by −1.

(c)
{Compute h ′ } Compute the B-spline coefficients of the constraint inverse polynomial, and then obtain h ′ , as minimum to maximum of these B-spline coefficients. The algorithm in [8] is suggested for the computation. if < + 1 then Compute

B-spline box consistency
First of all, we recall the procedure of interval box consistency described in [14,30]. The implementation of interval box consistency involves the application of a one dimensional Newton method, to solve a single equation for a single variable. Let us start by applying box consistency to the following equality constraint (4.4) We use box consistency to eliminate those subboxes of x that do not satisfy ℎ( ) = 0. If we replace all the variables except th variable by their interval bounds, we obtain the equation If 0 / ∈ ℎ( ) for in some subinterval x ′ of x , then we do not have consistency for ∈ x ′ and the subbox ( If 0 / ∈ h( ), then < so we use a interval Newton method to remove points from the lower end of x that are less than . Thus, the Newton result when expanding about the point is [14] ( The new contracted interval can be obtained by intersecting the interval values of x and . That is, Similarly, the upper bound can be contracted using the right narrowing operation. This procedure is repeated for x , = 1, . . . , , using the other equality constraints. We can also apply box consistency to an inequality constraints, by replacing an inequality of the form ( ) ≤ 0 by the equation ( ) = [−∞, 0] in [14], that is rewrite this an another function ℎ( ) given by Lets consider to apply box consistency for domain reduction of variable x ′ = [ , ] using equality constraint (4.9). To apply B-spline Newton contractor for ℎ(x ) for end point , we required interval range enclosure h( ) of equality constraint (4.9). To obtain h( ) we first find interval range enclosure g( ) for ( ) then the interval range enclosure h( ) is obtained by considering min g( ) and ∞ as lower and upper bound of h( ) respectively. Thus by applying B-spline Newton contractor for ℎ( ) for end point , we obtain )︁ .
When we apply the B-spline box consistency to a selected variable of a multivariate constraint, then only that variable's domain will be contracted. The domains of the remaining variables will remain unaffected. We apply B-spline box consistency to each variable in turn, to get a contracted box in all variables. Suppose ℎ( ) = 0 is the given equality constraint. We compute the B-spline coefficients array (x) for this constraint, and consider a variable direction, say the first one x 1 = [ , ]. In the B-spline box consistency, we try to increase the value of and decrease the value of , thus effectively reducing the width of x 1 . The procedure to increase the value of is listed in the below steps, (1) Compute interval ℎ( ). Corresponding to 1 = , find all B-spline coefficients in (x). The minimum to maximum of these B-spline coefficients gives an interval ℎ( ). As 0 ∈ ℎ( ), means ℎ( ) is not completely positive interval, then we cannot increase . We instead switch over to the other end point and try to decrease it in the same way as we try to increase .
We illustrate the BsBC method for equality constraints via an example.
Note that we are constructing B-spline with + 1 number of segments. Consider the application of BsBC along the first component direction, that is, along 1 . Along the direction 1 , the first row corresponds to 1 = = 0.2, and the sixth row corresponds to 1 = = 1. Along the first row, the minimum to maximum values of the B-spline coefficients, are −0.88 and 0.12 giving h( ) = [−0.88, 0.12]. Along the sixth row, the minimum and maximum values of the B-spline coefficients, respectively are 2 and 3 giving h( ) = [2,3]. Since 0 ∈ h( ) the left end point cannot be increased. However, 0 ̸ ∈ h( ), hence the right end point can be decreased. The partial derivative in the direction 1 , that is, h ′ 1 will be as follows Therefore, the updated value of x 1 is Next we present B-spline box consistency (BsBC) algorithm.
Algorithm 4: BsBC ( , , , x, , , ). Input : Here is a cell structure containing the coefficients array of all the constraints, is a cell structure, containing degree vector for all constraints. Where elements of degree vector defines the degree of each variable occurring in all constraints polynomial, is a cell structure containing vectors corresponding to all constraints, i.e. , ℎ . Where elements of this vector define the number of B-spline segments in each variable direction, the number of constraint variables , number of inequality constraints and number of equality constraints . Output: A box x that is contracted using B-spline box consistency technique for the given constraints. Begin Algorithm {Compute B-spline coefficient for constraint polynomial} ( ) ( ) using ( ( ) , ( ) , ( ) , x). The algorithm in [8] is suggested for the computation.
there is no zero of ℎ in entire interval x , and hence the constraint ℎ is infeasible over box x. Exit the algorithm in this case with x ′ = ∅.

Proposed improved algorithm for constrained global optimization
To apply the proposed B-spline box and B-spline hull consistency algorithms using the basic algorithm (see Sect. 3.4), we modify step 7 of basic algorithm to step * 7 given below. We refer to the resulting modified algorithm as improved algorithm.

Numerical test
In this section, we present the result and analysis of our tests. The computations are done on a PC Intel i3-370M 2.40 GHz processor, 6 GB RAM, while the algorithms are implemented in MATLAB [24]. An accuracy = 10 −6 is prescribed for computing the global minimum and minimizer(s) in each test problem. For the tests, we select 11 benchmark optimization problems taken from [1,4,5,10,16] (described in Appendix A). Table 1 reports the global minimum obtained with the proposed algorithms.

Comparisons between basic and improved proposed algorithms
First, we test and compare the performance of the basic and improved B-spline constrained global optimization algorithms on a set of 11 test problems. The performance metrics are taken as the number of subdivisions and computation time (in seconds) required to compute the global minimum for the problem. These values are reported in Tables 2 and 3 respectively. For these metrics, we give the computed as where PMBA is the performance metric with the basic algorithm, and PMIA is the performance metric with the improved algorithm. In Table 2, we compare the number of subdivisions required for obtaining global minima using the basic and improved algorithms. It is found that except problem P2 for all test problems the number of subdivisions required to compute the global minimum reduces in the range of 2.22% to 96%. Whereas algorithm is unable to solve problem P10. In Table 3, we compare the computation time required for obtaining the global minima using the basic and improved algorithm, it is observed that the improved algorithm computes the results slower than that of the basic algorithm except for the problems number 3, 6, 7 and 9. The use of B-spline box and B-spline hull consistency algorithms found to be very effective in pruning the search domain by discarding those subboxes that cannot contain global minimizers. As shown in Section 4.3 during each iteration improved algorithm requires to compute B-spline coefficients more than once (like steps * 7( ) and * 7( ) in improved algorithm) as compared to basic algorithm. In improved algorithm during each iteration B-spline coefficients are computed after domain pruning by the B-spline box and B-spline hull consistency. Also, in each iteration of B-spline hull consistency algorithm requires to compute B-spline coefficients after each constraint inversion (cf. 4.1) to contract the variable bounds. All the variable bounds are contracted using each constraint function individually. This is also evident from the results in the Table 3 for an improved algorithm. We would like to mention that, a more sophisticated implementation of B-spline box and hull consistencies can alleviate the computational bottle neck associated with them. In case of problem P2 algorithm is unable to minimize the number of subdivision still it is slow because B-spline box consistency (Newton method) is fast when the initial intervals are small. Unfortunately, the running time of algorithm increases linearly with the size of the initial interval [36].

Comparison with GloptiPoly and other NLP solvers
Next, we compare the optimal value of global minimum found using the proposed B-spline algorithms ( & ) with CENSO, GloptiPoly [16,20] and BARON, LINDO Global and SCIP NLP solvers. For the current tests, we consider the above set of 11 test functions. The idea is to investigate where same optimal value of global minimum can be found with the considered methods. The CENSO solver implements the algorithm in [13] and is available on Github: https://github.com/bgrimstad/censo. The GAMS and GloptiPoly source code for these problems are available at [https://bit.ly/2YHonrh].The GAMS interface for BARON, LINDO Global and SCIP is available through the NEOS server [28]. Table 4 report only those test functions from the above set of 11 test functions for which the proposed algorithms provides the optimal value of global minimum compare to CENSO, GloptiPoly and NLP solvers. The bold values in the table indicate the local minimum value. For P8 and STP1 BARON, LINDOGlobal and SCIP are able to capture the local minimum, where CENSO, GloptiPoly and proposed algorithms found the Notes. (⋆) Indicates that the algorithm did not give the result even after one hour and therefore terminated. Notes. (#) Indicates that the algorithm did not give the result even after one hour and therefore terminated. correct global minimum. In P10, BARON, LINDOGlobal and SCIP are able to capture the local minimum and GloptiPoly did not give the result even after an hour and is therefore terminated. In P8 and STP1 (STP1 is a problem constructed by authors based on P8), the relaxation order for GloptiPoly had to be systematically increased to a high order to obtain convergence to the final results. The proposed algorithms and CENSO find the optimal value of global minimum for all the test problem considered. Where as BARON and LINDO Global are unable to capture optimal value of global minimum in all the 3 test functions even after specifying a large value of the AbsConFeasTol option in GAMS. The SCIP solver is unable to find the optimal value of global minimum for P10 test problem.
For GloptiPoly, the relaxation order need to be systematically increased to a higher order to obtain convergence to the optimal value of global minimum. As the dimension and number of constraints in problem increases, to obtain convergence the relaxation order is gradually increased to a value greater than or equal to the dimension of the problem. For medium dimension problem (like 7 with 11 constraints) GloptiPoly exhausts with memory even with small relaxation order. Whereas other NLP solvers may be able to solve large dimension problems with non-optimal value of global minimum. Here, we want to mention that proposed algorithms are implemented in MATLAB, whereas except GloptiPoly all other NLP solvers mentioned above are implemented in different languages and therefore comparison based on computation time between proposed algorithms and these NLP solvers is not carried out.

Conclusion
We presented the basic and improved algorithm for constrained global optimization of multivariate polynomials using polynomial B-spline form. The performance of the basic and improved algorithm are tested and compared on 11 test problems. The test problems had dimensions ranging from 2 to 7 and number of constraints varying from 1 to 11. The results of the test show that improved algorithm is more efficient, in terms of a number of subdivisions with little extra computational time. We also compared the optimal value of global minimum obtained with the proposed algorithms with CENSO, GloptiPoly and some of the well known NLP solvers, on a set of 3 test problems. The result shows the superiority of the proposed algorithm and CENSO solver, in that it captures the global minimum to user specified accuracy. One possible extension of this research work is to investigate the performance of proposed B-spline approach for solving problems with more number of variables. This problem will be addressed in a future work.