diff --git a/docs/modules/models/knapsack.md b/docs/modules/models/knapsack.md
index 581829b9..f6802dda 100644
--- a/docs/modules/models/knapsack.md
+++ b/docs/modules/models/knapsack.md
@@ -1,5 +1,11 @@
 # Knapsack problem
 
+$$\text{Minimize } - \sum_{i = 1}^{N} c_{i}x_{i} + A_{\text{totalWeight}} \left ( \sum_{i = 1}^{N} w_{i}x_{i} - \sum_{k = 1}^{W} ky_{k} \right )^{2} + A_{\text{onlyOneWeight}} \left ( 1 - \sum_{k = 0}^{W} y_{k} \right )^{2}$$
+
+with:
+- $A_{\text{totalWeight}} = sum_{i = 1}^{N} c_{i} + 1$
+- $A_{\text{onlyOneWeight}} = \left ( sum_{i = 1}^{N} c_{i} + 1 \right ) \left ( W + sum_{i = 1}^{N} w_{i} + 1 \right)$
+
 ```{eval-rst}
 .. autoclass:: simulated_bifurcation.models.Knapsack
     :members:
diff --git a/src/simulated_bifurcation/models/__init__.py b/src/simulated_bifurcation/models/__init__.py
index 18aea333..0d611326 100644
--- a/src/simulated_bifurcation/models/__init__.py
+++ b/src/simulated_bifurcation/models/__init__.py
@@ -3,14 +3,12 @@
 solved using the Simulated Bifurcation algorithm. This module is a
 pre-defined models library to extend the use of SB.
 
-It also comes with an abstract class called `ABCModel` which acts as a
+It also comes with an abstract class called `ABCQuadraticModel` which acts as a
 basis for defining custom optimization models that are not implemented
 yet.
 
 """
 
-from .abc_model import ABCModel
-from .ising import Ising
-from .markowitz import Markowitz, SequentialMarkowitz
-from .number_partitioning import NumberPartitioning
-from .qubo import QUBO
+from .knapsack import knapsack
+from .markowitz import markowitz, sequential_markowitz
+from .number_partitioning import number_partitioning
diff --git a/src/simulated_bifurcation/models/abc_model.py b/src/simulated_bifurcation/models/abc_model.py
deleted file mode 100644
index 253bd9a9..00000000
--- a/src/simulated_bifurcation/models/abc_model.py
+++ /dev/null
@@ -1,78 +0,0 @@
-from abc import ABC
-from typing import List, Literal, Optional, Tuple, Union
-
-import torch
-
-from ..core import QuadraticPolynomial
-
-
-class ABCModel(ABC, QuadraticPolynomial):
-    """
-    Abstract class that serves as a base component to define quadratic
-    optimization models that can be solved using the Simulated
-    Bifurcation algorithm.
-
-    Attributes
-    ----------
-    domain : str
-        The optimization domain of the problem.
-
-    """
-
-    domain: Union[str, List[str]]
-
-    def minimize(
-        self,
-        *,
-        agents: int = 128,
-        max_steps: int = 10000,
-        best_only: bool = True,
-        mode: Literal["ballistic", "discrete"] = "ballistic",
-        heated: bool = False,
-        verbose: bool = True,
-        early_stopping: bool = True,
-        sampling_period: int = 50,
-        convergence_threshold: int = 50,
-        timeout: Optional[float] = None,
-    ) -> Tuple[torch.Tensor, torch.Tensor]:
-        return super().minimize(
-            domain=self.domain,
-            agents=agents,
-            max_steps=max_steps,
-            best_only=best_only,
-            mode=mode,
-            heated=heated,
-            verbose=verbose,
-            early_stopping=early_stopping,
-            sampling_period=sampling_period,
-            convergence_threshold=convergence_threshold,
-            timeout=timeout,
-        )
-
-    def maximize(
-        self,
-        *,
-        agents: int = 128,
-        max_steps: int = 10000,
-        best_only: bool = True,
-        mode: Literal["ballistic", "discrete"] = "ballistic",
-        heated: bool = False,
-        verbose: bool = True,
-        early_stopping: bool = True,
-        sampling_period: int = 50,
-        convergence_threshold: int = 50,
-        timeout: Optional[float] = None,
-    ) -> Tuple[torch.Tensor, torch.Tensor]:
-        return super().maximize(
-            domain=self.domain,
-            agents=agents,
-            max_steps=max_steps,
-            best_only=best_only,
-            mode=mode,
-            heated=heated,
-            verbose=verbose,
-            early_stopping=early_stopping,
-            sampling_period=sampling_period,
-            convergence_threshold=convergence_threshold,
-            timeout=timeout,
-        )
diff --git a/src/simulated_bifurcation/models/abc_quadratic_model.py b/src/simulated_bifurcation/models/abc_quadratic_model.py
new file mode 100644
index 00000000..5888d3b1
--- /dev/null
+++ b/src/simulated_bifurcation/models/abc_quadratic_model.py
@@ -0,0 +1,60 @@
+from abc import ABC, abstractmethod
+from typing import Any, Literal, Optional, Union
+
+import torch
+
+from ..core.quadratic_polynomial import QuadraticPolynomial
+
+
+class ABCQuadraticModel(ABC):
+    domain: str
+    sense: str
+
+    @abstractmethod
+    def _as_quadratic_polynomial(
+        self, dtype: Optional[torch.dtype], device: Optional[Union[str, torch.device]]
+    ) -> QuadraticPolynomial:
+        raise NotImplementedError() # pragma: no cover
+
+    @abstractmethod
+    def _from_optimized_tensor(
+        self, optimized_tensor: torch.Tensor, optimized_cost: torch.Tensor
+    ) -> Any:
+        raise NotImplementedError() # pragma: no cover
+
+    def solve(
+        self,
+        dtype: Optional[torch.dtype] = None,
+        device: Optional[Union[str, torch.device]] = None,
+        agents: int = 128,
+        max_steps: int = 10_000,
+        mode: Literal["ballistic", "discrete"] = "ballistic",
+        heated: bool = False,
+        verbose: bool = True,
+        early_stopping: bool = True,
+        sampling_period: int = 50,
+        convergence_threshold: int = 50,
+        timeout: Optional[float] = None,
+    ) -> Any:
+        quadratic_model = self._as_quadratic_polynomial(dtype=dtype, device=device)
+        kwargs = {
+            "domain": self.domain,
+            "agents": agents,
+            "max_steps": max_steps,
+            "mode": mode,
+            "heated": heated,
+            "best_only": True,
+            "verbose": verbose,
+            "early_stopping": early_stopping,
+            "sampling_period": sampling_period,
+            "convergence_threshold": convergence_threshold,
+            "timeout": timeout,
+        }
+        if self.sense == "minimize":
+            return self._from_optimized_tensor(*quadratic_model.minimize(**kwargs))
+        elif self.sense == "maximize":
+            return self._from_optimized_tensor(*quadratic_model.maximize(**kwargs))
+        else: # pragma: no cover
+            raise ValueError(
+                f'Unknown optimization sense {self.sense}. Expected "maximize" or "minimize".'
+            )
diff --git a/src/simulated_bifurcation/models/ising.py b/src/simulated_bifurcation/models/ising.py
deleted file mode 100644
index 62399a34..00000000
--- a/src/simulated_bifurcation/models/ising.py
+++ /dev/null
@@ -1,30 +0,0 @@
-from typing import Optional, Union
-
-import numpy as np
-import torch
-
-from .abc_model import ABCModel
-
-
-class Ising(ABCModel):
-    """
-    Implementation of the Ising model.
-
-    Solving an Ising problem means searching the spin vector S (with values in
-    {-1, 1}) such that, given a matrix J with zero diagonal and a
-    vector h, the following quantity - called Ising energy - is minimal (S is
-    then called the ground state): `-0.5 * ΣΣ J(i,j)s(i)s(j) + Σ h(i)s(i)`
-    """
-
-    domain = "spin"
-
-    def __init__(
-        self,
-        J: Union[torch.Tensor, np.ndarray],
-        h: Union[torch.Tensor, np.ndarray, None] = None,
-        dtype: Optional[torch.dtype] = None,
-        device: Optional[Union[str, torch.device]] = None,
-    ) -> None:
-        super().__init__(-0.5 * J, h, dtype=dtype, device=device)
-        self.J = J
-        self.h = h
diff --git a/src/simulated_bifurcation/models/knapsack.py b/src/simulated_bifurcation/models/knapsack.py
new file mode 100644
index 00000000..2a6a4fc1
--- /dev/null
+++ b/src/simulated_bifurcation/models/knapsack.py
@@ -0,0 +1,88 @@
+from typing import List, Optional, Union
+
+import numpy as np
+import torch
+
+from ..core.quadratic_polynomial import QuadraticPolynomial
+from .abc_quadratic_model import ABCQuadraticModel
+
+
+class Knapsack(ABCQuadraticModel):
+    """
+    Note
+    ----
+    It is advised to run `solve` in discrete mode with `torch.float32` dtype for an optimal behavior.
+    """
+
+    domain = "binary"
+    sense = "minimize"
+
+    def __init__(
+        self, costs: List[Union[int, float]], weights: List[int], max_weight: int
+    ):
+        self.__costs = list(map(float, costs))
+        self.__weights = list(map(float, weights))
+        self.__max_weight = max_weight
+
+    def _as_quadratic_polynomial(
+        self, dtype: Optional[torch.dtype], device: Optional[Union[str, torch.device]]
+    ) -> QuadraticPolynomial:
+        # Define penalties for quadratic constraints
+        weight_constraint_penalty = float(np.sum(self.__costs) + 1)
+        only_one_weight_constraint_penalty = weight_constraint_penalty * (
+            self.__max_weight + 1 + float(np.sum(self.__weights))
+        )
+
+        # Define arrays and tensors
+        costs_array = np.array(self.__costs).reshape(-1, 1)
+        weights_array = np.array(self.__weights).reshape(-1, 1)
+        possible_total_weights_array = np.arange(self.__max_weight + 1).reshape(-1, 1)
+
+        quadratic = np.block(
+            [
+                [
+                    weight_constraint_penalty * weights_array @ weights_array.T,
+                    -weight_constraint_penalty
+                    * weights_array
+                    @ possible_total_weights_array.T,
+                ],
+                [
+                    -weight_constraint_penalty
+                    * possible_total_weights_array
+                    @ weights_array.T,
+                    weight_constraint_penalty
+                    * (possible_total_weights_array @ possible_total_weights_array.T)
+                    + only_one_weight_constraint_penalty,
+                ],
+            ]
+        )
+
+        linear = np.concat(
+            (
+                -costs_array,
+                -2
+                * only_one_weight_constraint_penalty
+                * np.ones(self.__max_weight + 1).reshape(-1, 1),
+            )
+        ).reshape(
+            -1,
+        )
+
+        bias = only_one_weight_constraint_penalty
+
+        return QuadraticPolynomial(quadratic, linear, bias, dtype=dtype, device=device)
+
+    def _from_optimized_tensor(
+        self, optimized_tensor: torch.Tensor, optimized_cost: torch.Tensor
+    ) -> List[int]:
+        return [
+            index
+            for index in range(len(self.__costs))
+            if optimized_tensor[index].item() == 1.0
+        ]
+
+
+def knapsack(
+    costs: List[Union[int, float]], weights: List[int], max_weight: int
+) -> Knapsack:
+    return Knapsack(costs, weights, max_weight)
diff --git a/src/simulated_bifurcation/models/markowitz.py b/src/simulated_bifurcation/models/markowitz.py
index 58e34914..565fef0c 100644
--- a/src/simulated_bifurcation/models/markowitz.py
+++ b/src/simulated_bifurcation/models/markowitz.py
@@ -3,10 +3,11 @@
 import numpy as np
 import torch
 
-from .abc_model import ABCModel
+from ..core.quadratic_polynomial import QuadraticPolynomial
+from .abc_quadratic_model import ABCQuadraticModel
 
 
-class SequentialMarkowitz(ABCModel):
+class SequentialMarkowitz(ABCQuadraticModel):
     """
     Implementation of the Markowitz model for the integer
     trading trajectory optimization problem.
@@ -26,7 +27,7 @@ class SequentialMarkowitz(ABCModel):
     where:
 
     - :math:`1, ..., t, ... T` are the trading decision instants
-    - :math:`w_{t}` is the stocks portolio at time :math:`t`. All the stocks
+    - :math:`w_{t}` is the stocks portfolio at time :math:`t`. All the stocks
       are integer values.
     - :math:`\\Sigma_{t}` is the assets' covariance matrix at time :math:`t`
     - :math:`\\mu_{t}` is the vector of expected returns at time :math:`t`
@@ -35,6 +36,8 @@ class SequentialMarkowitz(ABCModel):
       resulting from the stock difference between time :math:`t - 1` and :math:`t`
     """
 
+    sense = "maximize"
+
     def __init__(
         self,
         covariances: Union[torch.Tensor, np.ndarray],
@@ -43,50 +46,44 @@ def __init__(
         initial_stocks: Optional[Union[torch.Tensor, np.ndarray]] = None,
         risk_coefficient: float = 1,
         number_of_bits: int = 1,
-        dtype: Optional[torch.dtype] = None,
-        device: Optional[Union[str, torch.device]] = None,
     ) -> None:
         """
         Instantiates a Markowitz problem model that can be optimized using
         the Simulated Bifurcation algorithm.
         """
-        self.covariances = self.__cast_to_tensor(
-            covariances, dtype=dtype, device=device
-        )
-        self.expected_returns = self.__cast_to_tensor(
-            expected_returns, dtype=dtype, device=device
-        )
-        self.rebalancing_costs = self.__cast_to_tensor(
-            rebalancing_costs, dtype=dtype, device=device
-        )
+        self.covariances = covariances
+        self.expected_returns = expected_returns
+        self.rebalancing_costs = rebalancing_costs
         self.risk_coefficient = risk_coefficient
         self.timestamps = covariances.shape[0]
         self.assets = covariances.shape[1]
 
         self.initial_stocks = (
-            torch.zeros(self.assets) if initial_stocks is None else initial_stocks
+            torch.zeros(
+                self.assets,
+                dtype=rebalancing_costs.dtype,
+                device=rebalancing_costs.device,
+            )
+            if initial_stocks is None
+            else initial_stocks
         )
 
         self.number_of_bits = number_of_bits
         self.domain = f"int{number_of_bits}"
 
-        matrix, vector, constant = self.compile_model()
-        super().__init__(
-            matrix,
-            vector.reshape(
-                -1,
-            ),
-            constant,
+        self.shape = self.timestamps, self.assets
+
+    def _as_quadratic_polynomial(
+        self, dtype: Optional[torch.dtype], device: Optional[Union[str, torch.device]]
+    ) -> QuadraticPolynomial:
+        return QuadraticPolynomial(
+            self.__compile_matrix(),
+            self.__compile_vector(),
+            self.__compile_offset(),
             dtype=dtype,
             device=device,
         )
 
-    def compile_model(self) -> Tuple[torch.Tensor, torch.Tensor, float]:
-        matrix = self.__compile_matrix()
-        vector = self.__compile_vector()
-        constant = self.__compile_constant()
-        return matrix, vector, constant
-
     def __compile_matrix(self) -> torch.Tensor:
         block_diagonal = self.__create_block_diagonal()
         upper_block_diagonal = self.__create_upper_block_diagonal()
@@ -116,40 +113,19 @@ def __compile_vector(self) -> torch.Tensor:
             first_timestamp_rebalancing_costs + first_timestamp_rebalancing_costs.t()
         ) @ self.initial_stocks.reshape(-1, 1)
         stacked_expected_returns[: self.assets] += initial_stocks_contribution
-        return stacked_expected_returns
+        return stacked_expected_returns.reshape(
+            -1,
+        )
 
-    def __compile_constant(self) -> float:
+    def __compile_offset(self) -> float:
         initial_portfolio = self.initial_stocks.reshape(-1, 1)
-        constant = (
-            -initial_portfolio.t() @ self.rebalancing_costs[0] @ initial_portfolio
-        )
-        return constant.item()
+        constant = -initial_portfolio.T @ self.rebalancing_costs[0] @ initial_portfolio
+        return float(constant)
 
-    def __cast_to_tensor(
-        self,
-        tensor_like: Union[torch.Tensor, np.ndarray],
-        dtype: torch.dtype,
-        device: Union[str, torch.device],
-    ) -> torch.Tensor:
-        if isinstance(tensor_like, torch.Tensor):
-            return tensor_like.to(dtype=dtype, device=device)
-        else:
-            return torch.from_numpy(tensor_like).to(dtype=dtype, device=device)
-
-    @property
-    def portfolio(self) -> Optional[torch.Tensor]:
-        if self.sb_result is None:
-            return None
-        best_agent = torch.argmax(self(self.sb_result.t())).item()
-        return self.sb_result[:, best_agent].reshape(self.timestamps, self.assets)
-
-    @property
-    def gains(self) -> float:
-        return (
-            self(self.portfolio.reshape(1, -1)).item()
-            if self.portfolio is not None
-            else None
-        )
+    def _from_optimized_tensor(
+        self, optimized_tensor: torch.Tensor, optimized_cost: torch.Tensor
+    ) -> Tuple[torch.Tensor, float]:
+        return optimized_tensor.reshape(*self.shape), optimized_cost.item()
 
 
 class Markowitz(SequentialMarkowitz):
@@ -164,8 +140,6 @@ def __init__(
         expected_return: Union[torch.Tensor, np.ndarray],
         risk_coefficient: float = 1,
         number_of_bits: int = 1,
-        dtype: Optional[torch.dtype] = None,
-        device: Optional[Union[str, torch.device]] = None,
     ) -> None:
         covariance = torch.unsqueeze(covariance, 0)
         expected_return = torch.unsqueeze(expected_return, 0)
@@ -177,25 +151,32 @@ def __init__(
             None,
             risk_coefficient,
             number_of_bits,
-            dtype,
-            device,
-        )
-
-    @property
-    def covariance(self) -> torch.Tensor:
-        return -(2 / self.risk_coefficient) * self._quadratic_coefficients
-
-    @property
-    def expected_return(self) -> torch.Tensor:
-        return self._linear_coefficients
-
-    @property
-    def portfolio(self) -> Optional[torch.Tensor]:
-        portfolio = super().portfolio
-        return (
-            portfolio.reshape(
-                -1,
-            )
-            if portfolio is not None
-            else None
         )
+        self.shape = (-1,)
+
+
+def markowitz(
+    covariance: Union[torch.Tensor, np.ndarray],
+    expected_return: Union[torch.Tensor, np.ndarray],
+    risk_coefficient: float = 1,
+    number_of_bits: int = 1,
+) -> Markowitz:
+    return Markowitz(covariance, expected_return, risk_coefficient, number_of_bits)
+
+
+def sequential_markowitz(
+    covariances: Union[torch.Tensor, np.ndarray],
+    expected_returns: Union[torch.Tensor, np.ndarray],
+    rebalancing_costs: Union[torch.Tensor, np.ndarray],
+    initial_stocks: Optional[Union[torch.Tensor, np.ndarray]] = None,
+    risk_coefficient: float = 1,
+    number_of_bits: int = 1,
+) -> SequentialMarkowitz:
+    return SequentialMarkowitz(
+        covariances,
+        expected_returns,
+        rebalancing_costs,
+        initial_stocks,
+        risk_coefficient,
+        number_of_bits,
+    )
diff --git a/src/simulated_bifurcation/models/number_partitioning.py b/src/simulated_bifurcation/models/number_partitioning.py
index 12ca26c2..0f0f1b5a 100644
--- a/src/simulated_bifurcation/models/number_partitioning.py
+++ b/src/simulated_bifurcation/models/number_partitioning.py
@@ -1,54 +1,52 @@
 from typing import Dict, List, Optional, Union
 
+import numpy as np
 import torch
-from numpy import sum
 
-from .abc_model import ABCModel
+from ..core.quadratic_polynomial import QuadraticPolynomial
+from .abc_quadratic_model import ABCQuadraticModel
 
 
-class NumberPartitioning(ABCModel):
+class NumberPartitioning(ABCQuadraticModel):
     """
     A solver that separates a set of numbers into two subsets, the
     respective sums of which are as close as possible.
     """
 
     domain = "spin"
+    sense = "minimize"
 
-    def __init__(
-        self,
-        numbers: list,
-        dtype: Optional[torch.dtype] = None,
-        device: Optional[Union[str, torch.device]] = None,
-    ) -> None:
-        self.numbers = numbers
-        tensor_numbers = torch.tensor(self.numbers, dtype=dtype, device=device).reshape(
-            -1, 1
-        )
-        super().__init__(
-            tensor_numbers @ tensor_numbers.t(), dtype=dtype, device=device
+    def __init__(self, numbers: List[Union[int, float]]) -> None:
+        self.__numbers = numbers[:]
+
+    def _as_quadratic_polynomial(
+        self, dtype: Optional[torch.dtype], device: Optional[Union[str, torch.device]]
+    ) -> QuadraticPolynomial:
+        tensor_numbers = np.array(self.__numbers).reshape(-1, 1)
+        return QuadraticPolynomial(
+            tensor_numbers @ tensor_numbers.T, dtype=dtype, device=device
         )
 
-    @property
-    def partition(self) -> Dict[str, Dict[str, Union[List[int], int, None]]]:
+    def _from_optimized_tensor(
+        self, optimized_tensor: torch.Tensor, optimized_cost: torch.Tensor
+    ) -> Dict:
         result = {
             "left": {"values": [], "sum": None},
             "right": {"values": [], "sum": None},
         }
-        if self.sb_result is None:
-            return result
-
-        i_min = torch.argmin(self(self.sb_result.t())).item()
-        best_vector = self.sb_result[:, i_min]
-
         left_subset = []
         right_subset = []
-        for elt in range(self._dimension):
-            if best_vector[elt].item() > 0:
-                left_subset.append(self.numbers[elt])
+        for elt in range(len(self.__numbers)):
+            if optimized_tensor[elt].item() > 0:
+                left_subset.append(self.__numbers[elt])
             else:
-                right_subset.append(self.numbers[elt])
+                right_subset.append(self.__numbers[elt])
         result["left"]["values"] = left_subset
-        result["left"]["sum"] = sum(left_subset)
+        result["left"]["sum"] = float(np.sum(left_subset))
         result["right"]["values"] = right_subset
-        result["right"]["sum"] = sum(right_subset)
+        result["right"]["sum"] = float(np.sum(right_subset))
         return result
+
+
+def number_partitioning(numbers: List[Union[int, float]]) -> NumberPartitioning:
+    return NumberPartitioning(numbers)
diff --git a/src/simulated_bifurcation/models/qubo.py b/src/simulated_bifurcation/models/qubo.py
deleted file mode 100644
index 7d727e62..00000000
--- a/src/simulated_bifurcation/models/qubo.py
+++ /dev/null
@@ -1,26 +0,0 @@
-from typing import Optional, Union
-
-import numpy as np
-import torch
-
-from .abc_model import ABCModel
-
-
-class QUBO(ABCModel):
-    """
-    Quadratic Unconstrained Binary Optimization
-
-    Given a matrix `Q` the value to minimize is the quadratic form
-    `ΣΣ Q(i,j)b(i)b(j)` where the `b(i)`'s values are either `0` or `1`.
-    """
-
-    domain = "binary"
-
-    def __init__(
-        self,
-        Q: Union[torch.Tensor, np.ndarray],
-        dtype: Optional[torch.dtype] = None,
-        device: Optional[Union[str, torch.device]] = None,
-    ) -> None:
-        super().__init__(Q, dtype=dtype, device=device)
-        self.Q = self._quadratic_coefficients
diff --git a/tests/models/__init__.py b/tests/models/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/tests/models/test_ising_model.py b/tests/models/test_ising_model.py
deleted file mode 100644
index d906b1f7..00000000
--- a/tests/models/test_ising_model.py
+++ /dev/null
@@ -1,17 +0,0 @@
-import torch
-
-from src.simulated_bifurcation.models import Ising
-
-
-def test_ising():
-    torch.manual_seed(42)
-    J = torch.tensor([[0, -2, 3], [-2, 0, 1], [3, 1, 0]])
-    h = torch.tensor([1, -4, 2])
-    model = Ising(J, h, dtype=torch.float32, device=torch.device("cpu"))
-    spin_vector, value = model.minimize(
-        agents=10, best_only=True, mode="ballistic", verbose=False
-    )
-    assert torch.equal(
-        torch.tensor([-1.0, 1.0, -1.0], dtype=torch.float32), spin_vector
-    )
-    assert -11.0 == value
diff --git a/tests/models/test_knapsack.py b/tests/models/test_knapsack.py
new file mode 100644
index 00000000..5382ad34
--- /dev/null
+++ b/tests/models/test_knapsack.py
@@ -0,0 +1,425 @@
+import pytest
+import torch
+
+from src.simulated_bifurcation.models import knapsack
+
+from ..test_utils import DEVICES
+
+
+@pytest.mark.parametrize("device", DEVICES)
+def test_quadratic_model_conversion(device: torch.device):
+    weights = [14, 2, 5, 3, 6]
+    costs = [10.0, 6.0, 7.0, 1.0, 3.0]
+    max_weight = 12
+
+    model = knapsack(costs, weights, max_weight)
+
+    quadratic_polynomial = model._as_quadratic_polynomial(
+        dtype=torch.float32, device=device
+    )
+
+    assert torch.equal(
+        torch.tensor(
+            [
+                [
+                    5488.0,
+                    784.0,
+                    1960.0,
+                    1176.0,
+                    2352.0,
+                    0.0,
+                    -392.0,
+                    -784.0,
+                    -1176.0,
+                    -1568.0,
+                    -1960.0,
+                    -2352.0,
+                    -2744.0,
+                    -3136.0,
+                    -3528.0,
+                    -3920.0,
+                    -4312.0,
+                    -4704.0,
+                ],
+                [
+                    784.0,
+                    112.0,
+                    280.0,
+                    168.0,
+                    336.0,
+                    0.0,
+                    -56.0,
+                    -112.0,
+                    -168.0,
+                    -224.0,
+                    -280.0,
+                    -336.0,
+                    -392.0,
+                    -448.0,
+                    -504.0,
+                    -560.0,
+                    -616.0,
+                    -672.0,
+                ],
+                [
+                    1960.0,
+                    280.0,
+                    700.0,
+                    420.0,
+                    840.0,
+                    0.0,
+                    -140.0,
+                    -280.0,
+                    -420.0,
+                    -560.0,
+                    -700.0,
+                    -840.0,
+                    -980.0,
+                    -1120.0,
+                    -1260.0,
+                    -1400.0,
+                    -1540.0,
+                    -1680.0,
+                ],
+                [
+                    1176.0,
+                    168.0,
+                    420.0,
+                    252.0,
+                    504.0,
+                    0.0,
+                    -84.0,
+                    -168.0,
+                    -252.0,
+                    -336.0,
+                    -420.0,
+                    -504.0,
+                    -588.0,
+                    -672.0,
+                    -756.0,
+                    -840.0,
+                    -924.0,
+                    -1008.0,
+                ],
+                [
+                    2352.0,
+                    336.0,
+                    840.0,
+                    504.0,
+                    1008.0,
+                    0.0,
+                    -168.0,
+                    -336.0,
+                    -504.0,
+                    -672.0,
+                    -840.0,
+                    -1008.0,
+                    -1176.0,
+                    -1344.0,
+                    -1512.0,
+                    -1680.0,
+                    -1848.0,
+                    -2016.0,
+                ],
+                [
+                    0.0,
+                    0.0,
+                    0.0,
+                    0.0,
+                    0.0,
+                    1204.0,
+                    1204.0,
+                    1204.0,
+                    1204.0,
+                    1204.0,
+                    1204.0,
+                    1204.0,
+                    1204.0,
+                    1204.0,
+                    1204.0,
+                    1204.0,
+                    1204.0,
+                    1204.0,
+                ],
+                [
+                    -392.0,
+                    -56.0,
+                    -140.0,
+                    -84.0,
+                    -168.0,
+                    1204.0,
+                    1232.0,
+                    1260.0,
+                    1288.0,
+                    1316.0,
+                    1344.0,
+                    1372.0,
+                    1400.0,
+                    1428.0,
+                    1456.0,
+                    1484.0,
+                    1512.0,
+                    1540.0,
+                ],
+                [
+                    -784.0,
+                    -112.0,
+                    -280.0,
+                    -168.0,
+                    -336.0,
+                    1204.0,
+                    1260.0,
+                    1316.0,
+                    1372.0,
+                    1428.0,
+                    1484.0,
+                    1540.0,
+                    1596.0,
+                    1652.0,
+                    1708.0,
+                    1764.0,
+                    1820.0,
+                    1876.0,
+                ],
+                [
+                    -1176.0,
+                    -168.0,
+                    -420.0,
+                    -252.0,
+                    -504.0,
+                    1204.0,
+                    1288.0,
+                    1372.0,
+                    1456.0,
+                    1540.0,
+                    1624.0,
+                    1708.0,
+                    1792.0,
+                    1876.0,
+                    1960.0,
+                    2044.0,
+                    2128.0,
+                    2212.0,
+                ],
+                [
+                    -1568.0,
+                    -224.0,
+                    -560.0,
+                    -336.0,
+                    -672.0,
+                    1204.0,
+                    1316.0,
+                    1428.0,
+                    1540.0,
+                    1652.0,
+                    1764.0,
+                    1876.0,
+                    1988.0,
+                    2100.0,
+                    2212.0,
+                    2324.0,
+                    2436.0,
+                    2548.0,
+                ],
+                [
+                    -1960.0,
+                    -280.0,
+                    -700.0,
+                    -420.0,
+                    -840.0,
+                    1204.0,
+                    1344.0,
+                    1484.0,
+                    1624.0,
+                    1764.0,
+                    1904.0,
+                    2044.0,
+                    2184.0,
+                    2324.0,
+                    2464.0,
+                    2604.0,
+                    2744.0,
+                    2884.0,
+                ],
+                [
+                    -2352.0,
+                    -336.0,
+                    -840.0,
+                    -504.0,
+                    -1008.0,
+                    1204.0,
+                    1372.0,
+                    1540.0,
+                    1708.0,
+                    1876.0,
+                    2044.0,
+                    2212.0,
+                    2380.0,
+                    2548.0,
+                    2716.0,
+                    2884.0,
+                    3052.0,
+                    3220.0,
+                ],
+                [
+                    -2744.0,
+                    -392.0,
+                    -980.0,
+                    -588.0,
+                    -1176.0,
+                    1204.0,
+                    1400.0,
+                    1596.0,
+                    1792.0,
+                    1988.0,
+                    2184.0,
+                    2380.0,
+                    2576.0,
+                    2772.0,
+                    2968.0,
+                    3164.0,
+                    3360.0,
+                    3556.0,
+                ],
+                [
+                    -3136.0,
+                    -448.0,
+                    -1120.0,
+                    -672.0,
+                    -1344.0,
+                    1204.0,
+                    1428.0,
+                    1652.0,
+                    1876.0,
+                    2100.0,
+                    2324.0,
+                    2548.0,
+                    2772.0,
+                    2996.0,
+                    3220.0,
+                    3444.0,
+                    3668.0,
+                    3892.0,
+                ],
+                [
+                    -3528.0,
+                    -504.0,
+                    -1260.0,
+                    -756.0,
+                    -1512.0,
+                    1204.0,
+                    1456.0,
+                    1708.0,
+                    1960.0,
+                    2212.0,
+                    2464.0,
+                    2716.0,
+                    2968.0,
+                    3220.0,
+                    3472.0,
+                    3724.0,
+                    3976.0,
+                    4228.0,
+                ],
+                [
+                    -3920.0,
+                    -560.0,
+                    -1400.0,
+                    -840.0,
+                    -1680.0,
+                    1204.0,
+                    1484.0,
+                    1764.0,
+                    2044.0,
+                    2324.0,
+                    2604.0,
+                    2884.0,
+                    3164.0,
+                    3444.0,
+                    3724.0,
+                    4004.0,
+                    4284.0,
+                    4564.0,
+                ],
+                [
+                    -4312.0,
+                    -616.0,
+                    -1540.0,
+                    -924.0,
+                    -1848.0,
+                    1204.0,
+                    1512.0,
+                    1820.0,
+                    2128.0,
+                    2436.0,
+                    2744.0,
+                    3052.0,
+                    3360.0,
+                    3668.0,
+                    3976.0,
+                    4284.0,
+                    4592.0,
+                    4900.0,
+                ],
+                [
+                    -4704.0,
+                    -672.0,
+                    -1680.0,
+                    -1008.0,
+                    -2016.0,
+                    1204.0,
+                    1540.0,
+                    1876.0,
+                    2212.0,
+                    2548.0,
+                    2884.0,
+                    3220.0,
+                    3556.0,
+                    3892.0,
+                    4228.0,
+                    4564.0,
+                    4900.0,
+                    5236.0,
+                ],
+            ],
+            dtype=torch.float32,
+            device=device,
+        ),
+        quadratic_polynomial.quadratic,
+    )
+    assert torch.equal(
+        torch.tensor(
+            [
+                -10.0,
+                -6.0,
+                -7.0,
+                -1.0,
+                -3.0,
+                -2408.0,
+                -2408.0,
+                -2408.0,
+                -2408.0,
+                -2408.0,
+                -2408.0,
+                -2408.0,
+                -2408.0,
+                -2408.0,
+                -2408.0,
+                -2408.0,
+                -2408.0,
+                -2408.0,
+            ],
+            dtype=torch.float32,
+            device=device,
+        ),
+        quadratic_polynomial.linear,
+    )
+    assert torch.equal(
+        torch.tensor(1204.0, dtype=torch.float32, device=device),
+        quadratic_polynomial.bias,
+    )
+
+    # Solve and check result
+    assert [1, 2, 3] == model.solve(
+        mode="discrete", agents=2000, device=device, dtype=torch.float32
+    )
diff --git a/tests/models/test_markowitz.py b/tests/models/test_markowitz.py
index a844dd1c..e85acfce 100644
--- a/tests/models/test_markowitz.py
+++ b/tests/models/test_markowitz.py
@@ -1,25 +1,19 @@
 import torch
 
-from src.simulated_bifurcation.models import Markowitz
+from src.simulated_bifurcation.models import markowitz
 
 
 def test_markowitz():
     torch.manual_seed(42)
     covariance = torch.tensor([[1.0, 1.2, 0.7], [1.2, 1.0, -1.9], [0.7, -1.9, 1.0]])
     expected_return = torch.tensor([0.2, 0.05, 0.17])
-    model = Markowitz(
+    model = markowitz(
         covariance,
         expected_return,
         risk_coefficient=1,
         number_of_bits=3,
-        dtype=torch.float32,
     )
 
-    assert torch.all(torch.isclose(covariance, model.covariance))
-    assert torch.equal(expected_return, model.expected_return)
-    assert model.portfolio is None
-    assert model.gains is None
-
-    model.maximize(agents=10, verbose=False)
-    assert torch.equal(torch.tensor([0.0, 7.0, 7.0]), model.portfolio)
-    assert 45.64 == round(model.gains, 2)
+    result = model.solve(agents=10, verbose=False)
+    assert torch.equal(torch.tensor([0.0, 7.0, 7.0]), result[0])
+    assert 45.64 == round(result[1], 2)
diff --git a/tests/models/test_number_partitioning.py b/tests/models/test_number_partitioning.py
index bb8c91d3..3f344869 100644
--- a/tests/models/test_number_partitioning.py
+++ b/tests/models/test_number_partitioning.py
@@ -1,29 +1,19 @@
 import torch
 
-from src.simulated_bifurcation.models import NumberPartitioning
+from src.simulated_bifurcation.models import number_partitioning
 
 
 def test_number_partitioning_with_even_sum():
     torch.manual_seed(42)
     numbers = [2, 7, 3, 8, 1, 5, 4, 6]
-    model = NumberPartitioning(numbers, dtype=torch.float32)
-    model.minimize(agents=10, verbose=False)
-    result = model.partition
+    model = number_partitioning(numbers)
+    result = model.solve(agents=10, verbose=False, dtype=torch.float32)
     assert result["left"]["sum"] == result["right"]["sum"]
 
 
 def test_number_partitioning_with_odd_sum():
     torch.manual_seed(42)
     numbers = [2, 7, 3, 9, 1, 5, 4, 6]
-    model = NumberPartitioning(numbers, dtype=torch.float32)
-    model.minimize(agents=10, verbose=False)
-    result = model.partition
+    model = number_partitioning(numbers)
+    result = model.solve(agents=10, verbose=False, dtype=torch.float32)
     assert abs(result["left"]["sum"] - result["right"]["sum"]) == 1
-
-
-def test_not_optimized_model():
-    model = NumberPartitioning([1, 2, 3])
-    assert model.partition == {
-        "left": {"values": [], "sum": None},
-        "right": {"values": [], "sum": None},
-    }
diff --git a/tests/models/test_qubo.py b/tests/models/test_qubo.py
deleted file mode 100644
index 2d3ff1c1..00000000
--- a/tests/models/test_qubo.py
+++ /dev/null
@@ -1,37 +0,0 @@
-import torch
-
-from src.simulated_bifurcation.models import QUBO
-
-
-def test_qubo():
-    torch.manual_seed(42)
-    Q = torch.tensor([[1, -2, 3], [0, -4, 1], [0, 0, 2]])
-    model = QUBO(Q, dtype=torch.float32)
-    binary_vector, value = model.minimize(agents=10, verbose=False, best_only=True)
-    assert torch.equal(torch.tensor([1.0, 1.0, 0.0]), binary_vector)
-    assert -5.0 == value
-    binary_vector, value = model.maximize(agents=10, verbose=False, best_only=True)
-    assert torch.equal(torch.tensor([1.0, 0.0, 1.0]), binary_vector)
-    assert 6.0 == value
-
-
-def test_lp_problem_formulated_as_qubo():
-    torch.manual_seed(42)
-    P = 15.5
-    Q = torch.tensor(
-        [
-            [2, -P, -P, 0, 0, 0],
-            [0, 2, -P, -P, 0, 0],
-            [0, 0, 2, -2 * P, 0, 0],
-            [0, 0, 0, 2, -P, 0],
-            [0, 0, 0, 0, 4.5, -P],
-            [0, 0, 0, 0, 0, 3],
-        ]
-    )
-    binary_vector, objective_value = QUBO(Q, dtype=torch.float32).maximize(
-        agents=10, verbose=False
-    )
-    assert torch.equal(
-        torch.tensor([1.0, 0.0, 0.0, 1.0, 0.0, 1.0], dtype=torch.float32), binary_vector
-    )
-    assert 7.0 == objective_value
diff --git a/tests/models/test_sequential_markowitz.py b/tests/models/test_sequential_markowitz.py
index 8467f685..bbd408c1 100644
--- a/tests/models/test_sequential_markowitz.py
+++ b/tests/models/test_sequential_markowitz.py
@@ -1,7 +1,6 @@
-import numpy as np
 import torch
 
-from src.simulated_bifurcation.models import SequentialMarkowitz
+from src.simulated_bifurcation.models import sequential_markowitz
 
 #  2 assets on 4 timestamps case test
 initial_stocks = torch.tensor([0.0, 1.0])
@@ -33,60 +32,18 @@
 
 def test_sequential_markowitz():
     torch.manual_seed(42)
-    model = SequentialMarkowitz(
+    model = sequential_markowitz(
         covariances,
         expected_returns,
-        np.array(
-            [
-                [[0.1, 0.0], [0.0, 0.2]],
-                [[0.11, 0.0], [0.0, 0.2]],
-                [[0.05, 0.0], [0.0, 0.4]],
-                [[0.01, 0.0], [0.0, 0.5]],
-            ]
-        ),
+        rebalancing_costs,
         initial_stocks,
         risk_coefficient=1,
         number_of_bits=1,
-        dtype=torch.float32,
     )
 
-    assert torch.equal(covariances, model.covariances)
-    assert torch.equal(expected_returns, model.expected_returns)
-    assert torch.equal(rebalancing_costs, model.rebalancing_costs)
-    assert torch.equal(initial_stocks, model.initial_stocks)
-    assert model.portfolio is None
-    assert model.gains is None
-
-    assert torch.all(
-        torch.isclose(
-            torch.tensor(
-                [
-                    [-0.71, -0.10, 0.11, 0.0, 0.0, 0.0, 0.0, 0.0],
-                    [-0.10, -0.90, 0.0, 0.20, 0.0, 0.0, 0.0, 0.0],
-                    [0.11, 0.0, -0.66, -0.15, 0.05, 0.0, 0.0, 0.0],
-                    [0.0, 0.20, -0.15, -1.10, 0.0, 0.40, 0.0, 0.0],
-                    [0.0, 0.0, 0.05, 0.0, -0.56, -0.175, 0.01, 0.0],
-                    [0.0, 0.0, 0.0, 0.40, -0.175, -1.40, 0.0, 0.50],
-                    [0.0, 0.0, 0.0, 0.0, 0.01, 0.0, -0.51, 0.06],
-                    [0.0, 0.0, 0.0, 0.0, 0.0, 0.50, 0.06, -1.0],
-                ]
-            ),
-            model.quadratic,
-        )
-    )
-    assert torch.all(
-        torch.isclose(
-            torch.tensor(
-                [0.5000, 1.0000, 0.5500, 0.5500, 0.7000, 0.5000, 1.4000, 0.2000]
-            ),
-            model.linear,
-        )
-    )
-    assert torch.equal(torch.tensor(-0.2), model.bias)
-
-    model.maximize(agents=128, early_stopping=False, verbose=False)
-    assert (4, 2) == model.portfolio.shape
+    result = model.solve(agents=128, early_stopping=False, verbose=False)
+    assert (4, 2) == result[0].shape
     assert torch.equal(
-        torch.tensor([[0.0, 1.0], [0.0, 0.0], [1.0, 0.0], [1.0, 0.0]]), model.portfolio
+        torch.tensor([[0.0, 1.0], [0.0, 0.0], [1.0, 0.0], [1.0, 0.0]]), result[0]
     )
-    assert 0.95 == round(model.gains, 4)
+    assert 0.95 == round(result[1], 2)