diff --git a/docs/Project.toml b/docs/Project.toml
index 55ac6c65..bd422b7d 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -5,7 +5,10 @@ Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
 ParametricOptInterface = "0ce4ce61-57bf-432b-a095-efac525d185e"
-Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 
 [compat]
 Documenter = "1"
+HiGHS = "1"
+Ipopt = "1"
+JuMP = "1"
+MathOptInterface = "1.13.2"
diff --git a/docs/make.jl b/docs/make.jl
index 7ec8d039..78432e32 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -3,8 +3,8 @@
 # Use of this source code is governed by an MIT-style license that can be found
 # in the LICENSE.md file or at https://opensource.org/licenses/MIT.
 
-using Documenter
-using ParametricOptInterface
+import Documenter
+import ParametricOptInterface
 
 Documenter.makedocs(;
     modules = [ParametricOptInterface],
@@ -19,12 +19,6 @@ Documenter.makedocs(;
     authors = "Tomás Gutierrez, and contributors",
     pages = [
         "Home" => "index.md",
-        "manual.md",
-        "Examples" => [
-            "Examples/example.md",
-            "Examples/benders.md",
-            "Examples/markowitz.md",
-        ],
         "reference.md",
     ],
     checkdocs = :none,
diff --git a/docs/src/Examples/benders.md b/docs/src/Examples/benders.md
deleted file mode 100644
index 2fb10ace..00000000
--- a/docs/src/Examples/benders.md
+++ /dev/null
@@ -1,404 +0,0 @@
-# Benders Quantile Regression
-
-We will apply Norm-1 regression to the [Linear Regression](https://en.wikipedia.org/wiki/Linear_regression) problem.
-Linear regression is a statistical tool to obtain the relation between one **dependent variable** and other **explanatory variables**.
-In other words, given a set of $n$ explanatory variables $X = \{ X_1, \dots, X_n \}$
-we would like to obtain the best possible estimate for $Y$.
-In order to accomplish such a task we make the hypothesis that $Y$
-is approximately linear function of $X$:
-
-$$Y = \sum_{j =1}^n \beta_j X_j + \varepsilon$$
-
-where $\varepsilon$ is some random error.
-
-The estimation of the $\beta$ values relies on observations of the variables:
-$\{y^i, x_1^i, \dots, x_n^i\}_i$.
-
-In this example we will solve a problem where the explanatory variables are sinusoids of differents frequencies.
-First, we define the number of explanatory variables and observations
-
-```julia
-using ParametricOptInterface,MathOptInterface,JuMP,HiGHS
-using TimerOutputs,LinearAlgebra,Random
-
-const POI = ParametricOptInterface
-const MOI = MathOptInterface
-const OPTIMIZER = HiGHS.Optimizer;
-
-const N_Candidates = 200
-const N_Observations = 2000
-const N_Nodes = 200
-
-const Observations = 1:N_Observations
-const Candidates = 1:N_Candidates
-const Nodes = 1:N_Nodes;
-```
-
-
-Initialize a random number generator to keep results deterministic
-```julia
-rng = Random.MersenneTwister(123);
-```
-
-Building regressors (explanatory) sinusoids
-
-```julia
-const X = zeros(N_Candidates, N_Observations)
-const time = [obs / N_Observations * 1 for obs in Observations]
-for obs in Observations, cand in Candidates
-    t = time[obs]
-    f = cand
-    X[cand, obs] = sin(2 * pi * f * t)
-end
-```
-
-Define coefficients
-
-```julia
-β = zeros(N_Candidates)
-for i in Candidates
-    if rand(rng) <= (1 - i / N_Candidates)^2 && i <= 100
-        β[i] = 4 * rand(rng) / i
-    end
-end
-```
-
-Create noisy observations
-
-```julia
-const y = X' * β .+ 0.1 * randn(rng, N_Observations)
-```
-
-### Benders Decomposition
-
-Benders decomposition is used to solve large optimization problems with some special characteristics.
-LP's can be solved with classical linear optimization methods
-such as the Simplex method or Interior point methods provided by
-solvers like HiGHS.
-However, these methods do not scale linearly with the problem size.
-In the Benders decomposition framework we break the problem in two pieces:
-A outer and a inner problem.
-Of course some variables will belong to both problems, this is where the
-cleverness of Benders kicks in:
-The outer problem is solved and passes the shared variables to the inner.
-The inner problem is solved with the shared variables FIXED to the values
-given by the outer problem. The solution of the inner problem can be used
-to generate a constraint to the outer problem to describe the linear
-approximation of the cost function of the shared variables.
-In many cases, like stochastic programming, the inner problems have a
-interesting structure and might be broken in smaller problem to be solved
-in parallel.
-
-We will descibe the decomposition similarly to what is done in:
-Introduction to Linear Optimization, Bertsimas & Tsitsiklis (Chapter 6.5):
-Where the problem in question has the form
-
-$$\begin{align}
-   & \min_{x, y_k}     &&  c^T x && + f_1^T y_1 && + \dots && + f_n^T y_n  &&  \notag \\
-   & \text{subject to} &&  Ax    &&             &&         &&              && = b \notag \\
-   &                   &&  B_1 x && + D_1 y_1   &&         &&              && = d_1 \notag \\
-   &                   &&  \dots &&             &&  \dots  &&              &&       \notag \\
-   &                   &&  B_n x &&             &&         && + D_n y_n    && = d_n \notag \\
-   &                   &&   x,   &&     y_1,    &&         &&       y_n    && \geq 0 \notag \\
-\end{align}$$
-
-### Inner Problem
-
-Given a solution for the $x$ variables we can define the inner problem as
-
-$$\begin{align}
- z_k(x) \ = \ & \min_{y_k}        &&  f_k^T y_k &&                  \notag \\
-              & \text{subject to} &&  D_k y_k   &&  = d_k - B_k x \notag \\
-              &                   &&  y_k       && \geq 0          \notag \\
-\end{align}$$
-
-The $z_k(x)$ function represents the cost of the subproblem given a
-solution for $x$. This function is a convex function because $x$
-affects only the right hand side of the problem (this is a standard
-results in LP theory).
-
-For the special case of the Norm-1 reggression the problem is written as:
-
-$$\begin{align}
-z_k(\beta) \ = \ & \min_{\varepsilon^{up}, \varepsilon^{dw}}         &&  \sum_{i \in ObsSet(k)} {\varepsilon^{up}}_i + {\varepsilon^{dw}}_i               && \notag \\
-                 & \text{subject to}     &&  {\varepsilon^{up}}_i \geq + y_i - \sum_{j \in Candidates} \beta_j x_{i,j} && \forall i \in ObsSet(k) \notag \\
-                 &                       &&  {\varepsilon^{dw}}_i \geq - y_i + \sum_{j \in Candidates} \beta_j x_{i,j} && \forall i \in ObsSet(k) \notag \\
-                 &                       &&  {\varepsilon^{up}}_i, {\varepsilon^{dw}}_i \geq 0                             && \forall i \in ObsSet(k) \notag \\
-\end{align}$$
-
-The collection $ObsSet(k)$ is a sub-set of the `N_Observations`.
-Any partition of the `N_Observations` collection is valid.
-In this example we will partition with the function:
-
-```julia
-function ObsSet(K)
-    obs_per_block = div(N_Observations, N_Nodes)
-    return (1+(K-1)*obs_per_block):(K*obs_per_block)
-end
-```
-
-Which can be written in POI as follows:
-
-```julia
-function inner_model(K)
-
-    # initialize the POI model
-    inner = direct_model(POI.Optimizer(OPTIMIZER()))
-
-    # Define local optimization variables for norm-1 error
-    @variables(inner, begin
-        ɛ_up[ObsSet(K)] >= 0
-        ɛ_dw[ObsSet(K)] >= 0
-    end)
-
-    # create the regression coefficient representation
-    # Create parameters
-    β = [@variable(inner, set = MOI.Parameter(0.0)) for i in 1:N_Candidates]
-    for (i, βi) in enumerate(β)
-        set_name(βi, "β[$i]")
-    end
-
-    # create local constraints
-    # Note that *parameter* algebra is implemented just like variables
-    # algebra. We can multiply parameters by constants, add parameters,
-    # sum parameters and variables and so on.
-    @constraints(
-        inner,
-        begin
-            ɛ_up_ctr[i in ObsSet(K)],
-            ɛ_up[i] >= +sum(X[j, i] * β[j] for j in Candidates) - y[i]
-            ɛ_dw_ctr[i in ObsSet(K)],
-            ɛ_dw[i] >= -sum(X[j, i] * β[j] for j in Candidates) + y[i]
-        end
-    )
-
-    # create local objective function
-    @objective(inner, Min, sum(ɛ_up[i] + ɛ_dw[i] for i in ObsSet(K)))
-
-    # return the correct group of parameters
-    return (inner, β)
-end
-```
-
-### Outer Problem
-
-Now that all pieces of the original problem can be representad by
-the convex $z_k(x)$ functions we can recast the problem in the the equivalent form:
-
-$$\begin{align}
-  & \min_{x}          &&  c^T x + z_1(x) + \dots + z_n(x) && \notag \\
-  & \text{subject to} &&  Ax = b                          && \notag \\
-  &                   &&  x \geq 0                        && \notag \\
-\end{align}$$
-
-However we cannot pass a problem in this form to a linear programming
-solver (it could be passed to other kinds of solvers).
-
-Another standart result of optimization theory is that a convex function
-can be represented by its supporting hyper-planes:
-
-$$\begin{align}
-  z_k(x) \ = \ & \min_{z, x}       &&  z && \notag \\
-               & \text{subject to} &&  z \geq \pi_k(\hat{x}) (x - \hat{x}) + z_k(\hat{x}), \ \forall \hat{x} \in dom(z_k) && \notag \\
-\end{align}$$
-
-Then we can re-write (again) the outer problem as
-
-$$\begin{align}
-  & \min_{x, z_k}     &&  c^T x + z_1 + \dots + z_n \notag \\
-  & \text{subject to} &&  z_i \geq \pi_i(\hat{x}) (x - \hat{x}) + z_i(\hat{x}), \ \forall \hat{x} \in dom(z_i), i \in \{1, \dots, n\} \notag \\
-  &                   &&  Ax = b \notag \\
-  &                   &&  x \geq 0 \notag \\
-\end{align}$$
-
-Which is a linear program! However, it has infinitely many constraints !!
-
-We can relax the infinite constraints and write:
-
-$$\begin{align}
-  & \min_{x, z_k}     &&  c^T x + z_1 + \dots + z_n \notag \\
-  & \text{subject to} &&  Ax = b \notag \\
-  &                   &&  x \geq 0 \notag \\
-\end{align}$$
-
-But now its only an underestimated problem.
-In the case of our problem it can be written as:
-
-$$\begin{align}
-  & \min_{\varepsilon, \beta} &&  \sum_{i \in Nodes} \varepsilon_i \notag \\
-  & \text{subject to} &&  \varepsilon_i \geq 0 \notag \\
-\end{align}$$
-
-This model can be written in JuMP:
-
-```julia
-function outer_model()
-    outer = Model(OPTIMIZER)
-    @variables(outer, begin
-        ɛ[Nodes] >= 0
-        β[1:N_Candidates]
-    end)
-    @objective(outer, Min, sum(ɛ[i] for i in Nodes))
-    sol = zeros(N_Candidates)
-    return (outer, ɛ, β, sol)
-end
-```
-
-The method to solve the outer problem and query its solution is given here:
-
-```julia
-function outer_solve(outer_model)
-    model = outer_model[1]
-    β = outer_model[3]
-    optimize!(model)
-    return (value.(β), objective_value(model))
-end
-```
-
-### Supporting Hyperplanes
-
-With these building blocks in hand, we can start building the algorithm.
-So far we know how to:
-- Solve the relaxed outer problem
-- Obtain the solution for the $\hat{x}$ (or $\beta$ in our case)
-
-
-Now we can:
-- Fix the values of $\hat{x}$ in the inner problems
-- Solve the inner problems
-- query the solution of the inner problems to obtain the supporting hyperplane
-
-the value of $z_k(\hat{x})$, which is the objective value of the inner problem
-
-and the derivative $\pi_k(\hat{x}) = \frac{d z_k(x)}{d x} \Big|_{x = \hat{x}}$
-The derivative is the dual variable associated to the variable $\hat{x}$,
-which results by applying the chain rule on the constraints duals.
-These new steps are executed by the function:
-
-```julia
-function inner_solve(model, outer_solution)
-    β0 = outer_solution[1]
-    inner = model[1]
-
-    # The first step is to fix the values given by the outer problem
-    @timeit "fix" begin
-        β = model[2]
-        MOI.set.(inner, POI.ParameterValue(), β, β0)
-    end
-
-    # here the inner problem is solved
-    @timeit "opt" optimize!(inner)
-
-    # query dual variables, which are sensitivities
-    # They represent the subgradient (almost a derivative)
-    # of the objective function for infinitesimal variations
-    # of the constants in the linear constraints
-    # POI: we can query dual values of *parameters*
-    π = MOI.get.(inner, POI.ParameterDual(), β)
-
-    # π2 = shadow_price.(β_fix)
-    obj = objective_value(inner)
-    rhs = obj - dot(π, β0)
-    return (rhs, π, obj)
-end
-```
-
-Now that we have cutting plane in hand we can add them to the outer problem
-
-```julia
-function outer_add_cut(outer_model, cut_info, node)
-    outer = outer_model[1]
-    ɛ = outer_model[2]
-    β = outer_model[3]
-
-    rhs = cut_info[1]
-    π = cut_info[2]
-
-    @constraint(outer, ɛ[node] >= sum(π[j] * β[j] for j in Candidates) + rhs)
-end
-```
-
-### Algorithm wrap up
-
-The complete algorithm is
-
-- Solve the relaxed master problem
-- Obtain the solution for the $\hat{x}$ (or $\beta$ in our case)
-- Fix the values of $\hat{x}$ in the slave problems
-- Solve the slave problem
-- query the solution of the slave problem to obtain the supporting hyperplane
-- add hyperplane to master problem
-- repeat
-
-Now we grab all the pieces that we built and we write the benders
-algorithm by calling the above function in a proper order.
-
-The macros `@timeit` are use to time each step of the algorithm.
-
-```julia
-function decomposed_model(;print_timer_outputs::Bool = true)
-    reset_timer!() # reset timer fo comparision
-    time_init = @elapsed @timeit "Init" begin
-        # Create the outer problem with no cuts
-        @timeit "outer" outer = outer_model()
-
-        # initialize solution for the regression coefficients in zero
-        @timeit "Sol" solution = (zeros(N_Candidates), Inf)
-        best_sol = deepcopy(solution)
-
-        # Create the inner problems
-        @timeit "inners" inners =
-            [inner_model(i) for i in Candidates]
-
-        # Save initial version of the inner problems and create
-        # the first set of cuts
-        @timeit "Cuts" cuts =
-            [inner_solve(inners[i], solution) for i in Candidates]
-    end
-
-    UB = +Inf
-    LB = -Inf
-
-    # println("Initialize Iterative step")
-    time_loop = @elapsed @timeit "Loop" for k in 1:80
-
-        # Add cuts generated from each inner problem to the outer problem
-        @timeit "add cuts" for i in Candidates
-            outer_add_cut(outer, cuts[i], i)
-        end
-
-        # Solve the outer problem with the new set of cuts
-        # Obtain new solution candidate for the regression coefficients
-        @timeit "solve outer" solution = outer_solve( outer)
-
-        # Pass the new candidate solution to each of the inner problems
-        # Solve the inner problems and obtain cutting planes
-        @timeit "solve nodes" for i in Candidates
-            cuts[i] = inner_solve( inners[i], solution)
-        end
-
-        LB = solution[2]
-        new_UB = sum(cuts[i][3] for i in Candidates)
-        if new_UB <= UB
-            best_sol = deepcopy(solution)
-        end
-        UB = min(UB, new_UB)
-
-        if abs(UB - LB) / (abs(UB) + abs(LB)) < 0.05
-            break
-        end
-    end
-
-    print_timer_outputs && print_timer()
-
-    return best_sol[1]
-end
-```
-
-Run benders decomposition with POI
-
-```julia
-β2 = decomposed_model(; print_timer_outputs = false);
-GC.gc()
-β2 = decomposed_model();
-```
\ No newline at end of file
diff --git a/docs/src/Examples/example.md b/docs/src/Examples/example.md
deleted file mode 100644
index 7debb5e6..00000000
--- a/docs/src/Examples/example.md
+++ /dev/null
@@ -1,455 +0,0 @@
-# Basic Examples
-
-## MOI example - step by step usage
-
-Let's write a step-by-step example of `POI` usage at the MOI level.
-
-First, we declare a [`ParametricOptInterface.Optimizer`](@ref) on top of a `MOI`
-optimizer. In the example, we consider `HiGHS` as the underlying solver:
-
-```@example moi1
-import HiGHS
-import MathOptInterface as MOI
-import ParametricOptInterface as POI
-optimizer = POI.Optimizer(HiGHS.Optimizer())
-```
-
-We declare the variable `x` as in a typical `MOI` model, and we add a
-non-negativity constraint:
-
-```@example moi1
-x = MOI.add_variables(optimizer, 2)
-for x_i in x
-    MOI.add_constraint(optimizer, x_i, MOI.GreaterThan(0.0))
-end
-```
-
-Now, let's consider 3 `MOI.Parameter`. Two of them, `y`, `z`, will be placed in
-the constraints and one, `w`, in the objective function. We'll start all three
-of them with a value equal to `0`:
-
-```@example moi1
-w, cw = MOI.add_constrained_variable(optimizer, MOI.Parameter(0.0))
-y, cy = MOI.add_constrained_variable(optimizer, MOI.Parameter(0.0))
-z, cz = MOI.add_constrained_variable(optimizer, MOI.Parameter(0.0))
-```
-
-Let's add the constraints. Notice that we treat parameters and variables in the
-same way when building the functions that will be placed in some set to create a
-constraint (`Function-in-Set`):
-
-```@example moi1
-cons1 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2.0, 1.0, 3.0], [x[1], x[2], y]), 0.0)
-ci1 = MOI.add_constraint(optimizer, cons1, MOI.LessThan(4.0))
-cons2 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 2.0, 0.5], [x[1], x[2], z]), 0.0)
-ci2 = MOI.add_constraint(optimizer, cons2, MOI.LessThan(4.0))
-```
-
-Finally, we declare and add the objective function, with its respective sense:
-
-```@example moi1
-obj_func = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([4.0, 3.0, 2.0], [x[1], x[2], w]), 0.0)
-MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), obj_func)
-MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)
-```
-
-Now we can optimize the model and assess its termination and primal status:
-
-```@example moi1
-MOI.optimize!(optimizer)
-MOI.get(optimizer, MOI.TerminationStatus())
-MOI.get(optimizer, MOI.PrimalStatus())
-```
-
-Given the optimized solution, we check that its value is, as expected, equal to
-`28/3`, and the solution vector `x` is `[4/3, 4/3]`:
-
-```@example moi1
-isapprox(MOI.get(optimizer, MOI.ObjectiveValue()), 28/3, atol = 1e-4)
-isapprox(MOI.get(optimizer, MOI.VariablePrimal(), x[1]), 4/3, atol = 1e-4)
-isapprox(MOI.get(optimizer, MOI.VariablePrimal(), x[2]), 4/3, atol = 1e-4)
-```
-
-We can also retrieve the dual values associated to each parameter,
-**as they are all additive**:
-
-```@example moi1
-MOI.get(optimizer, MOI.ConstraintDual(), cy)
-MOI.get(optimizer, MOI.ConstraintDual(), cz)
-MOI.get(optimizer, MOI.ConstraintDual(), cw)
-```
-
-Notice the direct relationship in this case between the parameters' duals and
-the associated constraints' duals.
-
-The  `y` parameter, for example, only appears in the `cons1`. If we compare
-their duals, we can check that the dual of `y` is equal to its coefficient in
-`cons1` multiplied by the constraint's dual itself, as expected:
-
-```@example moi1
-isapprox(
-    MOI.get(optimizer, MOI.ConstraintDual(), cy),
-    3*MOI.get(optimizer, MOI.ConstraintDual(), ci1);
-    atol = 1e-4,
-)
-```
-
-The same is valid for the remaining parameters. In case a parameter appears in
-more than one constraint, or both some constraints and in the objective
-function, its dual will be equal to the linear combination of the functions'
-duals multiplied by the respective coefficients.
-
-So far, we only added some parameters that had no influence at first in solving
-the model. Let's change the values associated to each parameter to assess its
-implications.
-
-First, we set the value of parameters `y` and `z` to `1.0`. Notice that we are
-changing the feasible set of the decision variables:
-
-```@example moi1
-MOI.set(optimizer, POI.ParameterValue(), y, 1.0)
-MOI.set(optimizer, POI.ParameterValue(), z, 1.0)
-```
-
-However, if we check the optimized model now, there will be no changes in the
-objective function value or the in the optimized decision variables:
-
-```@example moi1
-isapprox.(MOI.get(optimizer, MOI.ObjectiveValue()), 28/3, atol = 1e-4)
-isapprox.(MOI.get(optimizer, MOI.VariablePrimal(), x[1]), 4/3, atol = 1e-4)
-isapprox.(MOI.get(optimizer, MOI.VariablePrimal(), x[2]), 4/3, atol = 1e-4)
-```
-
-Although we changed the parameter values, we didn't optimize the model yet.
-Thus, **to apply the parameters' changes, the model must be optimized again**:
-
-```@example moi1
-MOI.optimize!(optimizer)
-```
-
-The `MOI.optimize!()` function handles the necessary updates, properly fowarding
-the new outer model (`POI` model) additions to the inner model (`MOI` model)
-which will be handled by the solver. Now we can assess the updated optimized
-information:
-
-```@example moi1
-isapprox.(MOI.get(optimizer, MOI.ObjectiveValue()), 3.0, atol = 1e-4)
-MOI.get.(optimizer, MOI.VariablePrimal(), x) == [0.0, 1.0]
-```
-
-If we update the parameter `w`, associated to the objective function, we are
-simply adding a constant to it. Notice how the new objective function is
-precisely equal to the previous one plus the new value of `w`. In addition, as
-we didn't update the feasible set, the optimized decision variables remain the
-same.
-
-```@example moi1
-MOI.set(optimizer, POI.ParameterValue(), w, 2.0)
-# Once again, the model must be optimized to incorporate the changes
-MOI.optimize!(optimizer)
-# Only the objective function value changes
-isapprox.(MOI.get(optimizer, MOI.ObjectiveValue()), 7.0, atol = 1e-4)
-MOI.get.(optimizer, MOI.VariablePrimal(), x) == [0.0, 1.0]
-```
-
-## JuMP Example - step by step usage
-
-Let's write a step-by-step example of `POI` usage at the JuMP level.
-
-First, we declare a `Model` on top of a `Optimizer` of an underlying solver. In
-the example, we consider `HiGHS` as the underlying solver:
-
-```@example jump1
-using HiGHS
-using JuMP
-
-using ParametricOptInterface
-const POI = ParametricOptInterface
-
-model = Model(() -> ParametricOptInterface.Optimizer(HiGHS.Optimizer()))
-```
-
-We declare the variable `x` as in a typical `JuMP` model:
-
-```@example jump1
-@variable(model, x[i = 1:2] >= 0)
-```
-
-Now, let's consider 3 `MOI.Parameter`. Two of them, `y`, `z`, will be placed in
-the constraints and one, `w`, in the objective function. We'll start all three
-of them with a value equal to `0`:
-
-```@example jump1
-@variable(model, y in MOI.Parameter(0.0))
-@variable(model, z in MOI.Parameter(0.0))
-@variable(model, w in MOI.Parameter(0.0))
-```
-
-Let's add the constraints. Notice that we treat parameters the same way we treat
-variables when writing the model:
-
-```@example jump1
-@constraint(model, c1, 2x[1] + x[2] + 3y <= 4)
-@constraint(model, c2, x[1] + 2x[2] + 0.5z <= 4)
-```
-
-Finally, we declare and add the objective function, with its respective sense:
-
-```@example jump1
-@objective(model, Max, 4x[1] + 3x[2] + 2w)
-```
-
-We can optimize the model and assess its termination and primal status:
-
-```@example jump1
-optimize!(model)
-termination_status(model)
-primal_status(model)
-```
-
-Given the optimized solution, we check that its value is, as expected, equal to
-`28/3`, and the solution vector `x` is `[4/3, 4/3]`:
-
-```@example jump1
-isapprox(objective_value(model), 28/3)
-isapprox(value.(x), [4/3, 4/3])
-```
-
-We can also retrieve the dual values associated to each parameter, **as they are all additive**:
-
-```@example jump1
-MOI.get(model, POI.ParameterDual(), y)
-MOI.get(model, POI.ParameterDual(), z)
-MOI.get(model, POI.ParameterDual(), w)
-```
-
-Notice the direct relationship in this case between the parameters' duals and the associated constraints' duals. The `y` parameter, for example, only appears in the `c1`. If we compare their duals, we can check that the dual of `y` is equal to its coefficient in `c1` multiplied by the constraint's dual itself, as expected:
-
-```@example jump1
-dual_of_y = MOI.get(model, POI.ParameterDual(), y)
-isapprox(dual_of_y, 3 * dual(c1))
-```
-
-The same is valid for the remaining parameters. In case a parameter appears in more than one constraint, or both some constraints and in the objective function, its dual will be equal to the linear combination of the functions' duals multiplied by the respective coefficients.
-
-So far, we only added some parameters that had no influence at first in solving the model. Let's change the values associated to each parameter to assess its implications. First, we set the value of parameters `y` and `z` to `1.0`. Notice that we are changing the feasible set of the decision variables:
-
-```@example jump1
-MOI.set(model, POI.ParameterValue(), y, 1)
-MOI.set(model, POI.ParameterValue(), z, 1)
-# We can also query the value in the parameters
-MOI.get(model, POI.ParameterValue(), y)
-MOI.get(model, POI.ParameterValue(), z)
-```
-
-To apply the parameters' changes, the model must be optimized again:
-
-```@example jump1
-optimize!(model)
-```
-
-The `optimize!` function handles the necessary updates, properly fowarding the new outer model (`POI` model) additions to the inner model (`MOI` model) which will be handled by the solver. Now we can assess the updated optimized information:
-
-```@example jump1
-isapprox(objective_value(model), 3)
-isapprox(value.(x), [0, 1])
-```
-
-If we update the parameter `w`, associated to the objective function, we are simply adding a constant to it. Notice how the new objective function is precisely equal to the previous one plus the new value of `w`. In addition, as we didn't update the feasible set, the optimized decision variables remain the same.
-
-```@example jump1
-MOI.set(model, POI.ParameterValue(), w, 2)
-# Once again, the model must be optimized to incorporate the changes
-optimize!(model)
-# Only the objective function value changes
-isapprox(objective_value(model), 7)
-isapprox(value.(x), [0, 1])
-```
-
-## JuMP Example - Declaring vectors of parameters
-
-Many times it is useful to declare a vector of parameters just like we declare a vector of variables, the JuMP syntax for variables works with parameters too:
-
-
-```@example jump2
-using HiGHS
-using JuMP
-using ParametricOptInterface
-const POI = ParametricOptInterface
-
-model = Model(() -> ParametricOptInterface.Optimizer(HiGHS.Optimizer()))
-@variable(model, x[i = 1:3] >= 0)
-@variable(model, p1[i = 1:3] in MOI.Parameter(0.0))
-@variable(model, p2[i = 1:3] in MOI.Parameter.([1, 10, 45]))
-@variable(model, p3[i = 1:3] in MOI.Parameter.(ones(3)))
-```
-
-## JuMP Example - Dealing with parametric expressions as variable bounds
-
-A very common pattern that appears when using ParametricOptInterface is to add variable and later add some expression with parameters that represent the variable bound. The following code illustrates the pattern:
-
-```@example jump3
-using HiGHS
-using JuMP
-using ParametricOptInterface
-const POI = ParametricOptInterface
-
-model = direct_model(POI.Optimizer(HiGHS.Optimizer()))
-@variable(model, x)
-@variable(model, p in MOI.Parameter(0.0))
-@constraint(model, x >= p)
-```
-
-Since parameters are treated like variables JuMP lowers this to MOI as `x - p >= 0` which is not a variable bound but a linear constraint.This means that the current representation of this problem at the solver level is:
-
-```math
-\begin{align}
-    & \min_{x} & 0
-    \\
-    & \;\;\text{s.t.} & x & \in \mathbb{R} \\
-    &   & x - p & \geq 0
-\end{align}
-```
-
-This behaviour might be undesirable because it creates extra rows in your problem. Users can set the [`ParametricOptInterface.ConstraintsInterpretation`](@ref) to control how the linear constraints should be interpreted. The pattern advised for users seeking the most performance out of ParametricOptInterface should use the followig pattern:
-
-```@example jump3
-using HiGHS
-using JuMP
-using ParametricOptInterface
-const POI = ParametricOptInterface
-
-model = direct_model(POI.Optimizer(HiGHS.Optimizer()))
-@variable(model, x)
-@variable(model, p in MOI.Parameter(0.0))
-
-# Indicate that all the new constraints will be valid variable bounds
-MOI.set(model, POI.ConstraintsInterpretation(), POI.ONLY_BOUNDS)
-@constraint(model, x >= p)
-# The name of this constraint was different to inform users that this is a
-# variable bound.
-
-# Indicate that all the new constraints will not be variable bounds
-MOI.set(model, POI.ConstraintsInterpretation(), POI.ONLY_CONSTRAINTS)
-# @constraint(model, ...)
-```
-
-This way the mathematical representation of the problem will be:
-
-```math
-\begin{align}
-    & \min_{x} & 0
-    \\
-    & \;\;\text{s.t.} & x & \geq p
-\end{align}
-```
-
-which might lead to faster solves.
-
-Users that just want everything to work can use the default value `POI.ONLY_CONSTRAINTS` or try to use `POI.BOUNDS_AND_CONSTRAINTS` and leave it to ParametricOptInterface to interpret the constraints as bounds when applicable and linear constraints otherwise.
-
-## MOI Example - Parameters multiplying Quadratic terms
-
-Let's start with a simple quadratic problem
-
-```@example moi2
-import Ipopt
-import MathOptInterface as MOI
-import ParametricOptInterface as POI
-model = POI.Optimizer(Ipopt.Optimizer())
-x, _ = MOI.add_constrained_variable(model, MOI.GreaterThan(0.0))
-y, _ = MOI.add_constrained_variable(model, MOI.GreaterThan(0.0))
-ci1 = MOI.add_constraint(model, 2.0 * x + 1.0 * y, MOI.LessThan(4.0))
-ci2 = MOI.add_constraint(model, 1.0 * x + 2.0 * y, MOI.LessThan(4.0))
-MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
-obj_func = 1.0 * x * x + 1.0 * y * y
-MOI.set(model, MOI.ObjectiveFunction{typeof(obj_func)}(), obj_func)
-```
-
-To multiply a parameter in a quadratic term, the user will
-need to use the `POI.QuadraticObjectiveCoef` model attribute.
-
-```@example moi2
-p, _ = MOI.add_constrained_variable.(model, MOI.Parameter(1.0))
-MOI.set(model, POI.QuadraticObjectiveCoef(), (x, y), p)
-```
-
-This function will add the term `p*xy` to the objective function.
-It's also possible to multiply a scalar affine function to the quadratic term.
-
-```@example moi2
-MOI.set(model, POI.QuadraticObjectiveCoef(), (x, y), 2 * p + 3)
-```
-
-This will set the term `(2p+3)*xy` to the objective function (it overwrites the
-last set). Then, just optimize the model.
-
-```@example moi2
-MOI.optimize!(model)
-isapprox(MOI.get(model, MOI.ObjectiveValue()), 32/3, atol=1e-4)
-isapprox(MOI.get(model, MOI.VariablePrimal(), x), 4/3, atol=1e-4)
-isapprox(MOI.get(model, MOI.VariablePrimal(), y), 4/3, atol=1e-4)
-```
-
-To change the parameter just set `POI.ParameterValue` and optimize again.
-
-```@example moi2
-MOI.set(model, POI.ParameterValue(), p, 2.0)
-MOI.optimize!(model)
-isapprox(MOI.get(model, MOI.ObjectiveValue()), 128/9, atol=1e-4)
-isapprox(MOI.get(model, MOI.VariablePrimal(), x), 4/3, atol=1e-4)
-isapprox(MOI.get(model, MOI.VariablePrimal(), y), 4/3, atol=1e-4)
-```
-
-## JuMP Example - Parameters multiplying Quadratic terms
-
-Let's get the same MOI example
-
-```@example jump4
-using Ipopt
-using JuMP
-using ParametricOptInterface
-const POI = ParametricOptInterface
-
-optimizer = POI.Optimizer(Ipopt.Optimizer())
-model = direct_model(optimizer)
-
-@variable(model, x >= 0)
-@variable(model, y >= 0)
-@variable(model, p in MOI.Parameter(1.0))
-@constraint(model, 2x + y <= 4)
-@constraint(model, x + 2y <= 4)
-@objective(model, Max, (x^2 + y^2)/2)
-```
-
-We use the same MOI function to add the parameter multiplied to the quadratic term.
-
-```@example jump4
-MOI.set(backend(model), POI.QuadraticObjectiveCoef(), (index(x),index(y)), 2index(p)+3)
-```
-
-If the user print the `model`, the term `(2p+3)*xy` won't show.
-It's possible to retrieve the parametric function multiplying the term `xy` with `MOI.get`.
-
-```@example jump4
-MOI.get(backend(model), POI.QuadraticObjectiveCoef(), (index(x),index(y)))
-```
-
-Then, just optimize the model
-
-```@example jump4
-optimize!(model)
-isapprox(objective_value(model), 32/3, atol=1e-4)
-isapprox(value(x), 4/3, atol=1e-4)
-isapprox(value(y), 4/3, atol=1e-4)
-```
-
-To change the parameter just set `POI.ParameterValue` and optimize again.
-
-```@example jump4
-MOI.set(model, POI.ParameterValue(), p, 2.0)
-optimize!(model)
-isapprox(objective_value(model), 128/9, atol=1e-4)
-isapprox(value(x), 4/3, atol=1e-4)
-isapprox(value(y), 4/3, atol=1e-4)
-```
diff --git a/docs/src/Examples/markowitz.md b/docs/src/Examples/markowitz.md
deleted file mode 100644
index 6215c2a3..00000000
--- a/docs/src/Examples/markowitz.md
+++ /dev/null
@@ -1,100 +0,0 @@
-# Markowitz Efficient Frontier
-
-In this example, we solve the classical portfolio problem where we introduce the
-weight parameter $\gamma$ and maximize $\gamma \text{ risk} - \text{expected return}$.
-
-By updating the values of $\gamma$ we trace the efficient frontier.
-
-Given the prices changes with mean $\mu$ and covariance $\Sigma$, we can
-construct the classical portfolio problem:
-
-$$\begin{array}{ll}
-     \text{maximize}   & \gamma* x^T \mu - x^T \Sigma x \\
-     \text{subject to} & \| x \|_1 = 1 \\
-     & x \succeq 0
-\end{array}$$
-
-The problem data was gotten from the example [portfolio optimization](https://jump.dev/Convex.jl/dev/examples/portfolio_optimization/portfolio_optimization2/)
-
-```@repl markowitz
-using JuMP
-import Ipopt
-import ParametricOptInterface as POI
-import Plots
-μ = [11.5; 9.5; 6] / 100
-Σ = [166 34 58; 34 64 4; 58 4 100] / 100^2
-```
-
-We first build the model with $\gamma$ as parameter in POI
-```@repl markowitz
-function first_model(μ, Σ)
-    portfolio = Model() do
-        inner = MOI.instantiate(Ipopt.Optimizer; with_cache_type = Float64)
-        return POI.Optimizer(inner)
-    end
-    set_silent(portfolio)
-    N = length(μ)
-    @variable(portfolio, x[1:N] >= 0)
-    @variable(portfolio, γ in Parameter(0.0))
-    @objective(portfolio, Max, γ * μ' * x - x' * Σ * x)
-    @constraint(portfolio, sum(x) == 1)
-    optimize!(portfolio)
-    return portfolio
-end
-```
-
-Then, we update the $\gamma$ value in the model
-
-```@repl markowitz
-function update_model!(portfolio, γ_value)
-    γ = portfolio[:γ]
-    set_parameter_value(γ, γ_value)
-    optimize!(portfolio)
-    return portfolio
-end
-```
-
-Collecting all the return and risk resuls for each $\gamma$
-
-```@repl markowitz
-function add_to_dict(portfolios_values, portfolio, μ, Σ)
-    γ = portfolio[:γ]
-    γ_value = value(γ)
-    x = portfolio[:x]
-    x_value = value.(x)
-    portfolio_return = μ' * x_value
-    portfolio_deviation = x_value' * Σ * x_value
-    portfolios_values[γ_value] = (portfolio_return, portfolio_deviation)
-    return
-end
-```
-
-Run the portfolio optimization for different values of $\gamma$
-
-```@repl markowitz
-portfolio = first_model(μ, Σ)
-portfolios_values = Dict()
-# Create a reference to the model to change it later
-portfolio_ref = [portfolio]
-add_to_dict(portfolios_values, portfolio, μ, Σ)
-for γ_value in 0.02:0.02:1.0
-    portfolio_ref[] = update_model!(portfolio_ref[], γ_value)
-    add_to_dict(portfolios_values,portfolio_ref[], μ, Σ)
-end
-```
-
-Plot the efficient frontier
-
-```@repl markowitz
-sort!(portfolios_values, by = first)
-portfolios_values_matrix =
-    hcat([[v[1], v[2]] for v in values(portfolios_values)]...)'
-Plots.plot(
-    portfolios_values_matrix[:,2],
-    portfolios_values_matrix[:,1];
-    legend = false,
-    xlabel = "Standard Deviation",
-    ylabel = "Return",
-    title = "Efficient Frontier",
-)
-```
diff --git a/docs/src/index.md b/docs/src/index.md
index 8c08d66a..bbf3d486 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -1,14 +1,285 @@
-# ParametricOptInterface.jl Documentation
+# ParametricOptInterface.jl
 
-ParametricOptInterface.jl (POI for short) is a package written on top of MathOptInterface.jl that allows users to add parameters to a MOI/JuMP problem explicitly.
+[ParametricOptInterface.jl](https://github.com/jump-dev/ParametricOptInterface.jl)
+is a package that adds parameters to models in JuMP and MathOptInterface.
+
+## License
+
+`ParametricOptInterface.jl` is licensed under the
+[MIT License](https://github.com/jump-dev/ParametricOptInterface.jl/blob/master/LICENSE.md).
 
 ## Installation
 
-To install the package you can use `Pkg.add` as follows:
+Install ParametricOptInterface using `Pkg.add`:
+
 ```julia
-pkg> add ParametricOptInterface
+import Pkg
+Pkg.add("ParametricOptInterface")
+```
+
+## Use with JuMP
+
+Use ParametricOptInterface with JuMP by following this brief example:
+
+```@repl
+using JuMP, HiGHS
+import ParametricOptInterface as POI
+model = Model(() -> POI.Optimizer(HiGHS.Optimizer()))
+@variable(model, x)
+@variable(model, p in Parameter(1))
+@constraint(model, p * x + p >= 3)
+@objective(model, Min, 2x + p)
+optimize!(model)
+value(x)
+set_parameter_value(p, 2.0)
+optimize!(model)
+value(x)
 ```
 
-## Contributing
+## How and why ParametricOptInterface is needed
+
+JuMP and MathOptInterface have support for _parameters_. Parameters are decision
+variables that belong to the `Parameter` set. The `Parameter` set is
+conceptually similar to the `EqualTo` set, except that solvers may treat a
+decision variable constrained to the `Parameter` set as a constant, and they
+need not add it as a decision variable to the model.
+
+In JuMP, a parameter can be added using the following syntax:
+```@repl
+using JuMP
+model = Model();
+@variable(model, p in Parameter(2))
+parameter_value(p)
+set_parameter_value(p, 3.0)
+parameter_value(p)
+```
+
+In MathOptInterface, a parameter can be added using the following syntax:
+```@repl
+import MathOptInterface as MOI
+model = MOI.Utilities.Model{Float64}();
+p, p_con = MOI.add_constrained_variable(model, MOI.Parameter(2.0))
+MOI.get(model, MOI.ConstraintSet(), p_con)
+new_set = MOI.Parameter(3.0)
+MOI.set(model, MOI.ConstraintSet(), p_con, new_set)
+MOI.get(model, MOI.ConstraintSet(), p_con)
+```
+
+### Some solvers have native support for parameters
+
+Some solvers have native support for parameters. example is Ipopt. To
+demonstrate, look at the following example. Even though there are two
+`@variable` calls, the log of Ipopt shows that it solved a problem with only one
+decision variable:
+```@repl
+using JuMP, Ipopt
+model = Model(Ipopt.Optimizer)
+@variable(model, x)
+@variable(model, p in Parameter(1))
+@constraint(model, x + p >= 3)
+@objective(model, Min, 2x)
+optimize!(model)
+```
+Internally, Ipopt replaced the parameter `p` with the constant `1.0`, and solved
+the problem:
+```@repl
+using JuMP, Ipopt
+model = Model(Ipopt.Optimizer)
+@variable(model, x)
+@constraint(model, x + 1 >= 3)
+@objective(model, Min, 2x)
+optimize!(model)
+```
+
+Parameters are most useful when you want to solve a sequence of problems in
+which some of the data changes between iterations:
+```@repl
+using JuMP, Ipopt
+model = Model(Ipopt.Optimizer)
+set_silent(model)
+@variable(model, x)
+@variable(model, p in Parameter(1))
+@constraint(model, x + p >= 3)
+@objective(model, Min, 2x)
+solution = Dict{Int,Float64}()
+for p_value in 0:5
+    set_parameter_value(p, p_value)
+    optimize!(model)
+    assert_is_solved_and_feasible(model)
+    solution[p_value] = value(x)
+end
+solution
+```
+
+### Some solvers do not have native support for parameters
+
+Even though solvers like Ipopt support parameters, many solvers do not. One
+example is HiGHS. Despite the fact that HiGHS doesn't support parameters, you
+can still build and solve a model with parameters:
+```@repl index_highs
+using JuMP, HiGHS
+model = Model(HiGHS.Optimizer)
+@variable(model, x)
+@variable(model, p in Parameter(1))
+@constraint(model, x + p >= 3)
+@objective(model, Min, 2x)
+optimize!(model)
+```
+This works because, behind the scenes, the bridges in MathOptInterface rewrote
+`p in Parameter(1)` to `p in MOI.EqualTo(1.0)`:
+```@repl index_highs
+print_active_bridges(model)
+```
+Thus, HiGHS solved the problem:
+```@repl
+using JuMP, HiGHS
+model = Model(HiGHS.Optimizer)
+@variable(model, x)
+@variable(model, p == 1)
+@constraint(model, x + p >= 3)
+@objective(model, Min, 2x)
+optimize!(model)
+```
+
+The downside to the bridge approach is that it adds a new decision variable with
+fixed bounds for every parameter in the problem. Moreover, the bridge approach
+does not handle `parameter * variable` terms, because the resulting problem is a
+quadratic constraint:
+
+```jldoctest
+julia> using JuMP, HiGHS
+
+julia> model = Model(HiGHS.Optimizer);
+
+julia> @variable(model, x);
 
-When contributing please note that the package follows the [JuMP style guide](https://jump.dev/JuMP.jl/stable/developers/style/).
+julia> @variable(model, p in Parameter(1));
+
+julia> @constraint(model, p * x >= 3)
+ERROR: Constraints of type MathOptInterface.ScalarQuadraticFunction{Float64}-in-MathOptInterface.GreaterThan{Float64} are not supported by the solver.
+
+If you expected the solver to support your problem, you may have an error in your formulation. Otherwise, consider using a different solver.
+
+The list of available solvers, along with the problem types they support, is available at https://jump.dev/JuMP.jl/stable/installation/#Supported-solvers.
+Stacktrace:
+[...]
+```
+
+### ParametricOptInterface
+
+ParametricOptInterface provides [`Optimizer`](@ref), which is a meta-optimizer
+that wraps another optimizer. Instead of adding fixed variables to the model,
+POI substitutes out the parameters with their value before passing the
+constraint or objective to the inner optimizer. When the parameter value is
+changed, POI efficiently modifies the inner optimizer to reflect the new
+parameter values.
+
+```@repl
+using JuMP, HiGHS
+import ParametricOptInterface as POI
+model = Model(() -> POI.Optimizer(HiGHS.Optimizer()));
+@variable(model, x);
+@variable(model, p in Parameter(1));
+@constraint(model, x + p >= 3);
+@objective(model, Min, 2x);
+optimize!(model)
+```
+Note how HiGHS now solves a problem with one decision variable.
+
+Because POI replaces parameters with their constant value, POI supports
+`parameter * variable` terms:
+```@repl
+using JuMP, HiGHS
+import ParametricOptInterface as POI
+model = Model(() -> POI.Optimizer(HiGHS.Optimizer()));
+@variable(model, x);
+@variable(model, p in Parameter(1));
+@constraint(model, p * x >= 3)
+@objective(model, Min, 2x)
+optimize!(model)
+```
+
+## When to use ParametricOptInterface
+
+You should use ParametricOptInterface when:
+
+ * you are using a solver that does not have native support for parameters
+ * you are solving a single problem for multiple values of the parameters.
+
+For problems with a small number of parameters, and in which the parameters
+appear additively in the constraints and the objective, it may be more efficient
+to use the bridge approach. In general, you should try with and without POI and
+choose the approach which works best for your model.
+
+## The dual of a parameter
+
+In some applications you may need the dual of a parameter. The dual can be
+computed only if the parameter appears additively in the problem. Query the dual
+assocaited with the parameter using [`ParameterDual`](@ref):
+```@repl
+using JuMP, HiGHS
+import ParametricOptInterface as POI
+model = Model(() -> POI.Optimizer(HiGHS.Optimizer()));
+set_silent(model)
+@variable(model, x);
+@variable(model, p in Parameter(1));
+@constraint(model, x + p >= 3);
+@objective(model, Min, 2x);
+optimize!(model)
+get_attribute(p, POI.ParameterDual())
+```
+
+Note how the dual is the same as the `reduced_cost` of an equivalent fixed variable:
+```@repl
+using JuMP, HiGHS
+model = Model(HiGHS.Optimizer);
+set_silent(model)
+@variable(model, x);
+@variable(model, p == 1);
+@constraint(model, x + p >= 3);
+@objective(model, Min, 2x);
+optimize!(model)
+reduced_cost(p)
+```
+
+## Variable bounds
+
+An ambiguity arises when the user writes a model like:
+```@repl
+using JuMP
+model = Model();
+@variable(model, x)
+@variable(model, p in Parameter(2))
+@constraint(model, x >= p)
+```
+Did they mean an affine constraint like `1.0 * x + 0.0 in GreaterThan(2.0)`, or
+did they mean a variable bound like `x in GreaterThan(2.0)`?
+
+By default, ParametricOptInterface does not attempt to simplify affine
+constraints involving parameters to variable bounds, but this behavior can be
+controlled using the [`ConstraintsInterpretation`](@ref) attribute.
+
+For example, the default is:
+```@repl
+using JuMP, HiGHS
+import ParametricOptInterface as POI
+model = Model(() -> POI.Optimizer(HiGHS.Optimizer()));
+set_attribute(model, POI.ConstraintsInterpretation(), POI.ONLY_CONSTRAINTS)
+@variable(model, x);
+@variable(model, p in Parameter(1));
+@constraint(model, x >= p)
+optimize!(model)
+```
+but by setting the [`ConstraintsInterpretation`](@ref) attribute to
+`POI.BOUNDS_AND_CONSTRAINTS`, we can solve a problem with one decision variable
+and zero constraint rows:
+```@repl
+using JuMP, HiGHS
+import ParametricOptInterface as POI
+model = Model(() -> POI.Optimizer(HiGHS.Optimizer()));
+set_attribute(model, POI.ConstraintsInterpretation(), POI.BOUNDS_AND_CONSTRAINTS)
+@variable(model, x);
+@variable(model, p in Parameter(1));
+@constraint(model, x >= p)
+optimize!(model)
+```
diff --git a/docs/src/manual.md b/docs/src/manual.md
deleted file mode 100644
index 8a00c73a..00000000
--- a/docs/src/manual.md
+++ /dev/null
@@ -1,114 +0,0 @@
-```@meta
-CurrentModule = ParametricOptInterface
-```
-
-# Manual
-
-## Why use parameters?
-
-A typical optimization model built using `MathOptInterface.jl` (`MOI`for short)
-has two main components:
-
-1. Variables
-2. Constants
-
-Using these basic elements, one can create functions and sets that, together,
-form the desired optimization model. The goal of `POI` is the implementation of
-a third type, parameters, which:
-
-* are declared similar to a variable, and inherits some functionalities (for
-  example, dual calculation)
-* acts like a constant, in the sense that it has a fixed value that will remain
-  the same unless explicitely changed by the user
-
-A main concern is to efficiently implement this new type, as one typical usage
-is to change its value to analyze the model behavior, without the need to build
-a new one from scratch.
-
-## How it works
-
-The main idea applied in POI is that the interaction between the solver, for
-example `HiGHS`, and the optimization model will be handled by `MOI` as usual.
-
-Because of that, `POI` is a higher level wrapper around `MOI`, responsible for
-receiving variables, constants and parameters, and forwarding to the lower level
-model only variables and constants.
-
-As `POI` receives parameters, it must analyze and decide how they should be
-handled on the lower level optimization model (the `MOI` model).
-
-## Usage
-
-In this manual we describe how to interact with the optimization model at the
-MOI level. In the _Examples_ section you can find some tutorials with the
-JuMP usage.
-
-### Supported constraints
-
-This is a list of supported `MOI` constraint functions that can handle
-parameters. If you try to add a parameter to a function that is not listed here,
-it will return an unsupported error.
-
-| MOI Function              |
-| :------------------------ |
-| `ScalarAffineFunction`    |
-| `ScalarQuadraticFunction` |
-| `VectorAffineFunction`    |
-
-
-### Supported objective functions
-
-| MOI Function              |
-| :------------------------ |
-| `ScalarAffineFunction`    |
-| `ScalarQuadraticFunction` |
-
-### Declare a Optimizer
-
-In order to use parameters, the user needs to declare an
-[`Optimizer`](@ref) on top of a `MOI` optimizer, such as `HiGHS.Optimizer()`.
-
-```@repl manual
-import ParametricOptInterface as POI
-import HiGHS
-optimizer = POI.Optimizer(HiGHS.Optimizer())
-```
-
-### Parameters
-
-A `MOI.Parameter` is a set used to define a variable with a fixed value that
-can be changed by the user. It is analogous to `MOI.EqualTo`, but can be used
-by special methods like the ones in this package to remove the fixed variable
-from the optimization problem. This permits the usage of multiplicative
-parameters in linear models and might speedup solves since the number of
-variables is reduced.
-
-### Adding a new parameter to a model
-
-To add a parameter to a model, we must use the `MOI.add_constrained_variable`
-function, passing as its arguments the model and a `MOI.Parameter` with its
-given value:
-
-```@repl manual
-y, cy = MOI.add_constrained_variable(optimizer, MOI.Parameter(0.0))
-```
-
-### Changing the parameter value
-
-To change a given parameter's value, access its `VariableIndex` and set it to
-the new value using the `MOI.Parameter` structure.
-
-```@repl manual
-MOI.set(optimizer, POI.ParameterValue(), y, MOI.Parameter(2.0))
-```
-
-### Retrieving the dual of a parameter
-
-Given an optimized model, one can compute the dual associated to a parameter,
-**as long as it is an additive term in the constraints or objective**.
-
-One can do so by getting the `MOI.ConstraintDual` attribute of the parameter's
-`MOI.ConstraintIndex`:
-```julia
-julia> MOI.get(optimizer, POI.ParameterDual(), y)
-```
diff --git a/src/MOI_wrapper.jl b/src/MOI_wrapper.jl
index 70b52a3d..6352876d 100644
--- a/src/MOI_wrapper.jl
+++ b/src/MOI_wrapper.jl
@@ -1908,24 +1908,42 @@ end
 
 Attribute to define how [`Optimizer`](@ref) should interpret constraints.
 
-- `POI.ONLY_CONSTRAINTS`: Only interpret `ScalarAffineFunction` constraints as linear constraints
-  If an expression such as `x >= p1 + p2` appears it will be trated like a new constraint.
+- `POI.ONLY_CONSTRAINTS`: always interpret `ScalarAffineFunction` constraints as
+  linear constraints. If an expression such as `x >= p1 + p2` appears, it will
+  be treated like an affine constraint.
   **This is the default behaviour of [`Optimizer`](@ref)**
 
-- `POI.ONLY_BOUNDS`: Only interpret `ScalarAffineFunction` constraints as a variable bound.
-  This is valid for constraints such as `x >= p` or `x >= p1 + p2`. If a constraint `x1 + x2 >= p` appears,
-  which is not a valid variable bound it will throw an error.
+- `POI.ONLY_BOUNDS`: always interpret `ScalarAffineFunction` constraints as a
+  variable bound. This is valid for constraints such as `x >= p` or
+  `x >= p1 + p2`. If a constraint `x1 + x2 >= p` appears which is not a valid
+  variable bound, an error will be thrown.
 
-- `POI.BOUNDS_AND_CONSTRAINTS`: Interpret `ScalarAffineFunction` constraints as a variable bound if they
-  are a valid variable bound, i.e., `x >= p` or `x >= p1 + p2` and interpret them as linear constraints
-  otherwise.
+- `POI.BOUNDS_AND_CONSTRAINTS`: interpret `ScalarAffineFunction` constraints as
+  a variable bound if they are a valid variable bound, for example, `x >= p` or
+  `x >= p1 + p2`, and interpret them as linear constraints otherwise.
 
 # Example
 
-```julia
-MOI.set(model, POI.ConstraintsInterpretation(), POI.ONLY_BOUNDS)
-MOI.set(model, POI.ConstraintsInterpretation(), POI.ONLY_CONSTRAINTS)
-MOI.set(model, POI.ConstraintsInterpretation(), POI.BOUNDS_AND_CONSTRAINTS)
+```jldoctest
+julia> import MathOptInterface as MOI
+
+julia> import ParametricOptInterface as POI
+
+julia> model = POI.Optimizer(MOI.Utilities.Model{Float64}())
+ParametricOptInterface.Optimizer{Float64, MOIU.Model{Float64}}
+├ ObjectiveSense: FEASIBILITY_SENSE
+├ ObjectiveFunctionType: MOI.ScalarAffineFunction{Float64}
+├ NumberOfVariables: 0
+└ NumberOfConstraints: 0
+
+julia> MOI.set(model, POI.ConstraintsInterpretation(), POI.ONLY_BOUNDS)
+ONLY_BOUNDS::ConstraintsInterpretationCode = 1
+
+julia> MOI.set(model, POI.ConstraintsInterpretation(), POI.ONLY_CONSTRAINTS)
+ONLY_CONSTRAINTS::ConstraintsInterpretationCode = 0
+
+julia> MOI.set(model, POI.ConstraintsInterpretation(), POI.BOUNDS_AND_CONSTRAINTS)
+BOUNDS_AND_CONSTRAINTS::ConstraintsInterpretationCode = 2
 ```
 """
 struct ConstraintsInterpretation <: MOI.AbstractOptimizerAttribute end