A Class of Continuous Separable Nonlinear Multidimensional Knapsack Problems

Bin Zhang¹, Zhe Lin¹, Yu Wang²

¹Lingnan College, Sun Yat-sen University, Guangzhou, China.
²School of Economic and Business Administration, Chongqing University, Chongqing, China.

DOI: 10.4236/ajor.2018.84015                

Abstract

The nonlinear multidimensional knapsack problem is defined as the minimization of a convex function with multiple linear constraints. The methods developed for nonlinear multidimensional programming problems are often applied to solve the nonlinear multidimensional knapsack problems, but they are inefficient or limited since most of them do not exploit the characteristics of the knapsack problems. In this paper, by establishing structural properties of the continuous separable nonlinear multidimensional knapsack problem, we develop a multi-tier binary solution method for solving the continuous nonlinear multidimensional knapsack problems with general structure. The computational complexity is polynomial in the number of variables. We presented two examples to illustrate the general application of our method and we used statistical results to show the effectiveness of our method.

Keywords

Nonlinear Programming, Convex Programming, Multidimensional Knapsack, Separable Knapsack, Lagrangian Relaxation

1. Introduction

The nonlinear multidimensional knapsack problem is defined as minimizing a convex function with multiple linear constraints. The nonlinear knapsack problem is a class of nonlinear programming, and some methods designed for nonlinear programming can be applied for solving the nonlinear multidimensional knapsack problems. The general nonlinear programming problems have been intensively studied in the last decades, and some different methods have been developed, such as Newton method [1] [2] [3], branch and bound method [4] [5], interior point method [6] [7], sequential quadratic programming method [8] [9] and the filter method [10] [11]. These methods are designed for nonlinear programming problems, and some of them are inefficient or limited for solving the nonlinear knapsack problems since they do not consider the characteristics of the knapsack problems.

Generally, it is much faster and more reliable to solve knapsack problems with specialized methods than with standard methods [12]. Many researchers studied the solution methods for the nonlinear knapsack problems based on the specialized knapsack structures. Most of the research studied the problems with single constraint. Two basic specialized methods are mainly applied for solving the single-constraint nonlinear knapsack problem. One is the multiplier search method [13], and another is the pegging method [14] [15]. Recently, some new methods are proposed for efficiently solving the single-constraint nonlinear knapsack problem. Zhang and Hua developed a united method for solving a class of continuous separable nonlinear knapsack problems [16]. Kiwiel developed the breakpoint searching method for the continuous quadratic knapsack problem [17]. Sharkey et al. studied a general class of nonlinear non-separable continuous knapsack problem [18].

Most research of nonlinear knapsack problems studied the one-dimensional problems with continuous or integer variables, and the proposed methods cannot be directly extended for solving multi-dimensional problems. Some researchers attempted to solve multi-dimensional problems with integer-valued variables. Morin and Marsten firstly studied the nonlinear multidimensional knapsack problems and developed the imbedded state space approach [19]. Some researchers investigate the efficiency of other methods, such as smart greedy method [20], cut method [21] [22], branch and bound method [23] and branch and cut method [24]. Other research studied different applications of multidimensional knapsack, e.g., multi-product newsvendor problems with multiple constraints [10] [25] [26] [27]. The continuous separable nonlinear multidimensional knapsack problems with general structure have not been well studied due to its complexity, and the specialized methods are very limited.

This paper establishes some structural properties of the continuous separable nonlinear multidimensional knapsack problem, and develops a multi-tier binary solution method for solving a class of continuous nonlinear multidimensional knapsack problems with general structure. The computational complexity is polynomial in the number of variables. We presented two examples to illustrate the application of our method, and the statistical study with the randomly generated instances for different problem sizes are reported to show the effectiveness of our method.

The paper is organized as follows. In Section 2, the nonlinear multidimensional knapsack problem is described. Section 3 studies the structural properties of the problem, and develops the algorithm. Section 4 presents the illustrative examples and the statistical results. Finally, the concluding remarks are given in Section 5. All proofs are listed in Appendix.

2. Problem Formulation

The continuous separable nonlinear multidimensional knapsack problem studied in this paper is as follows (denoted as problem P):

$\text{Min}f\left(x\right)={\displaystyle \underset{i=1}{\overset{N}{\sum }}{f}_i\left(x_i\right)}$ , (1)

Subject to

${\displaystyle \underset{i=1}{\overset{N}{\sum }}{c}_{i,j}{x}_i\le C_j,j=1,\cdots ,M},$ (2)

$l_i\le x_i\le u_i,i=1,\cdots ,N$. (3)

The notation used in this paper is listed in Table 1.

In problem P, all objective functions $f_i\left(x_i\right),i=1,\cdots ,n$ are convex and differentiable, the unit resource coefficient $c_{i,j}>0$ for all $i=1,\cdots ,N,j=1,\cdots ,M$ , the resource constraints $C_j>0$ for all $j=1,\cdots ,M$ , and the lower and upper bounds satisfy $0\le l_i<u_i$ for all $i=1,\cdots ,N$.

Since the objective functions and the feasible domain in problem P are all convex, the optimality condition for problem P can be characterized using KKT conditions. Let $\lambda =\left({\lambda }_1,\cdots ,{\lambda }_M\right)$ , ${\lambda }_j\ge 0,j=1,\cdots ,M$ , be the Lagrange multiplier vector for the constraints given in Equation (2), and $w=\left(w_1,\cdots ,w_N\right)$ , $w_i\ge 0,i=1,\cdots ,N$ , $v=\left(v_1,\cdots ,v_N\right)$ , $v_i\ge 0,i=1,\cdots ,N$ be the Lagrange multiplier vectors for the constraints in Equation (3). Thus, the Lagrange function for problem P can be written as:

$L\left(x,\lambda ,w,\nu \right)={\displaystyle \underset{i=1}{\overset{N}{\sum }}{f}_i\left(x_i\right)}-{\displaystyle \underset{j=1}{\overset{M}{\sum }}{\lambda }_j\left(C_j-{\displaystyle \underset{i=1}{\overset{N}{\sum }}{c}_{i,j}{x}_i}\right)}-{\displaystyle \underset{i=1}{\overset{N}{\sum }}{w}_i\left(x_i-l_i\right)}+{\displaystyle \underset{i=1}{\overset{N}{\sum }}{v}_i\left(x_i-u_i\right)}$. (4)

Table 1. Notation.

Let $g_i\left(x_i\right)=\text{d}{f}_i\left(x_i\right)/\text{d}{x}_i$ , $i=1,\cdots ,N$. The KKT conditions for problem P can be summarized as the following proposition.

