# Knapsack problem formulation

I would like to use a knapsack problem formulation based on Lucas paper. Namely, I try to implement the following math formula for Hamiltonian $$H =H_A + H_B$$, where $$H_{B} = -B\sum_{i=1}^{n}v_{i}x_{i}; H_{A} = A(1-\sum_{j=0}^{M-1}2^{j}y_{j}+(c+1-2^{M})y_{M}-\sum_{i=1}^{n}w_{i}x_{i})^{2}\,.$$

I implemented it using sympy (imported as sp):


def knapsack(values, weights, capacity, A, B):

n = len(values)
M = math.ceil(math.log2(capacity))

# Define the variables
x = [sp.Symbol(f'x{i}') for i in range(n)]
y = [sp.Symbol(f'y{i}') for i in range(M+1)]

# Define the objective function
HB = -B * sum(values[i] * x[i] for i in range(n))
HA = A * (1 - sum(2**j * y[j] for j in range(M)) + (capacity+1-2**M)*y[M] - sum(weights[i] * x[i] for i in range(n)))**2
H = HA + HB

# Expand the objective function
expanded_H = H.expand()

# Extract linear, quadratic, and constant terms from the expanded objective function
linear_terms = {str(symbol): coeff for symbol, coeff in expanded_H.as_coefficients_dict().items() if symbol.is_symbol}

for term, coeff in expanded_H.as_coefficients_dict().items():
if isinstance(term, sp.Mul) and len(term.args) == 2:
v1, v2 = sorted([str(a) for a in term.args])
if (v1, v2) not in quadratic_terms:
else:

constant = expanded_H.as_coefficients_dict().get(1, 0)

for var in x + y:
qp.binary_var(name=str(var))

# Set the objective in QuadraticProgram

Next, I convert it to the Ising model to get an operator to pass to VQE. Next, I set up VQE to compute the minimum eigenvalue. My question is if the above formulation of the knapsack is correct. When running VQE (on statevector_simulator), I get solutions that are more or less random, sometimes suboptimal, sometimes incorrect (breaking capacity constraint). I suppose I missed something in the formulation of the problem.