diff --git a/docs/Project.toml b/docs/Project.toml
index b10a7e207..c3f94d7a6 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,4 +1,5 @@
 [deps]
+Alpine = "07493b3f-dabb-5b16-a503-4139292d7dd4"
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 CSDP = "0a46da34-8e4b-519e-b418-48813639ff34"
 Clarabel = "61c947e1-3e6d-4ee4-985a-eec8c727bd6e"
@@ -10,6 +11,7 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 Dualization = "191a621a-6537-11e9-281d-650236a99e60"
 DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07"
 GroupsCore = "d5909c97-4eac-4ecc-a3dc-fdd0858a4120"
+HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b"
 HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327"
 ImplicitPlots = "55ecb840-b828-11e9-1645-43f4a9f9ace7"
 Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
@@ -22,6 +24,7 @@ MultivariateBases = "be282fd4-ad43-11e9-1d11-8bd9d7e43378"
 MultivariateMoments = "f4abf1af-0426-5881-a0da-e2f168889b5e"
 MultivariatePolynomials = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3"
 MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
+Pavito = "cd433a01-47d1-575d-afb7-6db927ee8d8f"
 PermutationGroups = "8bc5a954-2dfc-11e9-10e6-cd969bffa420"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 PolyJuMP = "ddf597a6-d67e-5340-b84c-e37d84115374"
diff --git a/docs/src/tutorials/Polynomial Optimization/polynomial_optimization.jl b/docs/src/tutorials/Polynomial Optimization/polynomial_optimization.jl
index a4af265a8..c2f946910 100644
--- a/docs/src/tutorials/Polynomial Optimization/polynomial_optimization.jl	
+++ b/docs/src/tutorials/Polynomial Optimization/polynomial_optimization.jl	
@@ -35,7 +35,7 @@ model = Model(Ipopt.Optimizer)
 @variable(model, a >= 0)
 @variable(model, b >= 0)
 @constraint(model, a + b >= 1)
-@NLobjective(model, Min, a^3 - a^2 + 2a*b - b^2 + b^3)
+@objective(model, Min, a^3 - a^2 + 2a*b - b^2 + b^3)
 optimize!(model)
 
 # As we can see below, the termination status is `LOCALLY_SOLVED` and not of `OPTIMAL`
@@ -93,6 +93,37 @@ solution_summary(gmodel)
 @test value(b) ≈ 0.5 rtol=1e-5 #src
 value(a), value(b)
 
+# ## QCQP approach
+
+import Alpine, HiGHS, Ipopt, Pavito
+ipopt = optimizer_with_attributes(
+    Ipopt.Optimizer,
+    MOI.Silent() => true,
+)
+highs = optimizer_with_attributes(
+    HiGHS.Optimizer,
+    "presolve" => "on",
+    "log_to_console" => false,
+)
+pavito = optimizer_with_attributes(
+    Pavito.Optimizer,
+    MOI.Silent() => true,
+    "mip_solver" => highs,
+    "cont_solver" => ipopt,
+    "mip_solver_drives" => false,
+)
+alpine = optimizer_with_attributes(
+    Alpine.Optimizer,
+    "nlp_solver" => ipopt,
+    "mip_solver" => pavito,
+)
+set_optimizer(model, () -> PolyJuMP.QCQP.Optimizer(MOI.instantiate(alpine)))
+optimize!(model)
+
+# We can see that it found the optimal solution
+
+termination_status(model), value(a), value(b)
+
 # ## Sum-of-Squares approach
 
 # We will now see how to find the optimal solution using Sum of Squares Programming.