Proposition 1: The KKT conditions for problem P are:

$g_i\left(x_i\right)+{\displaystyle \underset{j=1}{\overset{M}{\sum }}{\lambda }_j{c}_{i,j}-w_i+v_i}=0,i=1,\cdots ,N,$ (5)

${\displaystyle \underset{i=1}{\overset{N}{\sum }}{w}_i\left(x_i-l_i\right)}+{\displaystyle \underset{i=1}{\overset{N}{\sum }}{v}_i\left(x_i-u_i\right)}=0,$ (6)

${\lambda }_j\left({\displaystyle \underset{i=1}{\overset{N}{\sum }}{c}_{i,j}{x}_i-C_j}\right)=0,j=1,\cdots ,M.$ (7)

Since $f_i\left(x_i\right)$ is convex in $x_i$ , $g_i\left(x_i\right)$ is an increasing function of $x_i$. Let ${\bar{x}}_i$ be the point that satisfies $g_i\left(x_i\right)=0$ if $g_i\left(0\right)\le 0$ and $\underset{x_i\to +\infty }{\lim }{g}_i\left(x_i\right)\ge 0$. If $g_i\left(0\right)>0$ , we let ${\bar{x}}_i=0$. If $\underset{x_i\to +\infty }{\lim }{g}_i\left(x_i\right)<0$ , we set ${\bar{x}}_i=+\infty$. Then ${\bar{x}}_i$ is the optimal solution to the objective function in Equation (1) without any constraint. We summarize it as

$\begin{array}{l}{\bar{x}}_i=\arg \min \left\{f_i\left(x_i\right),0\le x_i\le +\infty \right\}\\ =\{\begin{array}{l}0,\text{if}{g}_i\left(0\right)>0,\\ \arg \left\{x_i|g_i\left(x_i\right)=0\right\},\text{if}{g}_i\left(0\right)\le 0\text{and}\underset{x_i\to +\infty }{\lim }{g}_i\left(x_i\right)\ge 0,\\ +\infty ,\text{if}\underset{x_i\to +\infty }{\lim }{g}_i\left(x_i\right)<0.\end{array}\end{array}$ (8)

3. Structural Properties and Solution Method

In this section, we first investigate the structural properties of the optimal solution to problem P. Then we develop a solution method based on the structural properties for solving problem P.

3.1. Structural Properties

We denote by problem PR the knapsack relaxation problem from problem P, in which the constraints in Equation (2) are relaxed. This implies that we do not consider Equation (2) in problem PR. By analyzing the solution to problem PR, we can find the way to construct the solution to problem P. We let ${\hat{x}}_i$ ( $i=1,\cdots ,N$ ) be the optimal solution to problem PR, then ${\hat{x}}_i$ ( $i=1,\cdots ,N$ ) has the following property.

Proposition 2: The optimal solution to problem PR is ${\hat{x}}_i=\min \left\{\max \left\{{\bar{x}}_i,l_i\right\},u_i\right\}$.

If ${\displaystyle \underset{i=1}{\overset{N}{\sum }}{c}_{i,j}{\hat{x}}_i\le C_j}$ holds for some $j=1,\cdots ,M$ , then the corresponding constraints in problem P are inactive, which can be removed from problem P. In the following, without loss of generality, we assume that ${\displaystyle \underset{i=1}{\overset{N}{\sum }}{c}_{i,j}{\hat{x}}_i>C_j}$ for all $j=1,\cdots ,M$. The KKT conditions in Equation (7) are met at either ${\lambda }_j=0$ , or ${\displaystyle \underset{i=1}{\overset{N}{\sum }}{c}_{i,j}{x}_i=C_j}$. The condition ${\lambda }_j=0$ implies that there is enough resource j at the optimal solution, and hence the j-th constraint is inactive. ${\displaystyle \underset{i=1}{\overset{N}{\sum }}{c}_{i,j}{x}_i=C_j}$

means that the j-th constraint is active, and knapsack space of the j-th constraint must be fully utilized at the optimal solution.

We denote by $x^{*}$ the optimal solution to problem P and ${\lambda }^{*}$ the corresponding Lagrange multiplier vector. Let $x_i\left(\lambda \right)$ be a solution of the KKT conditions in Equation (5) and Equation (6). We denote by $h_i(\cdot )=g_i^{-1}(\cdot )$ , then we have the following proposition.

Proposition 3. (a) $x_i\left(\lambda \right)=\min \left\{\max \left\{h_i\left(-{\displaystyle \underset{j=1}{\overset{M}{\sum }}{\lambda }_j{c}_{i,j}}\right),l_i\right\},u_i\right\}$ , $i=1,\cdots ,N$.

(b) If $\left(x\left(\lambda \right),\lambda \right)$ satisfies ${\lambda }_j=0$ or ${\displaystyle \underset{i=1}{\overset{N}{\sum }}{c}_{i,j}{x}_i=C_j}$ , $j=1,\cdots ,M$ , then we have $x^{*}=x\left(\lambda \right)$.

For any given ${\lambda }_M\ge 0$ , we let $x\left({\lambda }_M\right)$ and ${\lambda }_1,\cdots ,{\lambda }_{M-1}$ be the optimal so-

lution of Equations (5) and (6) and ${\lambda }_j\left({\displaystyle \underset{i=1}{\overset{N}{\sum }}{c}_{i,j}{x}_i-C_j}\right)=0$ , $j=1,\cdots ,M-1$. For

ease of exposition, we denote problem P as $P\left(f,M\right)$ , where $f=\left(f_i,\cdots ,f_N\right)$ is the objective function vector. Problem $P\left(\hat{f}\left({\lambda }_M\right),M-1\right)$ with ${\hat{f}}_i\left({\lambda }_M\right)=f_i+{\lambda }_M{c}_{i,M}{x}_i$ , $i=1,\cdots ,N$ , is an $M-1$ constraint problem with the objective function ${\hat{f}}_i\left({\lambda }_M\right)$ and the first $M-1$ knapsack constraints of problem P.

By analyzing the structural properties of $x\left({\lambda }_M\right)$ and $P\left(\hat{f}\left({\lambda }_M\right),M-1\right)$ , we can prove the following proposition.

Proposition 4. (a) If $\left(x\left({\lambda }_M\right),{\lambda }_M\right)$ satisfies ${\lambda }_M=0$ or ${\displaystyle \underset{i=1}{\overset{N}{\sum }}{c}_{i,M}{x}_i\left({\lambda }_M\right)=C_M}$ , then we have $x^{*}=x\left({\lambda }_M\right)$.

(b) $x\left({\lambda }_M\right)$ is the optimal solution to problem $P\left(\hat{f}\left({\lambda }_M\right),M-1\right)$ with ${\hat{f}}_i\left({\lambda }_M\right)=f_i+{\lambda }_M{c}_{i,M}{x}_i$ , $i=1,\cdots ,N$.

From Proposition 4(a), we know that the optimal solution to problem $P\left(\hat{f}\left({\lambda }_M\right),M-1\right)$ is obtained in two possible cases: 1) ${\lambda }_M=0$ , which means

that the constraint ${\displaystyle \underset{i=1}{\overset{N}{\sum }}{c}_{i,M}{x}_i\left({\lambda }_M\right)\le C_M}$ is not binding and it can be removed

from problem $P\left(f,M\right)$. Therefore, $x^{*}$ can be obtained by solving problem $P\left(f,M-1\right)$ , which has the same structure as problem $P\left(f,M\right)$ ; 2)

${\displaystyle \underset{i=1}{\overset{N}{\sum }}{c}_{i,M}{x}_i\left({\lambda }_M\right)=C_M}$ , which implies that ${\displaystyle \underset{i=1}{\overset{N}{\sum }}{c}_{i,M}{x}_i\left({\lambda }_M\right)\le C_M}$ is an active constraint, and the optimal solution must be obtained at ${\displaystyle \underset{i=1}{\overset{N}{\sum }}{c}_{i,M}{x}_i\left({\lambda }_M\right)=C_M}$ with ${\lambda }_M>0$.

Since problem $P\left(f,M\right)$ can be solved by solving problem $P\left(f,M-1\right)$ in the case of ${\lambda }_M=0$. In the following, we study the case of ${\lambda }_M>0$. Proposition 4(b) indicates that problem $P\left(\hat{f}\left({\lambda }_M\right),M-1\right)$ determines the optimal values of $x\left({\lambda }_M\right)$ and ${\lambda }_j$ , $j=1,\cdots ,M-1$. For any ${\lambda }_M>0$ , the $M-1$ resource constraints could be active or inactive, and the N decision variables could take bound values or non-bound values.

If ${\lambda }_j>0$ , $j=1,\cdots ,M-1$ , constraint j will be active, thus we denote by $J\left({\lambda }_M\right)=\left\{j|{\lambda }_j>0,j=1,\cdots ,M\right\}$ the active constraint set for the given ${\lambda }_M$. Note that $J\left({\lambda }_M\right)$ includes at least one active constraint for the case of ${\lambda }_M>0$.

From Equation (5), we know $x_i\left({\lambda }_M\right)>l_i$ if $-g_i\left(l_i\right)-{\displaystyle \underset{j=1}{\overset{M}{\sum }}{\lambda }_j{c}_{i,j}}>0$ , $i=1,\cdots ,N$ , and $x_i\left({\lambda }_M\right)<u_i$ if $-g_i\left(u_i\right)-{\displaystyle \underset{j=1}{\overset{M}{\sum }}{\lambda }_j{c}_{i,j}}>0$ , $i=1,\cdots ,N$. For the given ${\lambda }_M$ , we define the non-bound variable set $I\left({\lambda }_M\right)$ , and lower and upper bound variable sets $I_L\left({\lambda }_M\right)$ and $I_U\left({\lambda }_M\right)$ as

$I\left({\lambda }_M\right)=\left\{i|g_i\left(l_i\right)<-{\displaystyle \underset{j=1}{\overset{M}{\sum }}{\lambda }_j{c}_{i,j}<g_i\left(u_i\right)},i=1,\cdots ,N\right\}$ , (9)

$I_L\left({\lambda }_M\right)=\left\{i|-{\displaystyle \underset{j=1}{\overset{M}{\sum }}{\lambda }_j{c}_{i,j}\le g_i\left(l_i\right)},i=1,\cdots ,N\right\}$ , (10)

$I_U\left({\lambda }_M\right)=\left\{i|-{\displaystyle \underset{j=1}{\overset{M}{\sum }}{\lambda }_j{c}_{i,j}\ge g_i\left(u_i\right)},i=1,\cdots ,N\right\}$. (11)

Let $m=\left|J\left({\lambda }_M\right)\right|$ , $n=\left|I\left({\lambda }_M\right)\right|$ , $n_L=\left|I_L\left({\lambda }_M\right)\right|$ , and $n_U=\left|I_U\left({\lambda }_M\right)\right|$. For the given ${\lambda }_M>0$ , without changing the orders of indices j and i, we re-index the constraints in the active constraint set $J\left({\lambda }_M\right)$ as $j=1,\cdots ,m$ , and we re-index the variables in the non-bound variable set $I\left({\lambda }_M\right)$ as $i=1,\cdots ,n$ , and re-index the variables in $I_L\left({\lambda }_M\right)$ and $I_U\left({\lambda }_M\right)$ as $i=1,\cdots ,n_L$ , and $i=1,\cdots ,n_U$ , respectively. As a result, constraint M in the original problem is re-indexed as constraint m, and ${\lambda }_M$ is also restated as ${\lambda }_m$.

We define $G_j\left(\lambda {}_1,\cdots ,{\lambda }_m\right)\equiv {\displaystyle \underset{i=1}{\overset{N}{\sum }}{c}_{i,j}{x}_i\left(\lambda \right)-C_j}=0$ , $j=1,\cdots ,m-1$ , and substitute

$x_i\left(\lambda \right)=\min \left\{\max \left\{h_i\left(-{\displaystyle \underset{j=1}{\overset{M}{\sum }}{\lambda }_j{c}_{i,j}}\right),l_i\right\},u_i\right\}$

into $G_j\left({\lambda }_1,\cdots ,{\lambda }_m\right)$ , then we have

$G_j\left({\lambda }_1,\cdots ,{\lambda }_m\right)\equiv {\displaystyle {\sum }_{i=1}^n{c}_{i,j}{h}_i\left(-{\displaystyle {\sum }_{s=1}^m{\lambda }_s{c}_{i,s}}\right)}-\left(C_j-{\displaystyle {\sum }_{i=1}^{n_L}{c}_{i,j}{l}_i-{\displaystyle {\sum }_{i=1}^{n_U}{c}_{i,j}{u}_i}}\right)=0$ , (12)

Taking the derivative of Equation (12), we get

$\begin{array}{c}\frac{\text{d}{G}_j\left({\lambda }_1,\cdots ,{\lambda }_m\right)}{\text{d}{\lambda }_m}=-\left[{\displaystyle \underset{i=1}{\overset{n}{\sum }}\frac{c_{i,j}}{k_i\left(x_i\left({\lambda }_1,\cdots ,{\lambda }_m\right)\right)}{\displaystyle \underset{s=1}{\overset{m}{\sum }}\frac{\text{d}{\lambda }_s}{\text{d}{\lambda }_m}{c}_{i,s}}}\right]\\ =-{\displaystyle \underset{s=1}{\overset{m}{\sum }}{\displaystyle \underset{i=1}{\overset{n}{\sum }}\frac{c_{i,j}{c}_{i,s}}{k_i\left(x_i\left({\lambda }_1,\cdots ,{\lambda }_m\right)\right)}\frac{\text{d}{\lambda }_s}{\text{d}{\lambda }_m}}}\\ =0,j=1,\cdots ,m-1\end{array}$ , (13)

where $k_i\left(x_i\right)=\text{d}{g}_i\left(x_i\right)/\text{d}{x}_i$.

Since $f_i\left(x_i\right)$ , $i=1,\cdots ,n$ are differentiable convex, we know $g_i\left(x_i\right)$ is increasing and $k_i\left(x_i\left({\lambda }_1,\cdots ,{\lambda }_m\right)\right)>0$. Note that ${\hat{f}}_i\left({\lambda }_M\right)=f_i+{\lambda }_M{c}_{i,M}{x}_i$ has the same structure as $f_i\left(x_i\right)$. So we define

${\rho }_i=\frac{1}{k_i\left(x_i\left({\lambda }_1,\cdots ,{\lambda }_m\right)\right)}>0$ ,

$i=1,\cdots ,n$ , and $a_{js}={\displaystyle {\sum }_{i=1}^n{\rho }_i{c}_{i,j}{c}_{i,s}}$ , $j,s=1,\cdots ,m$ , then Equation (13) can be rewritten in matrix form:

$\left(\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1m}\\ a_{21}& a_{22}& \cdots & a_{2m}\\ ⋮& ⋮& \ddots & ⋮\\ a_{\left(m-1\right)1}& a_{\left(m-1\right)2}& \cdots & a_{\left(m-1\right)m}\end{array}\right)\left(\begin{array}{c}\text{d}{\lambda }_1/\text{d}{\lambda }_m\\ ⋮\\ \text{d}{\lambda }_{m-1}/\text{d}{\lambda }_m\\ 1\end{array}\right)=\left(\begin{array}{c}0\\ 0\\ ⋮\\ 0\end{array}\right).$ (14)

In order to solve $\frac{\text{d}{\lambda }_j}{\text{d}{\lambda }_m}$ , $j=1,\cdots ,m-1$ , from Equation (14), we further define

$H_m=\left|\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1m}\\ a_{21}& a_{22}& \cdots & a_{2m}\\ ⋮& ⋮& \ddots & ⋮\\ a_{m1}& a_{m2}& \cdots & a_{mm}\end{array}\right|$ , (15)

and denote by $H_{j\left(m-1\right)}$ , $j=1,\cdots ,m-1$ the m-1 dimensional determinant in which the j column of $H_{m-1}$ is replaced by ${\left(a_{1m},a_{2m},\cdots ,a_{\left(m-1\right)m}\right)}^{\text{T}}$. We have the following formula from Equation (14) and Equation (15):

$\frac{\text{d}{\lambda }_j}{\text{d}{\lambda }_m}=-\frac{H_{j\left(m-1\right)}}{H_{m-1}},j=1,\cdots ,m-1,m>1$. (16)

Notice that the above results have similar structures as the results in Zhang [27]. Using the similar way, we can prove that

$\begin{array}{l}\frac{\text{d}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{c}_{i,m}{x}_i\left({\lambda }_m\right)}}{\text{d}{\lambda }_m}=-{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\rho }_i{c}_{i,m}\left({\displaystyle \underset{j=1}{\overset{m-1}{\sum }}{c}_{i,j}\frac{\text{d}{\lambda }_j}{\text{d}{\lambda }_m}}+c_{i,m}\right)}\\ =-{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\rho }_i{c}_{i,m}\left(-{\displaystyle \underset{j=1}{\overset{m-1}{\sum }}{c}_{i,j}\frac{H_{j\left(m-1\right)}}{H_{m-1}}}+c_{i,m}\right)}=-\frac{H_m}{H_{m-1}}<0\end{array}$. (17)

Since constraint M in the original problem is re-indexed as constraint m, and ${\lambda }_M$ is also restated as ${\lambda }_m$ , then ${\displaystyle {\sum }_{i=1}^n{c}_{i,m}{x}_i\left({\lambda }_m\right)}$ is equivalent to ${\displaystyle {\sum }_{i=1}^N{c}_{i,M}{x}_i\left({\lambda }_M\right)}$ in problem P with the original index, thus we know that ${\displaystyle {\sum }_{i=1}^N{c}_{i,M}{x}_i\left({\lambda }_M\right)}$ is a decreasing in ${\lambda }_M$.

Therefore, there are three possible cases: 1) When ${\lambda }_M=0$ , we get the optimal solution to problem $P\left(f,M\right)$ by solving problem $P\left(f,M-1\right)$ ; 2) If ${\lambda }_M>0$ and $m=1$ , we obtain the optimal solution to problem $P\left(f,M\right)$ by setting $x_i\left({\lambda }_M\right)=\min \left\{\max \left\{h_i\left(-{\lambda }_M{c}_{i,M}\right),l_i\right\},u_i\right\}$ ; 3) When ${\lambda }_M>0$ and $m>1$ , we can solve problem $P\left(f,M\right)$ by studying problem $P\left(\hat{f}\left({\lambda }_M\right),M-1\right)$ , with ${\hat{f}}_i\left({\lambda }_M\right)=f_i+{\lambda }_M{c}_{i,M}{x}_i$.

3.2. Solution Method

According to Proposition 2, we can solve $x^{*}$ by searching the optimal value of $\lambda$. Before presenting the solution method, we first study the bounds for $\lambda$. The lower bound for $\lambda$ is 0, and the upper bound for $\lambda$ is given in the following proposition.

Proposition 5. The upper bound of ${\lambda }_i$ is $\max \left(0,{\max }_{i=1,\cdots ,N}\left\{-g_i\left(l_i\right)/c_{i,M}\right\}\right)$.

From Proposition 4, we get the optimal value of $x^{*}$ if the optimal solution $x\left({\lambda }_M\right)$ to problem $P\left(\hat{f}\left({\lambda }_M\right),M-1\right)$ satisfies

${\lambda }_M\left({\displaystyle \underset{i=1}{\overset{N}{\sum }}{c}_{i,M}{x}_i\left({\lambda }_M\right)-C_M}\right)=0$.

Since ${\displaystyle {\sum }_{i=1}^N{c}_{i,M}{x}_i\left({\lambda }_M\right)}$ is decreasing in ${\lambda }_M$ , the optimal solution can be found by applying the binary search over $\left[0,\max \left(0,{\max }_{i=1,\cdots ,N}\left\{-g_i\left(l_i\right)/c_{i,M}\right\}\right)\right]$. Since Problem $P\left(\hat{f}\left({\lambda }_M\right),M-1\right)$ has the same structure as problem $P\left(f,M\right)$ , we can use a multi-tier binary search method to solve problem P. Main steps of the multi-tier binary search method are given in Algorithm 1.

Algorithm 1: $\text{SloveP}\left(f,M\right)$

Step 1: If $M=0$ , then let $x_i^{*}=\min \left\{\max \left\{\begin{array}{c}\arg \end{array}\left\{g_i\left(x_i\right)=0\right\},l_i\right\},u_i\right\}$ , stop;

Step 2: Let ${\lambda }_M^L=0$ , ${\lambda }_M^U=\max \left(0,{\max }_{i=1,\cdots ,N}\left\{-g_i\left(l_i\right)/c_{i,M}\right\}\right)$ ;

Step 3: Let ${\lambda }_M=\left({\lambda }_M^L+{\lambda }_M^U\right)/2$ ;

Step 4: If ${\lambda }_M=0$ , then let $x_i^{*}=\text{SolveP}\left(f,M-1\right)$ and ${\lambda }_M^{*}=0$ , stop;

Step 5: If $M=1$ , then let $x_i\left({\lambda }_M\right)=\min \left\{\max \left\{h_i\left(-{\lambda }_M{c}_{i,M}\right),l_i\right\},u_i\right\}$ ;

If $M>1$ , then let $x_i\left({\lambda }_M\right)=\text{SolveP}\left(\hat{f}\left({\lambda }_M\right),M-1\right)$

Step 6: If ${\displaystyle {\sum }_{i=1}^N{c}_{i,M}{x}_i\left({\lambda }_M\right)}>C_M$ , then let ${\lambda }_M^L={\lambda }_M$ , go to Step 3;

If ${\displaystyle {\sum }_{i=1}^N{c}_{i,M}{x}_i\left({\lambda }_M\right)}<C_M$ , then let ${\lambda }_M^U={\lambda }_M$ , go to Step 3;

Step 7: Let $x_i^{*}=x_i\left({\lambda }_M\right)$ and ${\lambda }_M^{*}={\lambda }_M$ , stop.

In the algorithm, we first solve the unconstrained problem with bounded variables (Step 1) to obtain $x^{*}$. If the constraints are active, we apply the binary search procedure (Step 2 - 7) over interval $\left[{\lambda }_M^L,{\lambda }_M^U\right]$ to determine ${\lambda }_M^{*}$. If either ${\lambda }_M=0$ or ${\displaystyle {\sum }_{i=1}^N{c}_{i,M}{x}_i\left({\lambda }_M\right)}=C_M$ , the binary search procedure terminates. If ${\displaystyle {\sum }_{i=1}^N{c}_{i,M}{x}_i\left({\lambda }_M\right)}\le C_M$ is not binding, then the iterating process will end in Step 4 with ${\lambda }_M=0$. Therefore, we can get the optimal solution $x_i^{*}$ by solving problem $P\left(f,M-1\right)$. If the constraint ${\displaystyle {\sum }_{i=1}^N{c}_{i,M}{x}_i\left({\lambda }_M\right)}\le C_M$ is active, the so-

lution procedure will stop at Step 7 with ${\displaystyle {\sum }_{i=1}^N{c}_{i,M}{x}_i\left({\lambda }_M\right)}=C_M$. Step 5 derives

$x_i\left({\lambda }_M\right)$ by solving problem $P\left(\hat{f}\left({\lambda }_M\right),M-1\right)$ with ${\hat{f}}_i\left({\lambda }_M\right)=f_i+{\lambda }_M{c}_{i,M}{x}_i$ for the given ${\lambda }_M>0$. If $M=1$ , problem $P\left(\hat{f}\left({\lambda }_M\right),M-1\right)$ has no knapsack

constraint, and hence we have $x_i\left({\lambda }_M\right)=\min \left\{\max \left\{h_i\left(-{\lambda }_M{c}_{i,M}\right),l_i\right\},u_i\right\}$. If

$M>1$ , we can solve the problem recursively. Problem $P\left(\hat{f}\left({\lambda }_M\right),M-1\right)$ has the same structure as problem $P\left(f,M\right)$ , and hence the algorithm can call itself recursively to solve the problem $P\left(\hat{f}\left({\lambda }_M\right),M-1\right)$.

The algorithm is a recursive algorithm with M tiers of binary search loop. The computational complexity of M-tier binary search procedure is $O\left({\left({\log }_2\left(1/\epsilon \right)\right)}^M\right)$ , where $\epsilon$ is the error target for the binary search. The computational complexity of the last recursive step is $O\left(N\right)$. Therefore, the proposed algorithm has the computational complexity $O\left({\left({\log }_2\left(1/\epsilon \right)\right)}^{M}N\right)$ , which is polynomial in the number of decision variables N.

4. Numerical Study

The solution method developed in this paper can be used for solving the continuous nonlinear multidimensional knapsack problems with general structure, so many application problems with different objective functions summarized in Zhang and Hua with multiple constraints can be used to show the application of our method [16].

In our numerical study, we first show the application of our method using two examples: quadratic multidimensional knapsack problem (QMK) and the production planning problem presented in Bretthauer and Shetty [12]. Then we use the statistical study to show the efficiency of our method. All computational experiments are conducted on a laptop (dual processor 2.00 GHz, memory 2.96G) with Matlab R2011a.

4.1. The Illustrative Examples

The first illustrative example is a separable quadratic knapsack problem. We set the objective function as $f_i\left(x_i\right)=a_i{\left(x_i-b_i\right)}^2$ , $a_i>0,i=1,\cdots ,N$. It has two resource constraints: C₁ = 12,000 and C₂ = 10,000. Table 2 gives the relevant information for this example. $x_i^{*}$ is the optimal solution obtained by applying our algorithm. To show the efficiency of our method, we plot the values of ${\lambda }_M^L,{\lambda }_M^U$ and ${\lambda }_M$ in the iteration process for solving the example in Figure 1. Figure 1 shows our algorithm can solve the problem within very limited iterations.

In the second example, we solve the production planning problem in Bretthauer and Shetty [12]. The objective function was set as

$f_i\left(x_i\right)=\min {\displaystyle \underset{i=1}{\overset{n}{\sum }}\left(h_i+d_i{x}_i+\frac{e_i}{x_i}\right)}$ ,

$i=1,\cdots ,N$. There are three resource constraints: C₁ = 200, C₂ = 300, and C₃ = 500. We use the same parameters used in Bretthauer and Shetty [12]. The relevant

Table 2. Parameters and solution for the first example.

Figure 1. ${\lambda }_M^L,{\lambda }_M^U,{\lambda }_M$ in the iteration process for solving the first example.

information for this example is listed in Table 3. $x_i^{*}$ is the optimal solution obtained by applying our algorithm.

4.2. The Statistical Results

In this subsection, we present two numerical experiments to show the effectiveness of our method for solving problems with different scale and objective functions. In the first experiment, parameters of the QMK problems are all randomly generated. We use the notation $z~U\left(\alpha ,\beta \right)$ to denote that z is uniformly generated over $\left[\alpha ,\beta \right]$. The parameters of QMK instances are generated as follows: $a_i~U\left(1,2\right)$ , $b_i~U\left(5,10\right)$ , $c_{i,j}~U\left(1,10\right)$ , $l_i~U\left(5,15\right)$ , $u_i~U\left(20,30\right)$ and $C_j~N\times U\left(\text{100000},\text{200000}\right)$ , for $i=1,\cdots ,N;j=1,\cdots ,M$.

In this experiment, we set problems with different sizes, respectively with M = 4 and N = 10, M = 2 and N = 100, M = 3 and N = 100, M = 2 and N = 1000. For each problem size, 50 test instances are randomly generated. The statistical results on number of iterations and computation time (in seconds) are reported in Table 4.

In the Second experiment, we solve the production planning problem with randomly generated parameters. The parameters of the instances are generated as follows: $d_i~U\left(30,50\right)$ , $e_i~U\left(100,200\right)$ , $c_{i,j}~U\left(10,50\right)$ , $l_i~U\left(1,5\right)$ , $u_i~U\left(20,30\right)$ and $C_j~N\times U\left(\text{100000},\text{200000}\right)$ , for $i=1,\cdots ,N;j=1,\cdots ,M$.

In this experiment, we set problems with different sizes, respectively with M = 4 and N = 10, M = 2 and N = 100, M = 3 and N = 100, M = 2 and N = 1000. For each problem size, we randomly generated 50 test instances. The statistical results on number of iterations and computation time (in seconds) are presented in Table 5.

From Table 4 and Table 5, we observe that the standard deviations of number of iterations and computation times are quite low. It implies that our method is quite effective with different objective functions. We also observe that the

Table 3. Parameters and solutions for the second example.

Table 4. Statistical results for randomly generated QMK problems.

Table 5. Statistical results for randomly generated production planning problems.

computation time is more sensitive to the number of the resource constraints rather than the number of variables. Since the application problems often have much more variables than knapsack constraints, our algorithm is useful in practice.

5. Conclusions

In this paper, we study a class of continuous separable nonlinear multidimensional knapsack problems. By analyzing the structural properties of the optimal solution, we develop a multi-tier binary solution method. The proposed method has following advantages. 1) It is applicable for solving the nonlinear multidimensional knapsack problems with general structure. 2) It has computational complexity of polynomial in the number of variables.

This research can be further extended in several ways. One is to study non-separable multidimensional knapsack problems using the similar idea. Another way is to develop exact solution methods or heuristics for solving the integer multidimensional knapsack problems based on our method. Finally, the idea used in this study can be extended for investigating other complex optimization problems with multiple constraints.

Acknowledgements

This work is supported by National Natural Science Foundation of China (Grants No. 71672199).

Appendix

A.1 Proof of Proposition 2

It is defined that $0\le l_i<u_i$ for all $i=1,\cdots ,N$. The optimal solution to problem PR should satisfy Equation (1) and Equation (3). If $l_i\le {\bar{x}}_i\le u_i$ , it means that the bound constraint is inactive. Therefore, we have ${\hat{x}}_i={\bar{x}}_i$. Since $g_i\left(x_i\right)$ is increasing in $x_i$ , and $g_i\left({\bar{x}}_i\right)=0$ , we have $g_i\left(x_i\right)\ge 0$ if $x_i>{\bar{x}}_i$. If ${\bar{x}}_i<l_i<u_i$ , then $g_i\left(x_i\right)\ge 0$ if $l_i\le x_i\le u_i$. Thus for any $x_i\in \left[l_i,u_i\right]$ , we have $f_i\left(l_i\right)\le f_i\left(x_i\right)$ , and ${\hat{x}}_i=l_i$. If $l_i<u_i<{\bar{x}}_i$ , we have ${\hat{x}}_i=u_i$. It can be proved similar to the condition of ${\bar{x}}_i<l_i$.

A.2 Proof of Proposition 3

1) If $g_i\left(l_i\right)\le -{\displaystyle {\sum }_{j=1}^M{\lambda }_j}{c}_{i,j}\le g_i\left(u_i\right)$ , then we have $l_i\le h_i\left(-{\displaystyle {\sum }_{j=1}^M{\lambda }_j}{c}_{i,j}\right)\le u_i$ , and $w_i=v_i=0$ , which implies $x_i\left(\lambda \right)=h_i\left(-{\displaystyle {\sum }_{j=1}^M{\lambda }_j}{c}_{i,j}\right)$. If $-{\displaystyle {\sum }_{j=1}^M{\lambda }_j}{c}_{i,j}<g_i\left(l_i\right)$ , then we have $g_i\left(x_i\right)+{\displaystyle {\sum }_{j=1}^M{\lambda }_j}{c}_{i,j}\ge g_i\left(l_i\right)+{\displaystyle {\sum }_{j=1}^M{\lambda }_j}{c}_{i,j}>0$ , and hence $w_i>0$ , $x_i\left(\lambda \right)=l_i$.

If $-{\displaystyle {\sum }_{j=1}^M{\lambda }_j}{c}_{i,j}>g_i\left(u_i\right)$ , we have $g_i\left(x_i\right)+{\displaystyle {\sum }_{j=1}^M{\lambda }_j}{c}_{i,j}\le g_i\left(u_i\right)+{\displaystyle {\sum }_{j=1}^M{\lambda }_j}{c}_{i,j}<0$ , which means $v_i>0$ , and $x_i\left(\lambda \right)=u_i$. Therefore, we have

$x_i\left(\lambda \right)=\{\begin{array}{l}{l}_i,\text{if}-{\displaystyle \underset{j=1}{\overset{M}{\sum }}{\lambda }_j{c}_{i,j}}<g_i\left(l_i\right),\\ h_i\left(-{\displaystyle \underset{j=1}{\overset{M}{\sum }}{\lambda }_j{c}_{i,j}}\right),\text{if}{g}_i\left(l_i\right)\le -{\displaystyle \underset{j=1}{\overset{M}{\sum }}{\lambda }_j{c}_{i,j}}\le g_i\left(u_i\right),\\ u_i,\text{if}-{\displaystyle \underset{j=1}{\overset{M}{\sum }}{\lambda }_j{c}_{i,j}}>g_i\left(u_i\right).\end{array}$ (A1)

2) ${\lambda }_j=0$ or ${\displaystyle {\sum }_{i=1}^N{c}_{i,j}{x}_i}\left(\lambda \right)=C_j$ implies

${\lambda }_j\left({\displaystyle \underset{i=1}{\overset{N}{\sum }}{c}_{i,j}{x}_i\left(\lambda \right)-C_j}\right)=0,j=1,\cdots ,M$.

Because $x\left(\lambda \right)$ satisfies Equation (5) and Equation (6), $x\left(\lambda \right)$ will satisfy all KKT conditions. Therefore, $x^{*}=x\left(\lambda \right)$ if $\left(x\left(\lambda \right),\lambda \right)$ satisfies ${\lambda }_j=0$ or ${\displaystyle {\sum }_{i=1}^N{c}_{i,j}{x}_i}=C_j$ , $j=1,\cdots ,M$.

A.3 Proof of Proposition 4

1) ${\lambda }_M=0$ or ${\displaystyle {\sum }_{i=1}^N{c}_{i,M}{x}_i\left({\lambda }_M\right)}=C_M$ implies ${\lambda }_M\left({\displaystyle {\sum }_{i=1}^N{c}_{i,M}{x}_i-C_M}\right)=0$. Since $\left(x\left(\lambda {}_M\right),{\lambda }_M\right)$ satisfies Equation (5) and Equation (6), it will satisfy all KKT conditions. Therefore, $x^{*}=x\left({\lambda }_M\right)$ if $\left(x\left(\lambda {}_M\right),{\lambda }_M\right)$ satisfies ${\lambda }_M=0$ or ${\displaystyle {\sum }_{i=1}^N{c}_{i,M}{x}_i\left({\lambda }_M\right)}=C_M$.

2) KKT conditions for problem $P\left(\hat{f}\left({\lambda }_M\right),M-1\right)$ are

$\frac{\text{d}{\hat{f}}_i\left(x_i\right)}{\text{d}{x}_i}+{\displaystyle \underset{j=1}{\overset{M-1}{\sum }}{\lambda }_j{c}_{i,j}-w_i+v_i}=0,i=1,\cdots ,N,$ (A2)

${\displaystyle \underset{i=1}{\overset{N}{\sum }}{w}_i\left(x_i-l_i\right)+{\displaystyle \underset{i=1}{\overset{N}{\sum }}{v}_i\left(x_i-u_i\right)}=0},$ (A3)

${\lambda }_j\left({\displaystyle \underset{i=1}{\overset{N}{\sum }}{c}_{i,j}{x}_i-C_j}\right)=0,j=1,\cdots ,M-1$. (A4)

Notice that ${\hat{f}}_i\left(x_i\right)$ is a parameter-adjusted function of $f_i\left(x_i\right)$ , with ${\hat{f}}_i\left({\lambda }_M\right)=f_i+{\lambda }_M{c}_{i,M}{x}_i$. These conditions in Equations (A2)-(A3) are the same as KKT conditions given in Equations (5)-(7) without ${\lambda }_M\left({\displaystyle {\sum }_{i=1}^N{c}_{i,M}{x}_i-C_M}\right)=0$. Since $x\left({\lambda }_M\right)$ is the optimal solution of the KKT conditions in Equations (5)-(7) without ${\lambda }_M\left({\displaystyle {\sum }_{i=1}^N{c}_{i,M}{x}_i-C_M}\right)=0$ , it must be the optimal solution to problem $P\left(\hat{f}\left({\lambda }_M\right),M-1\right)$.

A.4 Proof of Proposition 5

Let ${\bar{\lambda }}_M=\max \left(0,{\max }_{i=1,\cdots ,N}\left\{-g_i\left(l_i\right)/c_{i,M}\right\}\right)$. If ${\lambda }_M^{*}>{\bar{\lambda }}_M$ , then we have ${\lambda }_M^{*}{c}_{i,M}>-g_i\left(l_i\right),i=1,\cdots ,N$. From Equation (5), we have

$w_i^{*}=g_i\left(x_i\right)+{\displaystyle \underset{j=1}{\overset{M}{\sum }}{\lambda }_j{c}_{i,j}+v_i}>g_i\left(x_i\right)+{\lambda }_M^{*}{c}_{i,M}+v_i>0,i=1,\cdots ,N$. (A5)

Since $w_i>0$ , from Equation (6), we know $x_i=l_i$ and $v_i=0$. Thus, we have

${\lambda }_M^{*}\left({\displaystyle \underset{i=1}{\overset{N}{\sum }}{c}_{i,M}{x}_i^{*}}-C_M\right)={\lambda }_M^{*}\left({\displaystyle \underset{i=1}{\overset{N}{\sum }}{c}_{i,M}{l}_i}-C_M\right)\ne 0$. (A6)

Equation (A6) violates the slackness condition ${\lambda }_M\left({\displaystyle \underset{i=1}{\overset{N}{\sum }}{c}_{i,M}{x}_i-C_M}\right)=0$ in Equation (7). Therefore, there must be ${\lambda }_M^{*}<{\bar{\lambda }}_M$.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Bonnans, J.F. (1994) Local Analysis of Newton-Type Methods for Variational Inequalities and Nonlinear Programming. Applied Mathematics and Optimization, 29, 161-186. https://doi.org/10.1007/BF01204181
[2]	Coleman, T.F. and Li, Y. (1994) On the Convergence of Interior-Reflective Newton Methods for Nonlinear Minimization Subject to Bounds. Mathematical Programming, 67, 189-224. https://doi.org/10.1007/BF01582221
[3]	Van Hentenryck, P., Michel, L. and Benhamou, F. (1998) Constraint Programming over Nonlinear Constraints. Science of Computer Programming, 30, 83-118. https://doi.org/10.1016/S0167-6423(97)00008-7
[4]	Borchers, B. and Mitchell, J.E. (1994) An Improved Branch and Bound Algorithm for Mixed Integer Nonlinear Programs. Computers & Operations Research, 21, 359-367. https://doi.org/10.1016/0305-0548(94)90024-8
[5]	Leyffer, S. (2001) Integrating SQP and Branch-and-Bound for Mixed Integer Nonlinear Programming. Computational Optimization and Applications, 18, 295-309. https://doi.org/10.1023/A:1011241421041
[6]	Byrd, R.H., Hribar, M.E. and Nocedal, J. (1999) An Interior Point Algorithm for Large-Scale Nonlinear Programming. SIAM Journal on Optimization, 9, 877-900. https://doi.org/10.1137/S1052623497325107
[7]	Benson, H.Y., Vanderbei, R.J. and Shanno, D.F. (2002) Interior-Point Methods for Nonconvex Nonlinear Programming: Filter Methods and Merit Functions. Computational Optimization and Applications, 23, 257-272. https://doi.org/10.1023/A:1020533003783
[8]	Spellucci, P. (1998) An SQP Method for General Nonlinear Programs Using only Equality Constrained Subproblems. Mathematical Programming, 82, 413-448. https://doi.org/10.1007/BF01580078
[9]	Zhu, Z. and Zhang, K. (2004) A New SQP Method of Feasible Directions for Nonlinear Programming. Applied Mathematics and Computation, 148, 121-134. https://doi.org/10.1016/S0096-3003(02)00832-9
[10]	Nie, P.Y. and Ma, C.F. (2006) A Trust Region Filter Method for General Non-Linear Programming. Applied Mathematics and Computation, 172, 1000-1017. https://doi.org/10.1016/j.amc.2005.03.004
[11]	Nie, P.Y. (2007) Sequential Penalty Quadratic Programming Filter Methods for Nonlinear Programming. Nonlinear Analysis: Real World Applications, 8, 118-129. https://doi.org/10.1016/j.nonrwa.2005.06.003
[12]	Bretthauer, K.M. and Shetty, B. (2002b) The Nonlinear Knapsack Problem-Algo- rithms and Applications. European Journal of Operational Research, 138, 459-472. https://doi.org/10.1016/S0377-2217(01)00179-5
[13]	Bretthauer, K.M. and Shetty, B. (1995) The Nonlinear Resource Allocation Problem. Operations Research, 43, 670-683. https://doi.org/10.1287/opre.43.4.670
[14]	Kodialam, M.S. and Luss, H. (1998) Algorithms for Separable Nonlinear Resource Allocation Problems. Operations Research, 46, 272-284. https://doi.org/10.1287/opre.46.2.272
[15]	Bretthauer, K.M. and Shetty, B. (2002a) A Pegging Algorithm for the Nonlinear Resource Allocation Problem. Computers & Operations Research, 29, 505-527. https://doi.org/10.1016/S0305-0548(00)00089-7
[16]	Zhang, B. and Hua, Z. (2008) A Unified Method for a Class of Convex Separable Nonlinear Knapsack Problems. European Journal of Operational Research, 191, 1-6. https://doi.org/10.1016/j.ejor.2007.07.005
[17]	Kiwiel, K.C. (2008) Breakpoint Searching Algorithms for the Continuous Quadratic Knapsack Problem. Mathematical Programming, 112, 473-491. https://doi.org/10.1007/s10107-006-0050-z
[18]	Sharkey, T.C., Romeijn, H.E. and Geunes, J. (2011) A Class of Nonlinear Nonseparable Continuous Knapsack and Multiple-Choice Knapsack Problems. Mathematical Programming, 126, 69-96. https://doi.org/10.1007/s10107-009-0274-9
[19]	Morin, T.L. and Marsten, R.E. (1976) An Algorithm for Nonlinear Knapsack Problems. Management Science, 22, 1147-1158. https://doi.org/10.1287/mnsc.22.10.1147
[20]	Ohtagaki, H., Iwasaki, A., Nakagawa, Y. and Narihisa, H. (2000) Smart Greedy Procedure for Solving a Multidimensional Nonlinear Knapsack Class of Reliability Optimization Problems. Mathematical and Computer Modelling, 31, 283-288. https://doi.org/10.1016/S0895-7177(00)00097-2
[21]	Li, D., Sun, X.L. and Wang, F.L. (2006) Convergent Lagrangian and Contour Cut Method for Nonlinear Integer Programming with a Quadratic Objective Function. SIAM Journal on Optimization, 17, 372-400. https://doi.org/10.1137/040606193
[22]	Li, D., Sun, X.L., Wang, J. and McKinnon, K.I. (2009) Convergent Lagrangian and Domain Cut Method for Nonlinear Knapsack Problems. Computational Optimization and Applications, 42, 67-104. https://doi.org/10.1007/s10589-007-9113-1
[23]	Quadri, D., Soutif, E. and Tolla, P. (2009) Exact Solution Method to Solve Large Scale Integer Quadratic Multidimensional Knapsack Problems. Journal of Combinatorial Optimization, 17, 157-167. https://doi.org/10.1007/s10878-007-9105-1
[24]	Wang, H., Kochenberger, G. and Glover, F. (2012) A Computational Study on the Quadratic Knapsack Problem with Multiple Constraints. Computers & Operations Research, 39, 3-11. https://doi.org/10.1016/j.cor.2010.12.017
[25]	Abdel-Malek, L.L. and Areeratchakul, N. (2007) A Quadratic Programming Approach to the Multi-Product Newsvendor Problem with Side Constraints. European Journal of Operational Research, 176, 1607-1619. https://doi.org/10.1016/j.ejor.2005.11.002
[26]	Abdel-Malek, L.L. and Otegbeye, M. (2013) Separable Programming/Duality Approach to Solving the Multi Constrained Newsboy/Gardener Problem. Applied Mathematical Modelling, 37, 4497-4508. https://doi.org/10.1016/j.apm.2012.09.059
[27]	Zhang, B. (2012) Multi-Tier Binary Solution Method for Multi-Product Newsvendor Problem with Multiple Constraints. European Journal of Operational Research, 218, 426-434. https://doi.org/10.1016/j.ejor.2011.10.053