diff --git a/CHANGELOG.md b/CHANGELOG.md
index 0bce43e..81dc371 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,11 +1,21 @@
 QuantumCircuitOpt.jl Change Log
 ===============================
 
+### v0.5.0
+- Minor: Fixed result `primal_status` issue in `log.jl`
+- Added helper functions for obtaining fixed indices of `U_var` unitary variables to zeros or constant values. Dropped linearization constraints for fixed `U_var` variables
+- Dropped support for `unit_magnitude_constraints`
+- Reformulated quadratic objective function of `slack_var` variables into a linear outer-approximation - increases code coverage due to MILP reformulation 
+- Fixed `slack_var` based on fixed `U_var` variables  
+- Minor updates in error messages for kron symboled gates
+- Changed Cbc test dependency to HiGHS (MIP) solver 
+- Updated docs and unit tests to reflect above changes
+
 ### v0.4.1
 - Added support for returning exact feasible decomposition (without optimality certificate) using `exact_feasible`
 - Added support for unitary constraints using binary linearizations (speeds up feasiblity problems)
-- Added `decompose_W_using_HSTCnot` in examples
-- Added `decompose_CNOT_using_GR` in examples (for Rydberg atom array-based simulators)
+- Added `W_using_HSTCnot` in examples
+- Added `CNOT_using_GR` in examples (for Rydberg atom array-based simulators)
 - Added support for `CSGate`, `CSdaggerGate`, `CTGate`, `CTdaggerGate` and `SSwapGate` (sqrt(`SwapGate`))
 - Dropped support for `CZRevGate` as it is invariant to qubit flip
 - Updated README and docs with the Youtube link for the JuliaCon's presentation
diff --git a/Project.toml b/Project.toml
index 3da2d25..ede1f12 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "QuantumCircuitOpt"
 uuid = "88128e30-b60a-4e54-ab02-1050a5f92a36"
 authors = ["Harsha Nagarajan"]
-version = "0.4.1"
+version = "0.5.0"
 
 [deps]
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
@@ -16,8 +16,8 @@ Memento = "~1.0, ~1.1, 1"
 julia = "^1"
 
 [extras]
-Cbc = "9961bab8-2fa3-5c5a-9d89-47fab24efd76"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b"
 
 [targets]
-test = ["Cbc", "Test"]
+test = ["HiGHS", "Test"]
diff --git a/README.md b/README.md
index 25ba7ef..7b87fa7 100644
--- a/README.md
+++ b/README.md
@@ -5,6 +5,7 @@
 Status: 
 [![CI](https://github.com/harshangrjn/QuantumCircuitOpt.jl/actions/workflows/ci.yml/badge.svg)](https://github.com/harshangrjn/QuantumCircuitOpt.jl/actions/workflows/ci.yml)
 [![codecov](https://codecov.io/gh/harshangrjn/QuantumCircuitOpt.jl/branch/master/graph/badge.svg?token=KGJWIV6QF4)](https://codecov.io/gh/harshangrjn/QuantumCircuitOpt.jl)
+[![version][https://juliahub.com/docs/QuantumCircuitOpt/version.svg]][https://juliahub.com/ui/Packages/QuantumCircuitOpt/dwSy1]
 
 <!-- Stable version: [![Documentation](https://github.com/harshangrjn/QuantumCircuitOpt.jl/actions/workflows/documentation.yml/badge.svg)](https://harshangrjn.github.io/QuantumCircuitOpt.jl/stable/) -->
 Dev version: [![Documentation](https://github.com/harshangrjn/QuantumCircuitOpt.jl/actions/workflows/documentation.yml/badge.svg)](https://harshangrjn.github.io/QuantumCircuitOpt.jl/dev/)
diff --git a/docs/src/index.md b/docs/src/index.md
index 55e8871..bcbec35 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -22,11 +22,11 @@ import Pkg
 Pkg.add("QuantumCircuitOpt")
 ```
 
-At least one mixed-integer programming solver is required for running QuantumCircuitOpt. The well-known [CPLEX](https://github.com/jump-dev/CPLEX.jl) or the [Gurobi](https://github.com/jump-dev/Gurobi.jl) solver is highly recommended, as it is fast, scaleable and can be used to solve on fairly large-scale graphs. However, open-source solvers such as [Cbc](https://github.com/jump-dev/Cbc.jl) or [GLPK](https://github.com/jump-dev/GLPK.jl) is also compatible with QuantumCircuitOpt which can be installed via the package manager with
+At least one mixed-integer programming (MIP) solver is required for running QuantumCircuitOpt. The well-known [Gurobi](https://github.com/jump-dev/Gurobi.jl) or IBM's [CPLEX](https://github.com/jump-dev/CPLEX.jl) solver is highly recommended, as it is fast, scaleable and can be used to solve on fairly large-scale circuits. However, the open-source MIP solver [HiGHS](https://github.com/jump-dev/HiGHS.jl) is also compatible with QuantumCircuitOpt. Gurobi (or any other MIP solver) can be installed via the package manager with
 
 ```julia
 import Pkg
-Pkg.add("Cbc")
+Pkg.add("Gurobi")
 ```
 
 ## Unit Tests
diff --git a/docs/src/quickguide.md b/docs/src/quickguide.md
index f00a8e3..ed1a045 100644
--- a/docs/src/quickguide.md
+++ b/docs/src/quickguide.md
@@ -34,7 +34,7 @@ Building on the recent success of [Julia](https://julialang.org), [JuMP](https:/
 ```
 
 ## Getting started
-To get started, install [QCOpt](https://github.com/harshangrjn/QuantumCircuitOpt.jl) and [JuMP](https://github.com/jump-dev/JuMP.jl), a modeling language layer for optimization. QCOpt also needs a MIP solver such as [CPLEX](https://github.com/jump-dev/CPLEX.jl) or [Gurobi](https://github.com/jump-dev/Gurobi.jl). If you prefer an open-source MIP solver, install [CBC](https://github.com/jump-dev/Cbc.jl) or [GLPK](https://github.com/jump-dev/GLPK.jl) from the Julia package manager, though be warned that the run times of QCOpt can be substantially slower using these open-source MIP solvers. 
+To get started, install [QCOpt](https://github.com/harshangrjn/QuantumCircuitOpt.jl) and [JuMP](https://github.com/jump-dev/JuMP.jl), a modeling language layer for optimization. QCOpt also needs a MIP solver such as [Gurobi](https://github.com/jump-dev/Gurobi.jl) or IBM's [CPLEX](https://github.com/jump-dev/CPLEX.jl). If you prefer an open-source MIP solver, install [HiGHS](https://github.com/jump-dev/HiGHS.jl) from the Julia package manager, though be warned that the run times of QCOpt can be substantially slower using any of the open-source MIP solvers. 
 
 # User inputs
 QCOpt takes two types of user-defined input specifications. The first type of input contains all the necessary circuit specifications. This is given by a dictionary in Julia, which is a collection of key-value pairs, where every key is of the type `String`, which admits values of various types. Below is the list of allowable keys for the dictionary, given in column 1, and it's respective values with descriptions, given in column 2. This input dictionary is represented as `params` in all the [example](https://github.com/harshangrjn/QuantumCircuitOpt.jl/tree/master/examples) circuit decompositions. 
@@ -88,7 +88,7 @@ The second set of inputs for QCOpt contains all the optional specifications for
 |`identity_gate_symmetry_constraints`| This option activates the valid constraints to eliminate symmetry in the Identity gate in the decomposition (default: `true`)| 
 |`idempotent_gate_constraints`| This option activates the valid constraints to eliminate idempotent gates in the elementary (native) gates set (default: `false`)| 
 |`convex_hull_gate_constraints`| This option activates the valid constraints to apply convex hull of complex entries in the elementary (native) gates set (default: `false`)| 
-|`unit_magnitude_constraints`| This option activates the valid outer-approximation constraints for unit-valued complex entries of the unitary gates (`U_var`) (default: `false`)|
+|`fix_unitary_variables`| This option evaluates all the fixed-valued indices of unitary matrix varaibles (`U_var`) at every depth, and appropriately builds the optimization model (default: `true`)|
 |`visualize_solution`| This option activates the visualization of the optimal circuit decomposition (default: `true`)| 
 | `relax_integrality` | This option transforms integer variables into continuous variables (default: `false`).  |
 | `optimizer_log` | This option enables or disables console logging for the `optimizer` (default: `true`).|
@@ -161,23 +161,25 @@ QCOpt.visualize_solution(results, data)
 For example, for the above controlled-Z gate decomposition, the processed output of QCOpt is as follows: 
 ```
 =============================================================================
+QuantumCircuitOpt version: v0.5.0
+
 Quantum Circuit Model Data
 
   Number of qubits: 2
   Total number of elementary gates (after presolve): 72
   Maximum depth of decomposition: 4
-  Input elementary gates: ["U3_1", "U3_2", "CNot_1_2", "Identity"]
+  Elementary gates: ["U3_1", "U3_2", "CNot_1_2", "Identity"]
     U3_θ discretization: [-180.0, -90.0, 0.0, 90.0, 180.0]
     U3_ϕ discretization: [-180.0, -90.0, 0.0, 90.0, 180.0]
     U3_λ discretization: [-180.0, -90.0, 0.0, 90.0, 180.0]
-  Type of decomposition: exact
+  Type of decomposition: exact_optimal
   MIP optimizer: Gurobi
 
 Optimal Circuit Decomposition
 
   U3_2(-90.0,0.0,0.0) * CNot_1_2 * U3_2(90.0,0.0,0.0) = Target gate
   Minimum optimal depth: 3
-  Optimizer run time: 2.64 sec.
+  Optimizer run time: 3.01 sec.
 =============================================================================
 ```
 
diff --git a/examples/2qubit_gates.jl b/examples/2qubit_gates.jl
index 8045d6e..dcb8d18 100644
--- a/examples/2qubit_gates.jl
+++ b/examples/2qubit_gates.jl
@@ -1,4 +1,4 @@
-function decompose_hadamard()
+function hadamard()
 
     println(">>>>> Hadamard Gate <<<<<")
  
@@ -17,7 +17,7 @@ function decompose_hadamard()
     )
 end
 
-function decompose_controlled_Z()
+function controlled_Z()
 
     println(">>>>> Controlled-Z Gate <<<<<")
 
@@ -37,7 +37,7 @@ function decompose_controlled_Z()
 
 end
 
-function decompose_controlled_V()
+function controlled_V()
 
     println(">>>>> Controlled-V Gate <<<<<")
 
@@ -52,7 +52,7 @@ function decompose_controlled_V()
     )
 end
 
-function decompose_controlled_H()
+function controlled_H()
 
     println(">>>>> Controlled-H Gate <<<<<")
 
@@ -71,7 +71,7 @@ function decompose_controlled_H()
     )
 end
 
-function decompose_controlled_H_with_R()
+function controlled_H_with_R()
 
     println(">>>>> Controlled-H with R Gate <<<<<")
 
@@ -90,7 +90,7 @@ function decompose_controlled_H_with_R()
     )
 end
 
-function decompose_magic()
+function magic()
     
     println(">>>>> M Gate <<<<<")
 
@@ -109,7 +109,7 @@ function decompose_magic()
         )
 end
 
-function decompose_magic_using_CNOT_1_2()
+function magic_using_CNOT_1_2()
     
     println(">>>>> M Gate <<<<<")
 
@@ -128,7 +128,7 @@ function decompose_magic_using_CNOT_1_2()
         )
 end
 
-function decompose_magic_using_SHCnot()
+function magic_using_SHCnot()
     
     println(">>>>> M gate using S, H and CNOT Gate <<<<<")
 
@@ -143,7 +143,7 @@ function decompose_magic_using_SHCnot()
         )
 end
 
-function decompose_S()
+function S_gate()
 
     println(">>>>> S Gate <<<<<")
 
@@ -162,7 +162,7 @@ function decompose_S()
         )
 end
 
-function decompose_revcnot()
+function revcnot()
 
     println(">>>>> Reverse CNOT Gate <<<<<")
 
@@ -173,11 +173,11 @@ function decompose_revcnot()
     "elementary_gates" => ["H_1", "H_2", "CNot_1_2", "Identity"],  
     "target_gate" => QCOpt.CNotRevGate(),
     "objective" => "minimize_depth", 
-    "decomposition_type" => "exact_optimal"     
+    "decomposition_type" => "exact_optimal"
     )
 end
 
-function decompose_revcnot_with_U()
+function revcnot_with_U()
     
     println(">>>>> Reverse CNOT using U3 and CNOT Gates <<<<<")
 
@@ -196,7 +196,7 @@ function decompose_revcnot_with_U()
         )
 end
 
-function decompose_swap()
+function swap()
 
     println(">>>>> SWAP Gate <<<<<")
 
@@ -211,7 +211,7 @@ function decompose_swap()
         )
 end
 
-function decompose_W()
+function W_gate()
 
     println(">>>>> W Gate <<<<<")
 
@@ -230,7 +230,7 @@ function decompose_W()
         )
 end
 
-function decompose_W_using_HSTCnot()
+function W_using_HSTCnot()
 
     println(">>>>> W Gate using H, S, T and CNOT gates <<<<<")
 
@@ -245,7 +245,7 @@ function decompose_W_using_HSTCnot()
     )
 end
 
-function decompose_W_using_HCnot()
+function W_using_HCnot()
 
     println(">>>>> W Gate using H and CNOT gates <<<<<")
 
@@ -260,7 +260,7 @@ function decompose_W_using_HCnot()
         )
 end
 
-function decompose_HCoinGate()
+function HCoinGate()
 
     println(">>>>> Hadamard Coin gate <<<<<")
 
@@ -275,7 +275,7 @@ function decompose_HCoinGate()
         )
 end
 
-function decompose_GroverDiffusion_using_HX()
+function GroverDiffusion_using_HX()
 
     println(">>>>> Grover's Diffusion Operator <<<<<")
 
@@ -292,7 +292,7 @@ function decompose_GroverDiffusion_using_HX()
         )
 end
 
-function decompose_GroverDiffusion_using_Clifford()
+function GroverDiffusion_using_Clifford()
 
     println(">>>>> Grover's Diffusion Operator <<<<<")
 
@@ -308,7 +308,7 @@ function decompose_GroverDiffusion_using_Clifford()
 end
 
 
-function decompose_GroverDiffusion_using_U3()
+function GroverDiffusion_using_U3()
 
     println(">>>>> Grover's Diffusion Operator using U3 gate <<<<<")
 
@@ -327,7 +327,7 @@ function decompose_GroverDiffusion_using_U3()
     )
 end
 
-function decompose_iSwap()
+function iSwap()
 
     println(">>>>> iSwap Gate <<<<<")
 
@@ -339,11 +339,11 @@ function decompose_iSwap()
         # "elementary_gates" => ["S_1", "S_2", "Sdagger_1", "Sdagger_2", "H_1", "H_2", "CNot_1_2", "CNot_2_1", "Identity"],
         "target_gate" => QCOpt.iSwapGate(),
         "objective" => "minimize_depth",
-        "decomposition_type" => "exact_optimal",                           
+        "decomposition_type" => "exact_optimal",  
         )
 end
 
-function decompose_qft2_using_R()
+function qft2_using_R()
     println(">>>>> QFT2 Gate using R and CNOT gates <<<<<")
 
     return Dict{String, Any}(
@@ -353,14 +353,14 @@ function decompose_qft2_using_R()
         "elementary_gates" => ["RX_1", "RZ_1", "RZ_2", "CNot_1_2", "CNot_2_1", "Identity"],
         "target_gate" => QCOpt.QFT2Gate(),
         "objective" => "minimize_depth",
-        "decomposition_type" => "approximate",
+        "decomposition_type" => "exact_feasible",
         
         "RX_discretization" => [π/2],
         "RZ_discretization" => [-π/4, π/2, 3*π/4, 7*π/4],                      
         )
 end
 
-function decompose_qft2_using_HT()
+function qft2_using_HT()
     println(">>>>> QFT2 Gate using H, T, CNOT gates <<<<<")
 
     return Dict{String, Any}(
@@ -374,11 +374,11 @@ function decompose_qft2_using_HT()
         )
 end
 
-function decompose_GRGate()
+function GRGate()
     println(">>>>> GRGate testing <<<<<")
 
     GR1 = QCOpt.get_full_sized_gate("GR", 2; angle = [π/6, π/3])
-    GR2 = QCOpt.get_full_sized_gate("GR", 2; angle = [π/3, π/6])
+    # GR2 = QCOpt.get_full_sized_gate("GR", 2; angle = [π/3, π/6])
     CNot_1_2 = QCOpt.get_full_sized_gate("CNot_1_2", 2)
     T = QCOpt.round_complex_values(GR1 * CNot_1_2)
 
@@ -395,7 +395,7 @@ function decompose_GRGate()
         )
 end
 
-function decompose_X_using_GR()
+function X_using_GR()
 
     println(">>>>> Pauli-X Gate using global rotation <<<<<")
  
@@ -415,7 +415,7 @@ function decompose_X_using_GR()
     )
 end
 
-function decompose_CNOT_using_GR()
+function CNOT_using_GR()
     println(">>>>> CNOT Gate using global rotation <<<<<")
     
     return Dict{String, Any}(
diff --git a/examples/3qubit_gates.jl b/examples/3qubit_gates.jl
index d12fcd1..9d3c029 100644
--- a/examples/3qubit_gates.jl
+++ b/examples/3qubit_gates.jl
@@ -1,4 +1,4 @@
-function decompose_RX_on_q3()
+function RX_on_q3()
 
     println(">>>>> RX Gate on third qubit using U3Gate <<<<<")
  
@@ -17,7 +17,7 @@ function decompose_RX_on_q3()
     )
 end
 
-function decompose_toffoli()
+function toffoli()
 
     function toffoli_circuit()
         # [(depth, gate)]
@@ -56,7 +56,7 @@ function decompose_toffoli()
     )
 end
 
-function decompose_toffoli_using_kronecker()
+function toffoli_using_kronecker()
 
     println(">>>>> Toffoli gate using Kronecker <<<<<")
  
@@ -73,7 +73,7 @@ function decompose_toffoli_using_kronecker()
     )
 end
 
-function decompose_toffoli_with_controlled_gates()
+function toffoli_with_controlled_gates()
 
     # Reference: https://doi.org/10.1109/TCAD.2005.858352
     println(">>>>> Toffoli gate with controlled gates <<<<<")
@@ -89,7 +89,7 @@ function decompose_toffoli_with_controlled_gates()
     )
 end
 
-function decompose_CNot_1_3()
+function CNot_1_3()
 
     return Dict{String, Any}(
 
@@ -102,7 +102,7 @@ function decompose_CNot_1_3()
     )
 end
 
-function decompose_FredkinGate()
+function FredkinGate()
 
     return Dict{String, Any}(
     "num_qubits" => 3,
@@ -115,7 +115,7 @@ function decompose_FredkinGate()
     )
 end
 
-function decompose_FredkinGate_using_SSwap()
+function FredkinGate_using_SSwap()
 
     return Dict{String, Any}(
     "num_qubits" => 3,
@@ -127,7 +127,7 @@ function decompose_FredkinGate_using_SSwap()
     )
 end
 
-function decompose_toffoli_left()
+function toffoli_left()
 
     # Reference: https://arxiv.org/pdf/0803.2316.pdf 
 
@@ -151,7 +151,7 @@ function decompose_toffoli_left()
     )
 end
 
-function decompose_toffoli_right()
+function toffoli_right()
 
     # Reference: https://arxiv.org/pdf/0803.2316.pdf 
     
@@ -177,7 +177,7 @@ function decompose_toffoli_right()
     )
 end
 
-function decompose_miller()
+function miller()
     # Reference: https://doi.org/10.1109/TCAD.2005.858352
 
     function target_gate()
@@ -201,7 +201,7 @@ function decompose_miller()
     )
 end
 
-function decompose_relative_toffoli()
+function relative_toffoli()
     #Reference: https://arxiv.org/pdf/1508.03273.pdf
 
     return Dict{String, Any}(
@@ -216,7 +216,7 @@ function decompose_relative_toffoli()
     )
 end
 
-function decompose_margolus()
+function margolus()
     #Reference: https://arxiv.org/pdf/quant-ph/0312225.pdf
 
     return Dict{String, Any}(
@@ -233,7 +233,7 @@ function decompose_margolus()
         )
 end
 
-function decompose_CiSwap()
+function CiSwap()
     #Reference: https://doi.org/10.1103/PhysRevResearch.2.033097
 
     return Dict{String, Any}(
diff --git a/examples/4qubit_gates.jl b/examples/4qubit_gates.jl
index 4f2bb09..2814a2c 100644
--- a/examples/4qubit_gates.jl
+++ b/examples/4qubit_gates.jl
@@ -1,4 +1,4 @@
-function decompose_CNot_41()
+function CNot_41()
     # Reference: https://doi.org/10.1109/DSD.2018.00005
 
     return Dict{String, Any}(
@@ -12,7 +12,7 @@ function decompose_CNot_41()
     )
 end
 
-function decompose_quantum_fulladder()
+function quantum_fulladder()
     #Reference-1: https://doi.org/10.1109/DATE.2005.249
     #Reference-2: https://doi.org/10.1109/TCAD.2005.858352 
     
@@ -42,7 +42,7 @@ function decompose_quantum_fulladder()
     )
 end
 
-function decompose_double_toffoli()
+function double_toffoli()
 
     # Reference: https://doi.org/10.1109/TCAD.2005.858352 
     
@@ -69,7 +69,7 @@ function decompose_double_toffoli()
     )
 end
 
-function decompose_double_peres()
+function double_peres()
 
     # Reference: https://doi.org/10.1109/TCAD.2005.858352 
     
diff --git a/examples/5qubit_gates.jl b/examples/5qubit_gates.jl
index 59c7220..a852bbe 100644
--- a/examples/5qubit_gates.jl
+++ b/examples/5qubit_gates.jl
@@ -1,4 +1,4 @@
-function decompose_RX_on_5qubits()
+function RX_on_5qubits()
 
     println(">>>>> RX Gate on fourth qubit using U3Gate <<<<<")
  
diff --git a/examples/Project.toml b/examples/Project.toml
index 3e411d7..5157379 100644
--- a/examples/Project.toml
+++ b/examples/Project.toml
@@ -1,7 +1,7 @@
 [deps]
 CPLEX = "a076750e-1247-5638-91d2-ce28b192dca0"
-Cbc = "9961bab8-2fa3-5c5a-9d89-47fab24efd76"
 Gurobi = "2e9cd046-0924-5485-92f1-d5272153d98b"
+HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 QuantumCircuitOpt = "88128e30-b60a-4e54-ab02-1050a5f92a36"
 Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
diff --git a/examples/README.md b/examples/README.md
new file mode 100644
index 0000000..897a99a
--- /dev/null
+++ b/examples/README.md
@@ -0,0 +1,2 @@
+# Examples for QuantumCircuitOpt
+`run_examples.jl` file contains the main driver file which can be executed to decompose a target gate into sequence of elementary gates of your choice. Multiple such gate implementations can be found inside files such as `2qubit_gates.jl`, `3qubit_gates.jl`, etc.
\ No newline at end of file
diff --git a/examples/decompose_all_gates.jl b/examples/decompose_all_gates.jl
new file mode 100644
index 0000000..b028582
--- /dev/null
+++ b/examples/decompose_all_gates.jl
@@ -0,0 +1,41 @@
+decompose_gates = [ 
+        #-- 2-qubit gates --#
+        "hadamard",
+        "controlled_Z",
+        "controlled_V",
+        "controlled_H",
+        "controlled_H_with_R",
+        "magic", 
+        "magic_using_CNOT_1_2",
+        "magic_using_SHCnot",
+        "S_gate", 
+        "revcnot", 
+        "revcnot_with_U", 
+        "swap",
+        "W_gate",
+        "W_using_HCnot",
+        "GroverDiffusion_using_Clifford",
+        "GroverDiffusion_using_U3",
+        "iSwap",
+        "qft2_using_HT",
+        "RX_on_q3",
+        "X_using_GR",
+        "CNOT_using_GR",
+        #-- 3-qubit gates --#
+        "toffoli_with_controlled_gates",
+        "CNot_1_3",
+        "FredkinGate",
+        "toffoli_left",
+        "toffoli_right",
+        "miller",
+        "relative_toffoli",
+        "margolus",
+        "CiSwap", # up to 1000s
+        #-- 4-qubit gates --#
+        "CNot_41",
+        "double_peres",
+        "quantum_fulladder",
+        "double_toffoli",
+        #-- 5-qubit gates --#
+        "RX_on_5qubits"
+        ]
\ No newline at end of file
diff --git a/examples/gates_sc21.jl b/examples/gates_sc21.jl
index 770ce38..f402512 100644
--- a/examples/gates_sc21.jl
+++ b/examples/gates_sc21.jl
@@ -7,7 +7,7 @@ Version of QCOpt: v0.3.0
 Last updated: Oct 12, 2021
 =#
 
-function decompose_controlled_Z()
+function controlled_Z()
 
     println(">>>>> Controlled-Z Gate <<<<<")
 
@@ -24,7 +24,7 @@ function decompose_controlled_Z()
     "U3_λ_discretization" => -π:π/2:π,)    
 end
 
-function decompose_controlled_V()
+function controlled_V()
 
     println(">>>>> Controlled-V Gate <<<<<")
 
@@ -37,7 +37,7 @@ function decompose_controlled_V()
     "decomposition_type" => "exact_optimal")
 end
 
-function decompose_controlled_H()
+function controlled_H()
 
     println(">>>>> Controlled-H Gate <<<<<")
 
@@ -55,7 +55,7 @@ function decompose_controlled_H()
     )
 end
 
-function decompose_magic_using_U3_CNot_2_1()
+function magic_using_U3_CNot_2_1()
     
     println(">>>>> Magic basis using U3 and CNot_2_1 <<<<<")
 
@@ -73,7 +73,7 @@ function decompose_magic_using_U3_CNot_2_1()
     )     
 end
 
-function decompose_iSwapGate()
+function iSwapGate()
 
     println(">>>>> iSwap Gate <<<<<")
 
@@ -89,7 +89,7 @@ function decompose_iSwapGate()
     
 end
 
-function decompose_GroverDiffusion_using_Clifford()
+function GroverDiffusion_using_Clifford()
 
     println(">>>>> Grover's Diffusion Operator <<<<<")
 
@@ -103,7 +103,7 @@ function decompose_GroverDiffusion_using_Clifford()
 end
 
 
-function decompose_GroverDiffusion_using_U3()
+function GroverDiffusion_using_U3()
 
     println(">>>>> Grover's Diffusion Operator using U3 gate <<<<<")
 
@@ -120,7 +120,7 @@ function decompose_GroverDiffusion_using_U3()
         "U3_λ_discretization" => [0])
 end
 
-function decompose_magic_using_SH_CNot_1_2()
+function magic_using_SH_CNot_1_2()
     
     println(">>>>> M gate using S, H and CNOT_1_2 Gate <<<<<")
 
@@ -134,7 +134,7 @@ function decompose_magic_using_SH_CNot_1_2()
 end
 
 
-function decompose_magic_using_SH_CNot_2_1()
+function magic_using_SH_CNot_2_1()
     
     println(">>>>> M gate using S, H and CNOT_2_1 Gate <<<<<")
 
@@ -147,7 +147,7 @@ function decompose_magic_using_SH_CNot_2_1()
         "decomposition_type" => "exact_optimal")
 end
 
-function decompose_magic_using_U3_CNot_1_2()
+function magic_using_U3_CNot_1_2()
     
     println(">>>>> Magic basis using U3 and CNot_1_2 <<<<<")
 
@@ -164,7 +164,7 @@ function decompose_magic_using_U3_CNot_1_2()
         "U3_λ_discretization" => -π:π/2:π,)
 end
 
-function decompose_toffoli_with_controlled_gates()
+function toffoli_with_controlled_gates()
 
     # Reference: https://doi.org/10.1109/TCAD.2005.858352
     println(">>>>> Toffoli gate using controlled gates <<<<<")
@@ -178,7 +178,7 @@ function decompose_toffoli_with_controlled_gates()
     "decomposition_type" => "exact_optimal")
 end
 
-function decompose_Fredkin()
+function Fredkin()
 
     println(">>>>> Fredkin gate using controlled gates <<<<<")
 
@@ -192,7 +192,7 @@ function decompose_Fredkin()
     "decomposition_type" => "exact_optimal")
 end
 
-function decompose_double_toffoli()
+function double_toffoli()
 
     println(">>>>> Double Toffoli using controlled gates <<<<<")
 
@@ -219,7 +219,7 @@ function decompose_double_toffoli()
     "decomposition_type" => "exact_optimal")
 end
 
-function decompose_quantum_fulladder()
+function quantum_fulladder()
     #Reference-1: https://doi.org/10.1109/DATE.2005.249
     #Reference-2: https://doi.org/10.1109/TCAD.2005.858352 
     
diff --git a/examples/optimizers.jl b/examples/optimizers.jl
index bbf1c20..c9b269a 100644
--- a/examples/optimizers.jl
+++ b/examples/optimizers.jl
@@ -1,7 +1,8 @@
-#=====================================#
-# MIP solvers (commercial, but fast)  #
-#=====================================#
+#===================================#
+# MIP solvers (commercial, faster)  #
+#===================================#
 
+# https://github.com/jump-dev/Gurobi.jl
 function get_gurobi()
     return JuMP.optimizer_with_attributes(Gurobi.Optimizer, 
                                           MOI.Silent() => false, 
@@ -9,6 +10,7 @@ function get_gurobi()
                                           "Presolve" => 1) 
 end
 
+# https://github.com/jump-dev/CPLEX.jl
 function get_cplex()
      return JuMP.optimizer_with_attributes(CPLEX.Optimizer, 
                                            MOI.Silent() => false, 
@@ -17,15 +19,25 @@ function get_cplex()
                                            "CPX_PARAM_PREIND" => 1)
 end
 
-#======================================#
-# MIP solvers (open-source, but slow)  #
-#======================================#
+#====================================#
+# MIP solvers (open-source, slower)  #
+#====================================#
+# https://github.com/jump-dev/HiGHS.jl
+function get_highs()
+    return JuMP.optimizer_with_attributes(
+        HiGHS.Optimizer,
+        "presolve" => "on",
+        "log_to_console" => true,
+    )
+end
 
+# https://github.com/jump-dev/Cbc.jl
 function get_cbc()
     return JuMP.optimizer_with_attributes(Cbc.Optimizer, 
                                           MOI.Silent() => false) 
 end
 
+# https://github.com/jump-dev/Glpk.jl
 function get_glpk()
     return JuMP.optimizer_with_attributes(GLPK.Optimizer, 
                                           MOI.Silent() => false)
@@ -34,7 +46,7 @@ end
 #========================================================#
 # Continuous nonlinear programming solver (open-source)  #
 #========================================================#
-
+# https://github.com/jump-dev/Ipopt.jl
 function get_ipopt()
      return JuMP.optimizer_with_attributes(Ipopt.Optimizer, 
                                        MOI.Silent() => true, 
@@ -46,6 +58,7 @@ end
 #=================================================================#
 # Local mixed-integer nonlinear programming solver (open-source)  #
 #=================================================================#
+# https://github.com/lanl-ansi/Juniper.jl
 function get_juniper()
      return JuMP.optimizer_with_attributes(Juniper.Optimizer, 
                                          MOI.Silent() => false, 
@@ -56,13 +69,11 @@ end
 #=================================================================#
 # Global mixed-integer nonlinear programming solver (open-source) #
 #=================================================================#
+# https://github.com/lanl-ansi/Alpine.jl
 function get_alpine()
         return JuMP.optimizer_with_attributes(Alpine.Optimizer, 
                                             "nlp_solver" => get_ipopt(),
                                             "minlp_solver" => get_juniper(),  
-                                            "mip_solver" => get_gurobi(),
-                                            "optimizer_presolve_bt" => false,
-                                            "optimizer_presolve_max_iter" => 10,
-                                            "optimizer_presolve_bp" => false,
-                                            "disc_ratio" => 10)
+                                            "mip_solver" => get_gurobi()
+        )
 end
\ No newline at end of file
diff --git a/examples/results_sc21.jl b/examples/results_sc21.jl
index 85e304b..773452e 100644
--- a/examples/results_sc21.jl
+++ b/examples/results_sc21.jl
@@ -19,12 +19,12 @@ qcm_optimizer = JuMP.optimizer_with_attributes(Gurobi.Optimizer, "Presolve" => 1
 #-----------------------------#
 #      Results of Table-I     #
 #-----------------------------#
-table_I_gates = ["decompose_controlled_Z",
-                 "decompose_controlled_V",
-                 "decompose_controlled_H",
-                 "decompose_magic_using_U3_CNot_2_1",
-                 "decompose_iSwapGate",
-                 "decompose_GroverDiffusion_using_U3"]
+table_I_gates = ["controlled_Z",
+                 "controlled_V",
+                 "controlled_H",
+                 "magic_using_U3_CNot_2_1",
+                 "iSwapGate",
+                 "GroverDiffusion_using_U3"]
 
 result_with_VC = Dict{String,Any}()
 result_without_VC = Dict{String,Any}()
@@ -54,10 +54,10 @@ end
 #      Results of Figure 3 (Magic basis)     #
 #--------------------------------------------#
 result = Dict{String,Any}()
-figure_3_gates = ["decompose_magic_using_SH_CNot_2_1",
-                  "decompose_magic_using_U3_CNot_2_1",
-                  "decompose_magic_using_SH_CNot_1_2",
-                  "decompose_magic_using_U3_CNot_1_2",
+figure_3_gates = ["magic_using_SH_CNot_2_1",
+                  "magic_using_U3_CNot_2_1",
+                  "magic_using_SH_CNot_1_2",
+                  "magic_using_U3_CNot_1_2",
                   ]
 
 run_times_fig3 = zeros(length(figure_3_gates),1)
@@ -76,8 +76,8 @@ end
 #      Results of Figure 4 (Grover operator)     #
 #------------------------------------------------#
 result = Dict{String,Any}()
-figure_4_gates = ["decompose_GroverDiffusion_using_Clifford",
-                  "decompose_GroverDiffusion_using_U3"]
+figure_4_gates = ["GroverDiffusion_using_Clifford",
+                  "GroverDiffusion_using_U3"]
 
 run_times_fig4 = zeros(length(figure_4_gates),1)
 
@@ -95,10 +95,10 @@ end
 #      Results of Figure 5 (larger circuits)     #
 #------------------------------------------------#
 result = Dict{String,Any}()
-figure_5_gates = ["decompose_toffoli_with_controlled_gates",
-                  "decompose_Fredkin",
-                  "decompose_double_toffoli",
-                  "decompose_quantum_fulladder",
+figure_5_gates = ["toffoli_with_controlled_gates",
+                  "Fredkin",
+                  "double_toffoli",
+                  "quantum_fulladder",
                   ]
 
 run_times_fig5 = zeros(length(figure_5_gates),1)
diff --git a/examples/run_examples.jl b/examples/run_examples.jl
index 5e68861..fecab93 100644
--- a/examples/run_examples.jl
+++ b/examples/run_examples.jl
@@ -1,50 +1,35 @@
 import QuantumCircuitOpt as QCOpt
 using JuMP
 using Gurobi
-# using Cbc
-# using CPLEX
+# using HiGHS
 
 include("optimizers.jl")
 include("2qubit_gates.jl")
 include("3qubit_gates.jl")
 include("4qubit_gates.jl")
 include("5qubit_gates.jl")
+include("decompose_all_gates.jl")
 
-decompose_gates = ["decompose_hadamard",
-                   "decompose_controlled_Z",
-                   "decompose_controlled_V",
-                   "decompose_controlled_H",
-                   "decompose_controlled_H_with_R",
-                   "decompose_magic", 
-                   "decompose_magic_using_SHCnot",
-                   "decompose_S", 
-                   "decompose_revcnot", 
-                   "decompose_revcnot_with_U", 
-                   "decompose_swap",
-                   "decompose_W",
-                   "decompose_W_using_HCnot",
-                   "decompose_GroverDiffusion_using_Clifford",
-                   "decompose_GroverDiffusion_using_U3",
-                   "decompose_iSwap",
-                   "decompose_qft2_using_HT",
-                   "decompose_RX_on_q3"]
-
-# decompose_gates = ["decompose_controlled_V"]
+# decompose_gates = ["iSwap"]
 
 #----------------------------------------------#
 #      Quantum Circuit Optimization model      #
 #----------------------------------------------#
+qcopt_optimizer = get_gurobi()
+
 result = Dict{String,Any}()
+times = zeros(length(decompose_gates), 1)
+
+for gates = 1:length(decompose_gates)
+    params = getfield(Main, Symbol(decompose_gates[gates]))()
 
-for gates in decompose_gates 
-    
-    params = getfield(Main, Symbol(gates))()
+    model_options = Dict{Symbol, Any}(
+        :convex_hull_gate_constraints => false,
+        :idempotent_gate_constraints  => false,
+        :unitary_constraints          => false,
+        :fix_unitary_variables        => true,
+    )
 
-    model_options = Dict{Symbol, Any}(:convex_hull_gate_constraints => false,
-                                      :unit_magnitude_constraints   => false,
-                                      :idempotent_gate_constraints  => false,
-                                      :unitary_constraints          => false)
-    
-    qcopt_optimizer = get_gurobi()
     global result = QCOpt.run_QCModel(params, qcopt_optimizer; options = model_options)
+    times[gates] = result["solve_time"]
 end
diff --git a/src/QuantumCircuitOpt.jl b/src/QuantumCircuitOpt.jl
index 02df84d..7df2385 100644
--- a/src/QuantumCircuitOpt.jl
+++ b/src/QuantumCircuitOpt.jl
@@ -8,6 +8,7 @@ const LA = LinearAlgebra
 const QCO = QuantumCircuitOpt
 const kron_symbol = 'x'
 const qubit_separator = '_'
+const _QCOPT_VERSION = "v0.5.0"
 
 # Create our module level logger (this will get precompiled)
 const _LOGGER = Memento.getlogger(@__MODULE__)
diff --git a/src/constraints.jl b/src/constraints.jl
index caa3c80..698ee01 100644
--- a/src/constraints.jl
+++ b/src/constraints.jl
@@ -5,9 +5,9 @@
 function constraint_single_gate_per_depth(qcm::QuantumCircuitModel)
 
     num_gates = size(qcm.data["gates_real"])[3]
-    depth     = qcm.data["maximum_depth"]
+    max_depth = qcm.data["maximum_depth"]
     
-    JuMP.@constraint(qcm.model, [d=1:depth], sum(qcm.variables[:z_bin_var][n,d] for n=1:num_gates) == 1)
+    JuMP.@constraint(qcm.model, [d=1:max_depth], sum(qcm.variables[:z_bin_var][n,d] for n=1:num_gates) == 1)
     
     return
 end
@@ -15,9 +15,9 @@ end
 function constraint_gates_onoff_per_depth(qcm::QuantumCircuitModel)
 
     num_gates = size(qcm.data["gates_real"])[3]
-    depth     = qcm.data["maximum_depth"]
+    max_depth = qcm.data["maximum_depth"]
 
-    JuMP.@constraint(qcm.model, [d=1:depth], qcm.variables[:G_var][:,:,d] .== 
+    JuMP.@constraint(qcm.model, [d=1:max_depth], qcm.variables[:G_var][:,:,d] .== 
                                     sum(qcm.variables[:z_bin_var][n,d] * qcm.data["gates_real"][:,:,n] for n=1:num_gates))
     
     return
@@ -36,15 +36,15 @@ end
 function constraint_gate_intermediate_products(qcm::QuantumCircuitModel)
 
     num_gates = size(qcm.data["gates_real"])[3]
-    depth     = qcm.data["maximum_depth"]
+    max_depth = qcm.data["maximum_depth"]
 
-    JuMP.@constraint(qcm.model, [d=1:(depth-1)], qcm.variables[:U_var][:,:,d] .== 
+    JuMP.@constraint(qcm.model, [d=1:(max_depth-1)], qcm.variables[:U_var][:,:,d] .== 
                                 sum(qcm.variables[:V_var][:,:,n,d] * qcm.data["gates_real"][:,:,n] for n=1:num_gates))
 
-    JuMP.@constraint(qcm.model, [d=2:depth], sum(qcm.variables[:V_var][:,:,n,d] for n=1:num_gates) .== 
+    JuMP.@constraint(qcm.model, [d=2:max_depth], sum(qcm.variables[:V_var][:,:,n,d] for n=1:num_gates) .== 
                                                 qcm.variables[:U_var][:,:,(d-1)])
 
-    JuMP.@constraint(qcm.model, [d=1:(depth-1)], qcm.variables[:V_var][:,:,:,(d+1)] .== 
+    JuMP.@constraint(qcm.model, [d=1:(max_depth-1)], qcm.variables[:V_var][:,:,:,(d+1)] .== 
                                          qcm.variables[:zU_var][:,:,:,d])
     
     return
@@ -52,16 +52,15 @@ end
 
 function constraint_gate_target_condition(qcm::QuantumCircuitModel)
 
-    depth              = qcm.data["maximum_depth"]
+    max_depth          = qcm.data["maximum_depth"]
     num_gates          = size(qcm.data["gates_real"])[3]
     decomposition_type = qcm.data["decomposition_type"]
     
-    # For correct implementation of this, use MutableArithmetics.jl >= v0.2.11
     if decomposition_type in ["exact_optimal", "exact_feasible"]
-        JuMP.@constraint(qcm.model, sum(qcm.variables[:V_var][:,:,n,depth] * qcm.data["gates_real"][:,:,n] for n=1:num_gates) .== qcm.data["target_gate"])
+        JuMP.@constraint(qcm.model, sum(qcm.variables[:V_var][:,:,n,max_depth] * qcm.data["gates_real"][:,:,n] for n=1:num_gates) .== qcm.data["target_gate"])
     
     elseif decomposition_type == "approximate"
-        JuMP.@constraint(qcm.model, sum(qcm.variables[:V_var][:,:,n,depth] * qcm.data["gates_real"][:,:,n] for n=1:num_gates) .== qcm.data["target_gate"][:,:] + qcm.variables[:slack_var][:,:])  
+        JuMP.@constraint(qcm.model, sum(qcm.variables[:V_var][:,:,n,max_depth] * qcm.data["gates_real"][:,:,n] for n=1:num_gates) .== qcm.data["target_gate"][:,:] + qcm.variables[:slack_var][:,:])  
 
     end
 
@@ -70,41 +69,47 @@ end
 
 function constraint_complex_to_real_symmetry(qcm::QuantumCircuitModel)
 
-    depth  = qcm.data["maximum_depth"]
+    max_depth  = qcm.data["maximum_depth"]
     n_r    = size(qcm.data["gates_real"])[1]
     n_c    = size(qcm.data["gates_real"])[2]
-
-    JuMP.@constraint(qcm.model, [i=1:2:n_r, j=1:2:n_c, d=1:(depth-1)], qcm.variables[:U_var][i,j,d]   ==  qcm.variables[:U_var][i+1,j+1,d])
-    JuMP.@constraint(qcm.model, [i=1:2:n_r, j=1:2:n_c, d=1:(depth-1)], qcm.variables[:U_var][i,j+1,d] == -qcm.variables[:U_var][i+1,j,d])
     
+    for i=1:2:n_r, j=1:2:n_c, d=1:(max_depth-1)
+        JuMP.@constraint(qcm.model, qcm.variables[:U_var][i,j,d] == qcm.variables[:U_var][i+1,j+1,d])
+      JuMP.@constraint(qcm.model, qcm.variables[:U_var][i,j+1,d] == -qcm.variables[:U_var][i+1,j,d]) # This seems to slow down the solution search
+    end
+    
+
     return
 end
 
 function constraint_gate_product_linearization(qcm::QuantumCircuitModel)
 
-    depth          = qcm.data["maximum_depth"]
+    max_depth      = qcm.data["maximum_depth"]
     n_r            = size(qcm.data["gates_real"])[1]
     n_c            = size(qcm.data["gates_real"])[2]
     num_gates      = size(qcm.data["gates_real"])[3]
     are_gates_real = qcm.data["are_gates_real"]
 
+    U_var  = qcm.variables[:U_var]
+    zU_var = qcm.variables[:zU_var]
+    z_bin_var = qcm.variables[:z_bin_var]
+
     i_val = 2 - are_gates_real
+    U_var_bound_tol = 1E-8
 
-    for i=1:i_val:n_r
-        for j=1:n_c
-            for n=1:num_gates
-                for d=1:(depth-1)
-                    
-                    QCO.relaxation_bilinear(qcm.model, qcm.variables[:zU_var][i,j,n,d], qcm.variables[:U_var][i,j,d], qcm.variables[:z_bin_var][n,(d+1)])
-                    if !are_gates_real
-                        if isodd(j)
-                            JuMP.@constraint(qcm.model, qcm.variables[:zU_var][i,j,n,d]   ==  qcm.variables[:zU_var][i+1,j+1,n,d])
-                            JuMP.@constraint(qcm.model, qcm.variables[:zU_var][i,j+1,n,d] == -qcm.variables[:zU_var][i+1,j,n,d])
-                        end
-                    end
+    for i=1:i_val:n_r, j=1:n_c, n=1:num_gates, d=1:(max_depth-1)
+        U_var_l = JuMP.lower_bound(U_var[i,j,d])
+        U_var_u = JuMP.upper_bound(U_var[i,j,d])
+        
+        if !(isapprox(abs(U_var_u - U_var_l), 0, atol = 2*U_var_bound_tol))
+            QCO.relaxation_bilinear(qcm.model, zU_var[i,j,n,d], U_var[i,j,d], z_bin_var[n,(d+1)])
+        else
+            JuMP.@constraint(qcm.model, zU_var[i,j,n,d] == (U_var_l + U_var_u)/2 * z_bin_var[n,(d+1)])
+        end
 
-                end
-            end
+        if !are_gates_real && isodd(j)
+            JuMP.@constraint(qcm.model, zU_var[i,j,n,d]   ==  zU_var[i+1,j+1,n,d])
+            JuMP.@constraint(qcm.model, zU_var[i,j+1,n,d] == -zU_var[i+1,j,n,d])
         end
     end
     
@@ -124,9 +129,9 @@ end
 function constraint_gate_intermediate_products_compact(qcm::QuantumCircuitModel)
 
     num_gates = size(qcm.data["gates_real"])[3]
-    depth     = qcm.data["maximum_depth"]
-    
-    JuMP.@constraint(qcm.model, [d=2:(depth-1)], qcm.variables[:U_var][:,:,d] .== 
+    max_depth = qcm.data["maximum_depth"]
+
+    JuMP.@constraint(qcm.model, [d=2:max_depth], qcm.variables[:U_var][:,:,d] .== 
                                         sum((qcm.variables[:zU_var][:,:,n,(d-1)]) * qcm.data["gates_real"][:,:,n] for n=1:num_gates))
 
     return
@@ -134,27 +139,52 @@ end
 
 function constraint_gate_target_condition_compact(qcm::QuantumCircuitModel)
 
-    depth   = qcm.data["maximum_depth"]
-    num_gates = size(qcm.data["gates_real"])[3]
+    max_depth = qcm.data["maximum_depth"]
     decomposition_type = qcm.data["decomposition_type"]
 
-    zU_var = qcm.variables[:zU_var]
+    U_var = qcm.variables[:U_var]
     
     # For correct implementation of this, use MutableArithmetics.jl >= v0.2.11
     if decomposition_type in ["exact_optimal", "exact_feasible"]
-        JuMP.@constraint(qcm.model, sum(zU_var[:,:,n,(depth-1)] * qcm.data["gates_real"][:,:,n] for n=1:num_gates) .== qcm.data["target_gate"][:,:])  
+        JuMP.@constraint(qcm.model, U_var[:,:,max_depth] .== qcm.data["target_gate"][:,:])
     
     elseif decomposition_type == "approximate"
-        JuMP.@constraint(qcm.model, sum(zU_var[:,:,n,(depth-1)] * qcm.data["gates_real"][:,:,n] for n=1:num_gates) .== qcm.data["target_gate"][:,:] + qcm.variables[:slack_var][:,:])    
+        JuMP.@constraint(qcm.model, U_var[:,:,max_depth] .== qcm.data["target_gate"][:,:] + qcm.variables[:slack_var][:,:])
     
     end
     
     return
 end
 
+function constraint_slack_var_outer_approximation(qcm::QuantumCircuitModel)
+    slack_var    = qcm.variables[:slack_var]
+    slack_var_oa = qcm.variables[:slack_var_oa]
+    
+    # Number of under-estimators for the quadratic function
+    num_points = 9 
+
+    for i=1:size(slack_var)[1], j=1:size(slack_var)[2]
+        lb = JuMP.lower_bound(slack_var[i,j])
+        ub = JuMP.upper_bound(slack_var[i,j])
+        if !(isapprox(lb, ub, atol = 1E-6))
+            oa_points = range(lb, ub, num_points)
+            JuMP.@constraint(qcm.model, [k=1:length(oa_points)], 
+                            slack_var_oa[i,j] >= 2*slack_var[i,j]*oa_points[k] - (oa_points[k])^2)
+            if !(0 in oa_points)
+                JuMP.@constraint(qcm.model, slack_var_oa[i,j] >= 0)
+            end
+        else
+            mid_point = (lb+ub)/2
+            JuMP.@constraint(qcm.model, slack_var_oa[i,j] >= 2*slack_var[i,j]*mid_point - (mid_point)^2)
+        end
+    end
+
+    return
+end
+
 function constraint_commutative_gate_pairs(qcm::QuantumCircuitModel)
     
-    depth     = qcm.data["maximum_depth"]
+    max_depth = qcm.data["maximum_depth"]
     z_bin_var = qcm.variables[:z_bin_var]
 
     commute_pairs, commute_pairs_prodIdentity = QCO.get_commutative_gate_pairs(qcm.data["gates_dict"])
@@ -164,7 +194,7 @@ function constraint_commutative_gate_pairs(qcm::QuantumCircuitModel)
         (length(commute_pairs) > 1)  && (Memento.info(_LOGGER, "Detected $(length(commute_pairs)) input elementary gate pairs which commute"))
 
         for i = 1:length(commute_pairs)
-            JuMP.@constraint(qcm.model, [d=1:(depth-1)], z_bin_var[commute_pairs[i][2], d] + z_bin_var[commute_pairs[i][1], d+1] <= 1)
+            JuMP.@constraint(qcm.model, [d=1:(max_depth-1)], z_bin_var[commute_pairs[i][2], d] + z_bin_var[commute_pairs[i][1], d+1] <= 1)
         end
     end
 
@@ -173,8 +203,8 @@ function constraint_commutative_gate_pairs(qcm::QuantumCircuitModel)
         (length(commute_pairs_prodIdentity) > 1)  && (Memento.info(_LOGGER, "Detected $(length(commute_pairs_prodIdentity)) input elementary gate pairs whose product is Identity"))
 
         for i = 1:length(commute_pairs_prodIdentity)
-            JuMP.@constraint(qcm.model, [d=1:(depth-1)], z_bin_var[commute_pairs_prodIdentity[i][2], d] + z_bin_var[commute_pairs_prodIdentity[i][1], d+1] <= 1)
-            JuMP.@constraint(qcm.model, [d=1:(depth-1)], z_bin_var[commute_pairs_prodIdentity[i][1], d] + z_bin_var[commute_pairs_prodIdentity[i][2], d+1] <= 1)
+            JuMP.@constraint(qcm.model, [d=1:(max_depth-1)], z_bin_var[commute_pairs_prodIdentity[i][2], d] + z_bin_var[commute_pairs_prodIdentity[i][1], d+1] <= 1)
+            JuMP.@constraint(qcm.model, [d=1:(max_depth-1)], z_bin_var[commute_pairs_prodIdentity[i][1], d] + z_bin_var[commute_pairs_prodIdentity[i][2], d+1] <= 1)
         end
     end
 
@@ -184,7 +214,7 @@ end
 function constraint_involutory_gates(qcm::QuantumCircuitModel)
 
     gates_dict = qcm.data["gates_dict"]
-    depth      = qcm.data["maximum_depth"]
+    max_depth  = qcm.data["maximum_depth"]
     z_bin_var  = qcm.variables[:z_bin_var]
 
     involutory_gates = QCO.get_involutory_gates(gates_dict)
@@ -194,7 +224,7 @@ function constraint_involutory_gates(qcm::QuantumCircuitModel)
         (length(involutory_gates) > 1)  && (Memento.info(_LOGGER, "Detected $(length(involutory_gates)) involutory elementary gates"))
 
         for i = 1:length(involutory_gates)
-            JuMP.@constraint(qcm.model, [d=1:(depth-1)], z_bin_var[involutory_gates[i], d] + z_bin_var[involutory_gates[i], d+1] <= 1)
+            JuMP.@constraint(qcm.model, [d=1:(max_depth-1)], z_bin_var[involutory_gates[i], d] + z_bin_var[involutory_gates[i], d+1] <= 1)
         end
     end
 
@@ -204,7 +234,7 @@ end
 function constraint_redundant_gate_product_pairs(qcm::QuantumCircuitModel)
 
     gates_dict = qcm.data["gates_dict"]
-    depth      = qcm.data["maximum_depth"]
+    max_depth  = qcm.data["maximum_depth"]
     z_bin_var  = qcm.variables[:z_bin_var]
 
     redundant_pairs = QCO.get_redundant_gate_product_pairs(gates_dict)
@@ -214,7 +244,7 @@ function constraint_redundant_gate_product_pairs(qcm::QuantumCircuitModel)
         (length(redundant_pairs) > 1)  && (Memento.info(_LOGGER, "Detected $(length(redundant_pairs)) redundant input elementary gate pairs"))
 
         for i = 1:length(redundant_pairs)
-            JuMP.@constraint(qcm.model, [d=1:(depth-1)], z_bin_var[redundant_pairs[i][2], d] + z_bin_var[redundant_pairs[i][1], d+1] <= 1)
+            JuMP.@constraint(qcm.model, [d=1:(max_depth-1)], z_bin_var[redundant_pairs[i][2], d] + z_bin_var[redundant_pairs[i][1], d+1] <= 1)
         end
     end
 
@@ -224,7 +254,7 @@ end
 function constraint_idempotent_gates(qcm::QuantumCircuitModel)
 
     gates_dict = qcm.data["gates_dict"]
-    depth      = qcm.data["maximum_depth"]
+    max_depth  = qcm.data["maximum_depth"]
     z_bin_var  = qcm.variables[:z_bin_var]
 
     idempotent_gates = QCO.get_idempotent_gates(gates_dict)
@@ -234,7 +264,7 @@ function constraint_idempotent_gates(qcm::QuantumCircuitModel)
         (length(idempotent_gates) > 1)  && (Memento.info(_LOGGER, "Detected $(length(idempotent_gates)) idempotent elementary gates"))
 
         for i = 1:length(idempotent_gates)
-            JuMP.@constraint(qcm.model, [d=1:(depth-1)], z_bin_var[idempotent_gates[i], d] + z_bin_var[idempotent_gates[i], d+1] <= 1)
+            JuMP.@constraint(qcm.model, [d=1:(max_depth-1)], z_bin_var[idempotent_gates[i], d] + z_bin_var[idempotent_gates[i], d+1] <= 1)
         end
     end
 
@@ -244,7 +274,7 @@ end
 function constraint_identity_gate_symmetry(qcm::QuantumCircuitModel)
 
     gates_dict = qcm.data["gates_dict"]
-    depth      = qcm.data["maximum_depth"]
+    max_depth  = qcm.data["maximum_depth"]
     z_bin_var  = qcm.variables[:z_bin_var]
 
     identity_idx = []
@@ -256,7 +286,7 @@ function constraint_identity_gate_symmetry(qcm::QuantumCircuitModel)
     
     if !isempty(identity_idx)
         for i = 1:length(identity_idx)
-            JuMP.@constraint(qcm.model, [d=1:(depth-1)], z_bin_var[identity_idx[i], d] <= z_bin_var[identity_idx[i], d+1])
+            JuMP.@constraint(qcm.model, [d=1:(max_depth-1)], z_bin_var[identity_idx[i], d] <= z_bin_var[identity_idx[i], d+1])
         end
     end
 
@@ -266,17 +296,17 @@ end
 function constraint_cnot_gate_bounds(qcm::QuantumCircuitModel)
 
     cnot_idx  = qcm.data["cnot_idx"]
-    depth     = qcm.data["maximum_depth"]
+    max_depth = qcm.data["maximum_depth"]
     z_bin_var = qcm.variables[:z_bin_var]
 
     if !isempty(cnot_idx)
         if "cnot_lower_bound" in keys(qcm.data)
-            JuMP.@constraint(qcm.model, sum(z_bin_var[n,d] for n in cnot_idx, d=1:depth) >= qcm.data["cnot_lower_bound"])
+            JuMP.@constraint(qcm.model, sum(z_bin_var[n,d] for n in cnot_idx, d=1:max_depth) >= qcm.data["cnot_lower_bound"])
             Memento.info(_LOGGER, "Applied CNot-gate lower bound constraint")
         end
         
         if "cnot_upper_bound" in keys(qcm.data)
-            JuMP.@constraint(qcm.model, sum(z_bin_var[n,d] for n in cnot_idx, d=1:depth) <= qcm.data["cnot_upper_bound"])
+            JuMP.@constraint(qcm.model, sum(z_bin_var[n,d] for n in cnot_idx, d=1:max_depth) <= qcm.data["cnot_upper_bound"])
             Memento.info(_LOGGER, "Applied CNot-gate upper bound constraint")
         end
     end
@@ -295,30 +325,28 @@ function constraint_convex_hull_complex_gates(qcm::QuantumCircuitModel)
         gates_real = qcm.data["gates_real"]
         gates_dict = qcm.data["gates_dict"]
         num_gates  = size(gates_real)[3]
-        depth      = qcm.data["maximum_depth"]
+        max_depth  = qcm.data["maximum_depth"]
         n_r        = size(gates_dict["1"]["matrix"])[1]
         n_c        = size(gates_dict["1"]["matrix"])[2]
 
         num_facets = 0
 
         vertices_dict = Dict{Set{}, Tuple{<:Number, <:Number}}()
-        for I=1:n_r
-            for J=1:n_c
-                vertices_coord_IJ = Set()
+        for I=1:n_r, J=1:n_c
+            vertices_coord_IJ = Set()
 
-                for K in keys(gates_dict)
-                    re = QCO.round_real_value(real(gates_dict[K]["matrix"][I,J]))
-                    im = QCO.round_real_value(imag(gates_dict[K]["matrix"][I,J]))
+            for K in keys(gates_dict)
+                re = QCO.round_real_value(real(gates_dict[K]["matrix"][I,J]))
+                im = QCO.round_real_value(imag(gates_dict[K]["matrix"][I,J]))
 
-                    push!(vertices_coord_IJ, (re, im))
-                end
-                
-                if (isapprox(minimum([x[1] for x in vertices_coord_IJ]), maximum([x[1] for x in vertices_coord_IJ]), atol = 1E-6)) || (isapprox(minimum([x[2] for x in vertices_coord_IJ]), maximum([x[2] for x in vertices_coord_IJ]), atol = 1E-6))
-                    continue
-                else 
-                    # Eliminates repeated set of extreme points
-                    vertices_dict[vertices_coord_IJ] = (I,J)
-                end
+                push!(vertices_coord_IJ, (re, im))
+            end
+            
+            if (isapprox(minimum([x[1] for x in vertices_coord_IJ]), maximum([x[1] for x in vertices_coord_IJ]), atol = 1E-6)) || (isapprox(minimum([x[2] for x in vertices_coord_IJ]), maximum([x[2] for x in vertices_coord_IJ]), atol = 1E-6))
+                continue
+            else 
+                # Eliminates repeated set of extreme points
+                vertices_dict[vertices_coord_IJ] = (I,J)
             end
         end
 
@@ -337,15 +365,15 @@ function constraint_convex_hull_complex_gates(qcm::QuantumCircuitModel)
 
                 if !isinf(slope)
                     if isapprox(abs(slope), 0, atol=1E-6)
-                        JuMP.@constraint(qcm.model, [d=1:depth], 
+                        JuMP.@constraint(qcm.model, [d=1:max_depth], 
                                         sum(gates_real[(2*I-1),(2*J), n_g] * z_bin_var[n_g,d] for n_g = 1:num_gates) - intercept == 0)
                     else
                         
-                        JuMP.@constraint(qcm.model, [d=1:depth], sum(gates_real[(2*I-1),(2*J), n_g] * z_bin_var[n_g,d] for n_g = 1:num_gates) 
+                        JuMP.@constraint(qcm.model, [d=1:max_depth], sum(gates_real[(2*I-1),(2*J), n_g] * z_bin_var[n_g,d] for n_g = 1:num_gates) 
                                                                 - slope*sum(gates_real[(2*I-1),(2*J-1), n_g] * z_bin_var[n_g,d] for n_g = 1:num_gates) - intercept == 0)
                     end
                 elseif isinf(slope)
-                    JuMP.@constraint(qcm.model, [d=1:depth], sum(gates_real[(2*I-1),(2*J-1), n_g] * z_bin_var[n_g,d] for n_g = 1:num_gates) == vertices[1][1])
+                    JuMP.@constraint(qcm.model, [d=1:max_depth], sum(gates_real[(2*I-1),(2*J-1), n_g] * z_bin_var[n_g,d] for n_g = 1:num_gates) == vertices[1][1])
                 end
                 
                 num_facets += 1
@@ -357,11 +385,6 @@ function constraint_convex_hull_complex_gates(qcm::QuantumCircuitModel)
                 end
                 
                 vertices_convex_hull = QCO.convex_hull(vertices)
-                
-                # if length(vertices_convex_hull) == 2
-                #     @show vertices_convex_hull
-                # end
-
                 num_ex_pt = size(vertices_convex_hull)[1]
                 
                 # Add a planar hull cut if num_ex_pt == 2
@@ -393,12 +416,10 @@ function constraint_convex_hull_complex_gates(qcm::QuantumCircuitModel)
 
                             if v3[2] - slope*v3[1] - intercept <= -1E-6
 
-                                if isapprox(abs(slope), 0, atol=1E-6)
-                                    
-                                    JuMP.@constraint(qcm.model, [d=1:depth], sum(gates_real[(2*I-1),(2*J), n_g] * z_bin_var[n_g,d] for n_g = 1:num_gates) - intercept <= 0)
-                                else
-                                    
-                                    JuMP.@constraint(qcm.model, [d=1:depth], sum(gates_real[(2*I-1),(2*J), n_g] * z_bin_var[n_g,d] for n_g = 1:num_gates) 
+                                if isapprox(abs(slope), 0, atol=1E-6)             
+                                    JuMP.@constraint(qcm.model, [d=1:max_depth], sum(gates_real[(2*I-1),(2*J), n_g] * z_bin_var[n_g,d] for n_g = 1:num_gates) - intercept <= 0)
+                                else           
+                                    JuMP.@constraint(qcm.model, [d=1:max_depth], sum(gates_real[(2*I-1),(2*J), n_g] * z_bin_var[n_g,d] for n_g = 1:num_gates) 
                                                                         - slope*(sum(gates_real[(2*I-1),(2*J-1), n_g] * z_bin_var[n_g,d] for n_g = 1:num_gates)) - intercept <= 0)
                                 end
                                 num_facets += 1
@@ -406,11 +427,9 @@ function constraint_convex_hull_complex_gates(qcm::QuantumCircuitModel)
                             elseif v3[2] - slope*v3[1] - intercept >= 1E-6
 
                                 if isapprox(abs(slope), 0, atol=1E-6)
-                                    
-                                    JuMP.@constraint(qcm.model, [d=1:depth], sum(gates_real[(2*I-1),(2*J), n_g] * z_bin_var[n_g,d] for n_g = 1:num_gates) - intercept >= 0)
+                                    JuMP.@constraint(qcm.model, [d=1:max_depth], sum(gates_real[(2*I-1),(2*J), n_g] * z_bin_var[n_g,d] for n_g = 1:num_gates) - intercept >= 0)
                                 else
-                                    
-                                    JuMP.@constraint(qcm.model, [d=1:depth], sum(gates_real[(2*I-1),(2*J), n_g] * z_bin_var[n_g,d] for n_g = 1:num_gates) 
+                                    JuMP.@constraint(qcm.model, [d=1:max_depth], sum(gates_real[(2*I-1),(2*J), n_g] * z_bin_var[n_g,d] for n_g = 1:num_gates) 
                                                                         - slope*(sum(gates_real[(2*I-1),(2*J-1), n_g] * z_bin_var[n_g,d] for n_g = 1:num_gates)) - intercept >= 0)
                                 end
                                 num_facets += 1
@@ -422,11 +441,9 @@ function constraint_convex_hull_complex_gates(qcm::QuantumCircuitModel)
                         else isinf(slope)
 
                             if v3[1] >= v1[1] + 1E-6
-                                
-                                JuMP.@constraint(qcm.model, [d=1:depth], sum(gates_real[(2*I-1),(2*J-1), n_g] * z_bin_var[n_g,d] for n_g = 1:num_gates) >= v1[1])
+                                JuMP.@constraint(qcm.model, [d=1:max_depth], sum(gates_real[(2*I-1),(2*J-1), n_g] * z_bin_var[n_g,d] for n_g = 1:num_gates) >= v1[1])
                             elseif v3[1] <= v1[1] - 1E-6
-                                
-                                JuMP.@constraint(qcm.model, [d=1:depth], sum(gates_real[(2*I-1),(2*J-1), n_g] * z_bin_var[n_g,d] for n_g = 1:num_gates) <= v1[1])
+                                JuMP.@constraint(qcm.model, [d=1:max_depth], sum(gates_real[(2*I-1),(2*J-1), n_g] * z_bin_var[n_g,d] for n_g = 1:num_gates) <= v1[1])
                             else
                                 Memento.warn(_LOGGER, "Indeterminate direction for the planar-hull cut")
                             end
@@ -439,7 +456,7 @@ function constraint_convex_hull_complex_gates(qcm::QuantumCircuitModel)
         end
         
         if num_facets > 0
-            Memento.info(_LOGGER, "Applied $num_facets planar-hull cuts per depth of the decomposition")
+            Memento.info(_LOGGER, "Applied $num_facets planar-hull cuts per max_depth of the decomposition")
         end
         
     end
@@ -447,22 +464,8 @@ function constraint_convex_hull_complex_gates(qcm::QuantumCircuitModel)
     return
 end
 
-function constraint_complex_unit_magnitude(qcm::QuantumCircuitModel)
-
-    depth = qcm.data["maximum_depth"]
-    n_r   = size(qcm.data["gates_real"])[1]
-    n_c   = size(qcm.data["gates_real"])[2]
-
-    JuMP.@constraint(qcm.model, [i=1:2:n_r, j=1:2:n_c, d=1:(depth-1)],  qcm.variables[:U_var][i,j,d] + qcm.variables[:U_var][i,j+1,d] <= sqrt(2))
-    JuMP.@constraint(qcm.model, [i=1:2:n_r, j=1:2:n_c, d=1:(depth-1)], -qcm.variables[:U_var][i,j,d] + qcm.variables[:U_var][i,j+1,d] <= sqrt(2))         
-    JuMP.@constraint(qcm.model, [i=1:2:n_r, j=1:2:n_c, d=1:(depth-1)],  qcm.variables[:U_var][i,j,d] - qcm.variables[:U_var][i,j+1,d] <= sqrt(2))
-    JuMP.@constraint(qcm.model, [i=1:2:n_r, j=1:2:n_c, d=1:(depth-1)], -qcm.variables[:U_var][i,j,d] - qcm.variables[:U_var][i,j+1,d] <= sqrt(2))
-
-    return
-end
-
 function constraint_unitary_property(qcm::QuantumCircuitModel)
-    depth      = qcm.data["maximum_depth"]
+    max_depth  = qcm.data["maximum_depth"]
     Z_var      = qcm.variables[:Z_var]
     z_bin_var  = qcm.variables[:z_bin_var]
     gates_real = qcm.data["gates_real"]
@@ -482,7 +485,7 @@ function constraint_unitary_property(qcm::QuantumCircuitModel)
         end
     end
 
-    for d = 1:depth
+    for d = 1:max_depth
         for i=1:num_gates
             for j=i:num_gates
                 if i == j 
@@ -500,7 +503,7 @@ function constraint_unitary_property(qcm::QuantumCircuitModel)
     end
 
     # Unitary constraint
-    JuMP.@constraint(qcm.model, [d=1:depth], sum(Z_var[i,j,d]* QCO.round_real_value.(gates_real[:,:,i] * gates_adjoint[:,:,j] + 
+    JuMP.@constraint(qcm.model, [d=1:max_depth], sum(Z_var[i,j,d]* QCO.round_real_value.(gates_real[:,:,i] * gates_adjoint[:,:,j] + 
                                                  gates_real[:,:,j] * gates_adjoint[:,:,i]) for i=1:(num_gates-1), 
                                                  j=(i+1):num_gates) .== zeros(n_r, n_c))
 
diff --git a/src/data.jl b/src/data.jl
index 9bc130e..9d0a72f 100644
--- a/src/data.jl
+++ b/src/data.jl
@@ -633,8 +633,11 @@ function get_full_sized_kron_symbol_gate(input::String, num_qubits::Int64)
 
     kron_gates = QCO._parse_gates_with_kron_symbol(input)
     
-    M = 1
+    if length(unique(kron_gates)) !== length(kron_gates)
+        Memento.error(_LOGGER, "Specify only a single gate per qubit within kron symbol gate $input")
+    end
 
+    M = 1
     for i = 1:length(kron_gates)
         
         gate_type, qubit_loc = QCO._parse_gate_string(kron_gates[i], type = true, qubits = true)
@@ -722,9 +725,11 @@ function _catch_input_gate_errors(gate_type::String, qubit_loc::Vector{Int64}, n
     end
     
     if (gate_type in QCO.TWO_QUBIT_GATES) && (length(qubit_loc) != 2)
-        Memento.error(_LOGGER, "Specify two qubits for $input_gate, which is a 2-qubit gate")
-    elseif (gate_type in QCO.ONE_QUBIT_GATES) && (length(qubit_loc) != 1)
-        Memento.error(_LOGGER, "Specify only one qubit for $input_gate, which is a 1-qubit gate")
+        Memento.error(_LOGGER, "Specify two qubits for 2-qubit gate $gate_type in elementary gates")
+    elseif (gate_type in QCO.ONE_QUBIT_GATES) && (length(qubit_loc) == 0)
+        Memento.error(_LOGGER, "Specify a qubit for the 1-qubit gate $gate_type in elementary gates")
+    elseif (gate_type in QCO.ONE_QUBIT_GATES) && (length(qubit_loc) >= 2)
+        Memento.error(_LOGGER, "Specify only one qubit for the 1-qubit gate $gate_type in elementary gates")
     end
 
     if isempty(qubit_loc) && (gate_type !== "GR")
diff --git a/src/log.jl b/src/log.jl
index b94c1b5..7635cf3 100644
--- a/src/log.jl
+++ b/src/log.jl
@@ -6,12 +6,15 @@ this function aids in visualizing the optimal circuit decomposition.
 """
 function visualize_solution(results::Dict{String, Any}, data::Dict{String, Any}; gate_sequence = false)
 
-    if results["primal_status"] != MOI.FEASIBLE_POINT 
+    _header_color = :cyan 
+    _main_color   = :White
+
+    if !(results["primal_status"] in [MOI.FEASIBLE_POINT, MOI.NEARLY_FEASIBLE_POINT])  
         
-        if results["termination_status"] != MOI.TIME_LIMIT
-            Memento.warn(_LOGGER, "Infeasible primal status. Gate decomposition may be inaccurate")
-        else 
+        if results["termination_status"]  == MOI.TIME_LIMIT 
             Memento.warn(_LOGGER, "Optimizer hits time limit with an infeasible primal status. Gate decomposition may be inaccurate")
+        else
+            Memento.warn(_LOGGER, "Infeasible primal status. Gate decomposition may be inaccurate")
         end
 
         return
@@ -21,58 +24,57 @@ function visualize_solution(results::Dict{String, Any}, data::Dict{String, Any};
 
     if !isempty(gates_sol_compressed)
 
-        printstyled("\n","=============================================================================","\n"; color = :cyan)
-        
-        printstyled("Quantum Circuit Model Data","\n"; color = :red)
+        printstyled("\n","=============================================================================","\n"; color = _main_color)
+        printstyled("QuantumCircuitOpt version: ", _QCOPT_VERSION, "\n"; color = _header_color)
+
+        printstyled("\n","Quantum Circuit Model Data","\n"; color = _header_color)
         
-        printstyled("\n","  ","Number of qubits: ", data["num_qubits"], "\n"; color = :cyan)
+        printstyled("\n","  ","Number of qubits: ", data["num_qubits"], "\n"; color = _main_color)
         
-        printstyled("  ","Total number of elementary gates (after presolve): ",size(data["gates_real"])[3],"\n"; color = :cyan)
+        printstyled("  ","Total number of elementary gates (after presolve): ",size(data["gates_real"])[3],"\n"; color = _main_color)
         
-        printstyled("  ","Maximum depth of decomposition: ", data["maximum_depth"],"\n"; color = :cyan)
+        printstyled("  ","Maximum depth of decomposition: ", data["maximum_depth"],"\n"; color = _main_color)
         
-        printstyled("  ","Input elementary gates: ", data["elementary_gates"],"\n"; color = :cyan)
+        printstyled("  ","Elementary gates: ", data["elementary_gates"],"\n"; color = _main_color)
 
         if "discretization" in keys(data)
             for i in keys(data["discretization"])
-                printstyled("    ","$i discretization: ", ceil.(rad2deg.(data["discretization"][i]), digits = 1),"\n"; color = :cyan)
+                printstyled("    ","$i discretization: ", ceil.(rad2deg.(data["discretization"][i]), digits = 1),"\n"; color = _main_color)
             end
         end
-        
-        # printstyled("  ","Input target gate: ", data["target_gate"]["type"],"\n"; color = :cyan)
-        
-        printstyled("  ","Type of decomposition: ", data["decomposition_type"],"\n"; color = :cyan)
+                
+        printstyled("  ","Type of decomposition: ", data["decomposition_type"],"\n"; color = _main_color)
 
-        printstyled("  ","MIP optimizer: ", results["optimizer"],"\n"; color = :cyan)
+        printstyled("  ","MIP optimizer: ", results["optimizer"],"\n"; color = _main_color)
 
-        printstyled("\n","Optimal Circuit Decomposition","\n","\n"; color = :red)
+        printstyled("\n","Optimal Circuit Decomposition","\n","\n"; color = _header_color)
         
         print("  ")
         
         for i=1:length(gates_sol_compressed)
             
             if i != length(gates_sol_compressed)
-                printstyled(gates_sol_compressed[i], " * "; color = :cyan)
+                printstyled(gates_sol_compressed[i], " * "; color = _main_color)
             else    
                 if data["decomposition_type"] in ["exact_optimal", "exact_feasible"]
-                    printstyled(gates_sol_compressed[i], " = ", "Target gate","\n"; color = :cyan)
+                    printstyled(gates_sol_compressed[i], " = ", "Target gate","\n"; color = _main_color)
                 elseif data["decomposition_type"] == "approximate"
-                    printstyled(gates_sol_compressed[i], " ≈ ", "Target gate","\n"; color = :cyan)
+                    printstyled(gates_sol_compressed[i], " ≈ ", "Target gate","\n"; color = _main_color)
                 end
             end
 
         end
 
         if data["decomposition_type"] == "approximate"
-            printstyled("  ","||Decomposition error||₂: ", round(LA.norm(results["solution"]["slack_var"]), digits = 10),"\n"; color = :cyan)
+            printstyled("  ","||Decomposition error||₂: ", round(LA.norm(results["solution"]["slack_var"]), digits = 10),"\n"; color = _main_color)
         end
 
         if data["objective"] == "minimize_depth"
 
             if length(data["identity_idx"]) >= 1 && (data["decomposition_type"] !== "exact_feasible") && !(results["termination_status"] == MOI.TIME_LIMIT)
-                printstyled("  ","Minimum optimal depth: ", length(gates_sol_compressed),"\n"; color = :cyan)
+                printstyled("  ","Minimum optimal depth: ", length(gates_sol_compressed),"\n"; color = _main_color)
             else 
-                printstyled("  ","Decomposition depth: ", length(gates_sol_compressed),"\n"; color = :cyan)
+                printstyled("  ","Decomposition depth: ", length(gates_sol_compressed),"\n"; color = _main_color)
             end
 
         elseif data["objective"] == "minimize_cnot"
@@ -80,24 +82,23 @@ function visualize_solution(results::Dict{String, Any}, data::Dict{String, Any};
             if !isempty(data["cnot_idx"])
                 
                 if data["decomposition_type"] in ["exact_optimal", "exact_feasible"]
-                    printstyled("  ","Minimum number of CNOT gates: ", round(results["objective"], digits = 6),"\n"; color = :cyan)
+                    printstyled("  ","Minimum number of CNOT gates: ", round(results["objective"], digits = 6),"\n"; color = _main_color)
                 
                 elseif data["decomposition_type"] == "approximate"
-                    printstyled("  ","Minimum number of CNOT gates: ", round((results["objective"] - results["objective_slack_penalty"]*LA.norm(results["solution"]["slack_var"])^2), digits = 6),"\n"; color = :cyan)
+                    printstyled("  ","Minimum number of CNOT gates: ", round((results["objective"] - results["objective_slack_penalty"]*LA.norm(results["solution"]["slack_var"])^2), digits = 6),"\n"; color = _main_color)
                 end
             
             end
 
         end
 
-        printstyled("  ","Optimizer run time: ", ceil(results["solve_time"], digits=2)," sec.","\n"; color = :cyan)
+        printstyled("  ","Optimizer run time: ", ceil(results["solve_time"], digits=2)," sec.","\n"; color = _main_color)
             
         if results["termination_status"] == MOI.TIME_LIMIT
-            printstyled("  ","Termination status: MOI.TIME_LIMIT", "\n"; color = :cyan)
-            printstyled("  ","Decomposition may not be optimal", "\n"; color = :cyan)
+            printstyled("  ","Termination status: TIME_LIMIT", "\n"; color = _main_color)
         end
 
-        printstyled("=============================================================================","\n"; color = :cyan)      
+        printstyled("=============================================================================","\n"; color = _main_color)      
 
     else
         Memento.warn(_LOGGER, "Valid integral feasible solutions could not be found to visualize the solution")
diff --git a/src/objective.jl b/src/objective.jl
index ba17fb7..faec2a8 100644
--- a/src/objective.jl
+++ b/src/objective.jl
@@ -4,9 +4,7 @@
 
 function objective_minimize_total_depth(qcm::QuantumCircuitModel)
     
-    n_r          = size(qcm.data["gates_real"])[1]
-    n_c          = size(qcm.data["gates_real"])[2]
-    depth        = qcm.data["maximum_depth"]
+    max_depth    = qcm.data["maximum_depth"]
     identity_idx = qcm.data["identity_idx"]
     num_gates    = size(qcm.data["gates_real"])[3]
 
@@ -14,13 +12,14 @@ function objective_minimize_total_depth(qcm::QuantumCircuitModel)
 
     if !isempty(identity_idx)
         if decomposition_type == "exact_optimal"
-            JuMP.@objective(qcm.model, Min, sum(qcm.variables[:z_bin_var][n,d] for n = 1:num_gates, d=1:depth if !(n in identity_idx)))
+            JuMP.@objective(qcm.model, Min, sum(qcm.variables[:z_bin_var][n,d] for n = 1:num_gates, d=1:max_depth if !(n in identity_idx)))
 
         elseif decomposition_type == "exact_feasible"
             QCO.objective_feasibility(qcm)
 
         elseif decomposition_type == "approximate"
-            JuMP.@objective(qcm.model, Min, sum(qcm.variables[:z_bin_var][n,d] for n = 1:num_gates, d=1:depth if !(n in identity_idx)) + qcm.options.objective_slack_penalty * sum(qcm.variables[:slack_var][i,j]^2 for i=1:n_r, j=1:n_c))
+            JuMP.@objective(qcm.model, Min, sum(qcm.variables[:z_bin_var][n,d] for n = 1:num_gates, d=1:max_depth if !(n in identity_idx)) 
+            + qcm.options.objective_slack_penalty * sum(qcm.variables[:slack_var_oa]))
         end
     else
         QCO.objective_feasibility(qcm)
@@ -37,9 +36,7 @@ function objective_feasibility(qcm::QuantumCircuitModel)
     if decomposition_type in ["exact_optimal", "exact_feasible"]
         # Feasibility objective
     elseif decomposition_type == "approximate"
-        n_r     = size(qcm.data["gates_real"])[1]
-        n_c     = size(qcm.data["gates_real"])[2]
-        JuMP.@objective(qcm.model, Min, sum(qcm.variables[:slack_var][i,j]^2 for i=1:n_r, j=1:n_c))
+        JuMP.@objective(qcm.model, Min, sum(qcm.variables[:slack_var_oa]))
     end 
 
     return
@@ -47,20 +44,19 @@ end
 
 function objective_minimize_cnot_gates(qcm::QuantumCircuitModel)
     
-    n_r      = size(qcm.data["gates_real"])[1]
-    n_c      = size(qcm.data["gates_real"])[2]
-    depth    = qcm.data["maximum_depth"] 
-    cnot_idx = qcm.data["cnot_idx"]
+    max_depth = qcm.data["maximum_depth"] 
+    cnot_idx  = qcm.data["cnot_idx"]
 
     decomposition_type = qcm.data["decomposition_type"]
     
     if !isempty(qcm.data["cnot_idx"])
 
         if decomposition_type in ["exact_optimal", "exact_feasible"]
-            JuMP.@objective(qcm.model, Min, sum(qcm.variables[:z_bin_var][n,d] for n in cnot_idx, d=1:depth))
+            JuMP.@objective(qcm.model, Min, sum(qcm.variables[:z_bin_var][n,d] for n in cnot_idx, d=1:max_depth))
             
         elseif decomposition_type == "approximate"            
-            JuMP.@objective(qcm.model, Min, sum(qcm.variables[:z_bin_var][n,d] for n in cnot_idx, d=1:depth) + (qcm.options.objective_slack_penalty * sum(qcm.variables[:slack_var][i,j]^2 for i=1:n_r, j=1:n_c)))
+            JuMP.@objective(qcm.model, Min, sum(qcm.variables[:z_bin_var][n,d] for n in cnot_idx, d=1:max_depth) 
+            + (qcm.options.objective_slack_penalty * sum(qcm.variables[:slack_var_oa])))
 
         end
     else
diff --git a/src/qc_model.jl b/src/qc_model.jl
index 9516909..eb866c4 100644
--- a/src/qc_model.jl
+++ b/src/qc_model.jl
@@ -44,6 +44,7 @@ function variable_QCModel_balas(qcm::QuantumCircuitModel)
 
     if qcm.data["decomposition_type"] == "approximate"
         QCO.variable_slack_for_feasibility(qcm)
+        QCO.variable_slack_var_outer_approximation(qcm)
     end
 
     QCO.variable_QCModel_valid(qcm)
@@ -60,6 +61,7 @@ function constraint_QCModel_balas(qcm::QuantumCircuitModel)
     QCO.constraint_gate_product_linearization(qcm)
     QCO.constraint_gate_target_condition(qcm)
     QCO.constraint_cnot_gate_bounds(qcm)
+    (qcm.data["decomposition_type"] == "approximate") && (QCO.constraint_slack_var_outer_approximation(qcm))
     (!qcm.data["are_gates_real"]) && (QCO.constraint_complex_to_real_symmetry(qcm))
 
     QCO.constraint_QCModel_valid(qcm)
@@ -75,6 +77,7 @@ function variable_QCModel_compact(qcm::QuantumCircuitModel)
 
     if qcm.data["decomposition_type"] == "approximate"
         QCO.variable_slack_for_feasibility(qcm)
+        QCO.variable_slack_var_outer_approximation(qcm)
     end
 
     QCO.variable_QCModel_valid(qcm)
@@ -106,7 +109,8 @@ function constraint_QCModel_compact(qcm::QuantumCircuitModel)
     QCO.constraint_gate_product_linearization(qcm)
     QCO.constraint_gate_target_condition_compact(qcm)
     QCO.constraint_cnot_gate_bounds(qcm)
-    (!qcm.data["are_gates_real"]) && (QCO.constraint_complex_to_real_symmetry(qcm))
+    (qcm.data["decomposition_type"] == "approximate") && (QCO.constraint_slack_var_outer_approximation(qcm))
+    # (!qcm.data["are_gates_real"]) && (QCO.constraint_complex_to_real_symmetry(qcm))
 
     QCO.constraint_QCModel_valid(qcm)
 
@@ -124,7 +128,6 @@ function constraint_QCModel_valid(qcm::QuantumCircuitModel)
             QCO.constraint_idempotent_gates(qcm)
             QCO.constraint_identity_gate_symmetry(qcm)
             QCO.constraint_convex_hull_complex_gates(qcm)
-            QCO.constraint_complex_unit_magnitude(qcm)
             QCO.constraint_unitary_property(qcm)
             
         elseif qcm.options.all_valid_constraints == 0 
@@ -134,7 +137,6 @@ function constraint_QCModel_valid(qcm::QuantumCircuitModel)
             qcm.options.idempotent_gate_constraints         && QCO.constraint_idempotent_gates(qcm)
             qcm.options.identity_gate_symmetry_constraints  && QCO.constraint_identity_gate_symmetry(qcm)
             qcm.options.convex_hull_gate_constraints        && QCO.constraint_convex_hull_complex_gates(qcm)
-            qcm.options.unit_magnitude_constraints          && QCO.constraint_complex_unit_magnitude(qcm)
             qcm.options.unitary_constraints                 && QCO.constraint_unitary_property(qcm)
         end
 
diff --git a/src/solution.jl b/src/solution.jl
index ff5de05..bbde968 100644
--- a/src/solution.jl
+++ b/src/solution.jl
@@ -11,7 +11,7 @@ function build_QCModel_result(qcm::QuantumCircuitModel, solve_time::Number)
     solution = Dict{String,Any}()
 
     if result_count > 0
-        solution = build_QCModel_solution(qcm)
+        solution = QCO.build_QCModel_solution(qcm)
     else
         Memento.warn(_LOGGER, "Quantum circuit model has no results, solution cannot be built")
     end
@@ -20,8 +20,8 @@ function build_QCModel_result(qcm::QuantumCircuitModel, solve_time::Number)
         "optimizer" => JuMP.solver_name(qcm.model),
         "termination_status" => JuMP.termination_status(qcm.model),
         "primal_status" => JuMP.primal_status(qcm.model),
-        "objective" => get_objective_value(qcm.model),
-        "objective_lb" => get_objective_bound(qcm.model),
+        "objective" => QCO.get_objective_value(qcm.model),
+        "objective_lb" => QCO.get_objective_bound(qcm.model),
         "solve_time" => solve_time,
         "solution" => solution,
     )
@@ -34,7 +34,6 @@ function build_QCModel_result(qcm::QuantumCircuitModel, solve_time::Number)
     return result
 end
 
-
 ""
 function get_objective_value(model::JuMP.Model)
     obj_val = NaN
diff --git a/src/types.jl b/src/types.jl
index c365bec..d510d40 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -13,10 +13,10 @@ mutable struct QCModelOptions
     involutory_gate_constraints        :: Bool
     redundant_gate_pair_constraints    :: Bool 
     identity_gate_symmetry_constraints :: Bool
+    fix_unitary_variables              :: Bool
     visualize_solution                 :: Bool
     idempotent_gate_constraints        :: Bool
     convex_hull_gate_constraints       :: Bool
-    unit_magnitude_constraints         :: Bool
     unitary_constraints                :: Bool
 
     time_limit                         :: Float64
@@ -30,18 +30,18 @@ end
 This function returns the default options for building the struct `QCModelOptions`.
 """
 function get_default_options()
-    model_type                         = "compact_formulation" # compact_formulation, balas_formulation
+    model_type                         = "compact_formulation" # check MODEL_TYPES in src/qc_model.jl for options
     
     all_valid_constraints              = 0    # -1, 0, 1
     commute_gate_constraints           = true # true, false
     involutory_gate_constraints        = true # true, false
     redundant_gate_pair_constraints    = true # true, false
     identity_gate_symmetry_constraints = true # true, false
+    fix_unitary_variables              = true # true, false
     visualize_solution                 = true # true, false
 
     idempotent_gate_constraints        = false # true, false
     convex_hull_gate_constraints       = false # true, false
-    unit_magnitude_constraints         = false # true, false
     unitary_constraints                = false # true, false
 
     time_limit                         = 10800 # float value
@@ -55,10 +55,10 @@ function get_default_options()
                           involutory_gate_constraints,
                           redundant_gate_pair_constraints,
                           identity_gate_symmetry_constraints,
+                          fix_unitary_variables,
                           visualize_solution,
                           idempotent_gate_constraints,
                           convex_hull_gate_constraints,
-                          unit_magnitude_constraints,
                           unitary_constraints,
                           time_limit,
                           relax_integrality,
diff --git a/src/utility.jl b/src/utility.jl
index 561fba9..8fb69a5 100644
--- a/src/utility.jl
+++ b/src/utility.jl
@@ -42,17 +42,14 @@ function gate_element_bounds(M::Array{Float64,3})
     M_l = zeros(size(M)[1], size(M)[2])
     M_u = zeros(size(M)[1], size(M)[2])
 
-    for i = 1:size(M)[1]
-        for j = 1:size(M)[2]
-
-            M_l[i,j] = minimum(M[i,j,:])
-            M_u[i,j] = maximum(M[i,j,:])
-
-            if M_l[i,j] > M_u[i,j]
-                Memento.error(_LOGGER, "Lower and upper bound conflict in the elements of input elementary gates")
-            elseif (M_l[i,j] in [-Inf, Inf]) || (M_u[i,j] in [-Inf, Inf])
-                Memento.error(_LOGGER, "Unbounded entry detected in the input elementary gates")
-            end
+    for i = 1:size(M)[1], j = 1:size(M)[2]
+        M_l[i,j] = minimum(M[i,j,:])
+        M_u[i,j] = maximum(M[i,j,:])
+
+        if M_l[i,j] > M_u[i,j]
+            Memento.error(_LOGGER, "Lower and upper bound conflict in the elements of input elementary gates")
+        elseif (M_l[i,j] in [-Inf, Inf]) || (M_u[i,j] in [-Inf, Inf])
+            Memento.error(_LOGGER, "Unbounded entry detected in the input elementary gates")
         end
     end
 
@@ -71,28 +68,25 @@ function get_commutative_gate_pairs(M::Dict{String,Any}; identity_in_pairs = tru
     commute_pairs = Array{Tuple{Int64,Int64},1}()
     commute_pairs_prodIdentity = Array{Tuple{Int64,Int64},1}()
 
-    for i = 1:(num_gates-1)
-        for j = (i+1):num_gates
-            M_i = M["$i"]["matrix"]
-            M_j = M["$j"]["matrix"]
-
-            if ("Identity" in M["$i"]["type"]) || ("Identity" in M["$j"]["type"])
-                continue
-            end
-            
-            M_ij = M_i*M_j
-            M_ji = M_j*M_i
-            Id = Matrix(LA.I, size(M_ij)[1], size(M_ij)[2])
+    for i = 1:(num_gates-1), j = (i+1):num_gates
+        M_i = M["$i"]["matrix"]
+        M_j = M["$j"]["matrix"]
 
-            # Commuting pairs == Identity 
-            if isapprox(M_ij, Id, atol=1E-6)
-                push!(commute_pairs_prodIdentity, (i,j))
-            
-            # Commuting pairs != Identity 
-            elseif isapprox(M_ij, M_ji, atol = 1E-4)
-                push!(commute_pairs, (i, j))
+        if ("Identity" in M["$i"]["type"]) || ("Identity" in M["$j"]["type"])
+            continue
+        end
+        
+        M_ij = M_i*M_j
+        M_ji = M_j*M_i
+        Id = Matrix(LA.I, size(M_ij)[1], size(M_ij)[2])
 
-            end
+        # Commuting pairs == Identity 
+        if isapprox(M_ij, Id, atol=1E-6)
+            push!(commute_pairs_prodIdentity, (i,j))
+        
+        # Commuting pairs != Identity 
+        elseif isapprox(M_ij, M_ji, atol = 1E-4)
+            push!(commute_pairs, (i, j))
 
         end
     end
@@ -133,29 +127,26 @@ function get_redundant_gate_product_pairs(M::Dict{String,Any})
     redundant_pairs_idx = Array{Tuple{Int64,Int64},1}()
 
     # Non-Identity redundant pairs
-    for i = 1:(num_gates-1)
-        for j = (i+1):num_gates
-            M_i = M["$i"]["matrix"]
-            M_j = M["$j"]["matrix"]
+    for i = 1:(num_gates-1), j = (i+1):num_gates
+        M_i = M["$i"]["matrix"]
+        M_j = M["$j"]["matrix"]
 
-            if ("Identity" in M["$i"]["type"]) || ("Identity" in M["$j"]["type"])
-                continue
-            end
+        if ("Identity" in M["$i"]["type"]) || ("Identity" in M["$j"]["type"])
+            continue
+        end
+        
+        for k = 1:num_gates 
+            if (k != i) && (k != j) && !("Identity" in M["$k"]["type"])                
+                M_k = M["$k"]["matrix"]
             
-            for k = 1:num_gates 
-                if (k != i) && (k != j) && !("Identity" in M["$k"]["type"])                
-                    M_k = M["$k"]["matrix"]
-                
-                    if isapprox(M_i*M_j, M_k, atol = 1E-4)
-                        push!(redundant_pairs_idx, (i, j))
-                    end
+                if isapprox(M_i*M_j, M_k, atol = 1E-4)
+                    push!(redundant_pairs_idx, (i, j))
+                end
 
-                    if isapprox(M_j*M_i, M_k, atol = 1E-4)
-                        push!(redundant_pairs_idx, (j, i))
-                    end
+                if isapprox(M_j*M_i, M_k, atol = 1E-4)
+                    push!(redundant_pairs_idx, (j, i))
                 end
             end
-
         end
     end
 
@@ -563,7 +554,7 @@ which the input gate is located. For example, if the input string is `CRX_2_3`,
 function _parse_gate_string(s::String; type=false, qubits=false)
 
     if occursin(kron_symbol, s)
-        Memento.error(_LOGGER, "Kron symbol is not supported for parsing qubit numbers in $s")
+        Memento.error(_LOGGER, "Kron symbol is not supported for parsing qubit locations in $s. Instead try `get_full_sized_kron_symbol_gate` function.")
     end
     
     gates = Vector{String}()
@@ -724,4 +715,101 @@ function _determinant_test_for_infeasibility(data::Dict{String,Any})
 
     end
 
+end
+
+"""
+    _get_elementary_gates_fixed_indices(M::Array{T,3} where T <: Number)
+
+Given the set of input elementary gates in real form, 
+this function returns a dictionary of tuples of indices wholse values are fixed in `sum_k (z_k*M[:,:,k])`. 
+"""
+function _get_elementary_gates_fixed_indices(M::Array{T,3} where T <: Number)
+    N = size(M[:,:,1])[1]
+    M_l, M_u = QCO.gate_element_bounds(M)
+
+    G_fixed_idx = Dict{Tuple{Int64, Int64}, Any}()
+    for i=1:N, j=1:N
+        if isapprox(M_l[i,j], M_u[i,j], atol = 1E-6)
+            G_fixed_idx[(i,j)] = Dict{String, Any}("value" => M_l[i,j])
+        end
+    end
+
+    return G_fixed_idx
+end
+
+"""
+    _get_unitary_variables_fixed_indices(M::Array{T,3} where T <: Number, 
+                                         maximum_depth::Int64)
+
+Given a 3D array of real square matrices (representing gates), and a maximum alowable depth, 
+this function returns a dictionary of tuples of indices wholse values are fixed in the unitary matrices for every depth 
+of the circuit. 
+"""
+function _get_unitary_variables_fixed_indices(M::Array{T,3} where T <: Number, maximum_depth::Int64)
+    
+    N = size(M)[1]
+    G_fixed_idx = QCO._get_elementary_gates_fixed_indices(M)
+
+    U_fixed_idx = Dict{Int64, Any}()
+    for depth = 1:(maximum_depth-1)
+        U_fixed_idx[depth] = Dict{Tuple{Int64, Int64}, Any}()
+        if depth == 1 
+            # Assuming data["initial_gate"] == "Identity"
+            U_fixed_idx[depth] = G_fixed_idx
+        else 
+            U_fixed_idx[depth] = QCO._get_matrix_product_fixed_indices(U_fixed_idx[depth-1], G_fixed_idx, N)
+        end
+    end
+    
+    return U_fixed_idx
+end
+
+"""
+    _get_matrix_product_fixed_indices(left_matrix_fixed_idx::Dict{Tuple{Int64, Int64}, Any}, 
+                                  right_matrix_fixed_idx::Dict{Tuple{Int64, Int64}, Any}, 
+                                  N::Int64)
+
+Given left and right square matrices of size `NxN`, in a dictionary format with tuples of indices whose values are fixed, 
+this function returns a dictionary of tuples of indices wholse values are fixed in `left_matrix * right_matrix`. 
+"""
+function _get_matrix_product_fixed_indices(left_matrix_fixed_idx::Dict{Tuple{Int64, Int64}, Any}, 
+                                           right_matrix_fixed_idx::Dict{Tuple{Int64, Int64}, Any}, 
+                                           N::Int64)
+    
+    product_fixed_idx = Dict{Tuple{Int64, Int64}, Any}()
+
+    for row = 1:N
+        left_matrix_constants_row = filter(p -> (p.first[1] == row), left_matrix_fixed_idx)
+        left_matrix_zeros_row = filter(p -> (isapprox(p.second["value"], 0, atol = 1E-6)), left_matrix_constants_row)
+
+        v_constants_row = keys(left_matrix_constants_row) |> collect .|> last
+        v_zeros_row = keys(left_matrix_zeros_row) |> collect .|> last
+    
+        for col = 1:N
+            right_matrix_constants_col = filter(p -> (p.first[2] == col), right_matrix_fixed_idx)
+            right_matrix_zeros_col = filter(p -> (isapprox(p.second["value"], 0, atol=1E-6)), right_matrix_constants_col)
+
+            v_constants_col = keys(right_matrix_constants_col) |> collect .|> first
+            v_zeros_col = keys(right_matrix_zeros_col) |> collect .|> first
+
+            all_zero_idx = union(v_zeros_row, v_zeros_col)
+            
+            # Keep track of zero value indices in matrix product
+            if sort(all_zero_idx) == 1:N
+                product_fixed_idx[(row,col)] = Dict{String, Any}("value" => 0)
+            end
+
+            # Keep track of non-zero constant value indices in matrix product
+            if (sort(v_constants_row) == 1:N) && (sort(v_constants_col) == 1:N)
+                value = 0
+                for i = 1:N
+                    value += left_matrix_constants_row[(row, i)]["value"] * right_matrix_constants_col[(i, col)]["value"]
+                end
+                product_fixed_idx[(row,col)] = Dict{String, Any}("value" => value)
+            end
+
+        end 
+    end
+
+    return product_fixed_idx
 end
\ No newline at end of file
diff --git a/src/variables.jl b/src/variables.jl
index b5a2d7c..cf1596d 100644
--- a/src/variables.jl
+++ b/src/variables.jl
@@ -7,24 +7,20 @@ function variable_gates_per_depth(qcm::QuantumCircuitModel)
     tol_0 = 1E-6
     n_r     = size(qcm.data["gates_real"])[1]
     n_c     = size(qcm.data["gates_real"])[2]
-    depth   = qcm.data["maximum_depth"]
+    max_depth   = qcm.data["maximum_depth"]
 
     M_l, M_u = QCO.gate_element_bounds(qcm.data["gates_real"])
 
-    qcm.variables[:G_var] = JuMP.@variable(qcm.model, M_l[i,j] <= G_var[i=1:n_r, j=1:n_c, 1:depth] <= M_u[i,j])
+    qcm.variables[:G_var] = JuMP.@variable(qcm.model, M_l[i,j] <= G_var[i=1:n_r, j=1:n_c, 1:max_depth] <= M_u[i,j])
 
     num_vars_fixed = 0
 
-    for i=1:n_r
-        for j=1:n_c
-
-            if isapprox(M_l[i,j], M_u[i,j], atol=tol_0)
-                for d=1:depth
-                    JuMP.fix(G_var[i,j,d], M_l[i,j]; force=true)
-                end
-                num_vars_fixed += 1
+    for i=1:n_r, j=1:n_c
+        if isapprox(M_l[i,j], M_u[i,j], atol=tol_0)
+            for d=1:max_depth
+                JuMP.fix(G_var[i,j,d], M_l[i,j]; force=true)
             end
-
+            num_vars_fixed += 1
         end
     end
 
@@ -37,9 +33,9 @@ end
 
 function variable_gates_onoff(qcm::QuantumCircuitModel)
     num_gates = size(qcm.data["gates_real"])[3]
-    depth   = qcm.data["maximum_depth"]
+    max_depth   = qcm.data["maximum_depth"]
 
-    qcm.variables[:z_bin_var] = JuMP.@variable(qcm.model, z_bin_var[1:num_gates,1:depth], Bin)
+    qcm.variables[:z_bin_var] = JuMP.@variable(qcm.model, z_bin_var[1:num_gates,1:max_depth], Bin)
 
     if "input_circuit" in keys(qcm.data)
         
@@ -62,22 +58,58 @@ function variable_gates_onoff(qcm::QuantumCircuitModel)
 end
 
 function variable_sequential_gate_products(qcm::QuantumCircuitModel)
-    depth  = qcm.data["maximum_depth"]
+    max_depth  = qcm.data["maximum_depth"]
     n_r    = size(qcm.data["gates_real"])[1]
     n_c    = size(qcm.data["gates_real"])[2]
 
-    qcm.variables[:U_var] = JuMP.@variable(qcm.model, -1 <= U_var[1:n_r, 1:n_c, 1:(depth-1)] <= 1)
+    if !qcm.options.fix_unitary_variables
+        qcm.variables[:U_var] = JuMP.@variable(qcm.model, -1 <= U_var[1:n_r, 1:n_c, 1:max_depth] <= 1)
+        return 
+    end
+
+    qcm.variables[:U_var] = JuMP.@variable(qcm.model, U_var[1:n_r, 1:n_c, 1:max_depth])
+    U_fixed_idx = QCO._get_unitary_variables_fixed_indices(qcm.data["gates_real"], max_depth)
+    # U_var_bound_tol = 1E-8
+
+    for depth = 1:(max_depth-1)
+        U_fix = U_fixed_idx[depth]
+        for ii = 1:n_r, jj = 1:n_c
+            if ((ii, jj) in keys(U_fix))
+                if isapprox(U_fix[(ii, jj)]["value"], -1, atol = 1E-6)
+                    JuMP.set_lower_bound(U_var[ii, jj, depth], -1)
+                    JuMP.set_upper_bound(U_var[ii, jj, depth], -1)
+                elseif isapprox(U_fix[(ii, jj)]["value"], 0, atol = 1E-6)
+                    JuMP.set_lower_bound(U_var[ii, jj, depth], 0)
+                    JuMP.set_upper_bound(U_var[ii, jj, depth], 0)
+                elseif isapprox(U_fix[(ii, jj)]["value"], 1, atol = 1E-6)
+                    JuMP.set_lower_bound(U_var[ii, jj, depth], 1)
+                    JuMP.set_upper_bound(U_var[ii, jj, depth], 1)
+                else 
+                    JuMP.set_lower_bound(U_var[ii, jj, depth], U_fix[(ii, jj)]["value"])
+                    JuMP.set_upper_bound(U_var[ii, jj, depth], U_fix[(ii, jj)]["value"])
+                end
+            else
+                JuMP.set_lower_bound(U_var[ii, jj, depth], -1)
+                JuMP.set_upper_bound(U_var[ii, jj, depth],  1)
+            end
+        end
+    end
+
+    for ii = 1:n_r, jj = 1:n_c
+        JuMP.set_lower_bound(U_var[ii, jj, max_depth], -1)
+        JuMP.set_upper_bound(U_var[ii, jj, max_depth],  1)
+    end
     
     return
 end
 
 function variable_gate_products_copy(qcm::QuantumCircuitModel)
-    depth   = qcm.data["maximum_depth"]
+    max_depth   = qcm.data["maximum_depth"]
     n_r     = size(qcm.data["gates_real"])[1]
     n_c     = size(qcm.data["gates_real"])[2]
     num_gates = size(qcm.data["gates_real"])[3]
 
-    qcm.variables[:V_var] = JuMP.@variable(qcm.model, -1 <= V_var[1:n_r, 1:n_c, 1:num_gates, 1:depth] <= 1)
+    qcm.variables[:V_var] = JuMP.@variable(qcm.model, -1 <= V_var[1:n_r, 1:n_c, 1:num_gates, 1:max_depth] <= 1)
     
     return
 end
@@ -85,10 +117,29 @@ end
 function variable_gate_products_linearization(qcm::QuantumCircuitModel)
     n_r     = size(qcm.data["gates_real"])[1]
     n_c     = size(qcm.data["gates_real"])[2]
-    depth   = qcm.data["maximum_depth"]
+    max_depth   = qcm.data["maximum_depth"]
     num_gates = size(qcm.data["gates_real"])[3]
 
-    qcm.variables[:zU_var] = JuMP.@variable(qcm.model, -1 <= zU_var[1:n_r, 1:n_c, 1:num_gates, 1:(depth-1)] <= 1)
+    if !qcm.options.fix_unitary_variables
+        qcm.variables[:zU_var] = JuMP.@variable(qcm.model, -1 <= zU_var[1:n_r, 1:n_c, 1:num_gates, 1:(max_depth-1)] <= 1)
+        return 
+    end
+
+    U_var = qcm.variables[:U_var]
+    qcm.variables[:zU_var] = JuMP.@variable(qcm.model, zU_var[1:n_r, 1:n_c, 1:num_gates, 1:(max_depth-1)])
+
+    for d = 1:(max_depth-1), i = 1:n_r, j = 1:n_c, k = 1:num_gates
+        U_var_l = JuMP.lower_bound(U_var[i,j,d])
+        U_var_u = JuMP.upper_bound(U_var[i,j,d])
+
+        if isapprox(U_var_l, 0, atol = 1E-6) && isapprox(U_var_u, 0, atol = 1E-6)
+            JuMP.set_lower_bound(zU_var[i,j,k,d], 0)
+            JuMP.set_upper_bound(zU_var[i,j,k,d], 0)
+        else
+            JuMP.set_lower_bound(zU_var[i,j,k,d], -1)
+            JuMP.set_upper_bound(zU_var[i,j,k,d],  1)
+        end
+    end 
     
     return
 end
@@ -96,16 +147,37 @@ end
 function variable_slack_for_feasibility(qcm::QuantumCircuitModel)
     n_r     = size(qcm.data["gates_real"])[1]
     n_c     = size(qcm.data["gates_real"])[2]
-
-    qcm.variables[:slack_var] = JuMP.@variable(qcm.model, -1 <= slack_var[1:n_r, 1:n_c] <= 1)
+    max_depth   = qcm.data["maximum_depth"]
+    U_var = qcm.variables[:U_var]
     
+    qcm.variables[:slack_var] = JuMP.@variable(qcm.model, slack_var[1:n_r, 1:n_c])
+
+    for i=1:n_r, j=1:n_c
+        lb = JuMP.lower_bound(U_var[i,j,max_depth-1])
+        ub = JuMP.upper_bound(U_var[i,j,max_depth-1])
+        if isapprox(lb, 0, atol = 1E-6) && isapprox(ub, 0, atol = 1E-6)
+            JuMP.set_lower_bound(slack_var[i,j], 0)
+            JuMP.set_upper_bound(slack_var[i,j], 0)
+        else 
+            JuMP.set_lower_bound(slack_var[i,j], -1)
+            JuMP.set_upper_bound(slack_var[i,j],  1)
+        end
+    end
+    return
+end
+
+function variable_slack_var_outer_approximation(qcm::QuantumCircuitModel)
+    n_r     = size(qcm.data["gates_real"])[1]
+    n_c     = size(qcm.data["gates_real"])[2]
+
+    qcm.variables[:slack_var_oa] = JuMP.@variable(qcm.model, -1 <= slack_var_oa[1:n_r, 1:n_c] <= 1)
     return
 end
 
 function variable_binary_products(qcm::QuantumCircuitModel)
     num_gates = size(qcm.data["gates_real"])[3]
-    depth     = qcm.data["maximum_depth"]
+    max_depth = qcm.data["maximum_depth"]
     
-    qcm.variables[:Z_var] = JuMP.@variable(qcm.model, 0 <= Z[1:num_gates, 1:num_gates, 1:depth] <= 1)
+    qcm.variables[:Z_var] = JuMP.@variable(qcm.model, 0 <= Z[1:num_gates, 1:num_gates, 1:max_depth] <= 1)
     return
 end
\ No newline at end of file
diff --git a/test/data_tests.jl b/test/data_tests.jl
index 604f92d..031fbaa 100644
--- a/test/data_tests.jl
+++ b/test/data_tests.jl
@@ -1,6 +1,6 @@
 # Unit tests for functions in data.jl
 
-@testset "Tests: building elementary universal gate" begin
+@testset "Data tests: building elementary universal gate" begin
     test_angle = π/3
     pauli_Y = QCO.YGate()
     H = QCO.HGate()
@@ -28,7 +28,7 @@
     @test isapprox(QCO.RZGate(test_angle), test_U1_1)
 end 
 
-@testset "Tests: get_full_sized_gate" begin
+@testset "Data tests: get_full_sized_gate" begin
 
     # 2-qubit gates
     params = Dict{String, Any}(
@@ -90,7 +90,7 @@ end
     @test length(keys(data["gates_dict"])) == 17
 end
 
-@testset "Tests: get_input_circuit_dict" begin
+@testset "Data tests: get_input_circuit_dict" begin
     
     function input_circuit_1()
         # [(depth, gate)]
diff --git a/test/gates_tests.jl b/test/gates_tests.jl
index 9b3d0b4..4d3dcd9 100644
--- a/test/gates_tests.jl
+++ b/test/gates_tests.jl
@@ -1,6 +1,6 @@
 # Unit tests for functions in gates.jl
 
-@testset "Tests: elementary universal gates in gates.jl" begin
+@testset "Gates tests: elementary universal gates in gates.jl" begin
     @test isapprox(QCO.U2Gate(0,-π/4), QCO.U3Gate(π/2,0,-π/4), atol=tol_0)
     @test isapprox(QCO.U1Gate(-π/4), QCO.U3Gate(0,0,-π/4), atol=tol_0)
 
diff --git a/test/qc_model_tests.jl b/test/qc_model_tests.jl
index 971367b..1fd814a 100644
--- a/test/qc_model_tests.jl
+++ b/test/qc_model_tests.jl
@@ -1,4 +1,4 @@
-@testset "Tests: Minimize depth for controlled-NOT gate" begin
+@testset "QC_model Tests: Minimize depth for controlled-NOT gate" begin
 
     params = Dict{String, Any}(
     "num_qubits" => 2, 
@@ -15,7 +15,7 @@
     model_options = Dict{Symbol, Any}(:model_type => "compact_formulation_1", # Testing incorrect model_type
                                       :all_valid_constraints => 2)            # Testing incorrect all_valid_constraints
 
-    result_qc = QCO.run_QCModel(params, CBC; options = model_options)
+    result_qc = QCO.run_QCModel(params, MIP_SOLVER; options = model_options)
 
     @test result_qc["termination_status"] == MOI.OPTIMAL
     @test result_qc["primal_status"]      == MOI.FEASIBLE_POINT
@@ -28,7 +28,7 @@
     
 end
 
-@testset "Tests: Minimum CNOT swap gate decomposition" begin
+@testset "QC_model Tests: Minimum CNOT swap gate decomposition" begin
 
     params = Dict{String, Any}(
         "num_qubits" => 2,
@@ -42,10 +42,9 @@ end
         )
     
     model_options = Dict{Symbol, Any}(:model_type => "balas_formulation",
-                                      :all_valid_constraints => 1,
-                                      :unit_magnitude_constraints => true)
+                                      :all_valid_constraints => 1)
 
-    result_qc = QCO.run_QCModel(params, CBC; options = model_options)
+    result_qc = QCO.run_QCModel(params, MIP_SOLVER; options = model_options)
 
     @test result_qc["termination_status"] == MOI.OPTIMAL
     @test result_qc["primal_status"]      == MOI.FEASIBLE_POINT
@@ -54,7 +53,7 @@ end
     
 end
 
-@testset "Tests: Minimum depth U3(0,0,π/4) gate" begin
+@testset "QC_model Tests: Minimum depth U3(0,0,π/4) gate" begin
 
     params = Dict{String, Any}(
     "num_qubits" => 2, 
@@ -68,7 +67,7 @@ end
     "decomposition_type" => "exact_optimal"                  
     )
 
-    result_qc = QCO.run_QCModel(params, CBC)
+    result_qc = QCO.run_QCModel(params, MIP_SOLVER)
 
     data = QCO.get_data(params)
 
@@ -81,7 +80,7 @@ end
     @test data["gates_dict"]["4"]["qubit_loc"] == "qubit_1"
 end
 
-@testset "Tests: Minimum depth CU3(0,0,π/4) gate" begin
+@testset "QC_model Tests: Minimum depth CU3(0,0,π/4) gate" begin
 
     params = Dict{String, Any}(
     "num_qubits" => 2, 
@@ -95,7 +94,7 @@ end
     "decomposition_type" => "exact_optimal"                  
     )
 
-    result_qc = QCO.run_QCModel(params, CBC)
+    result_qc = QCO.run_QCModel(params, MIP_SOLVER)
 
     data = QCO.get_data(params)
 
@@ -108,7 +107,7 @@ end
     @test data["gates_dict"]["4"]["qubit_loc"] == "qubit_1_2"
 end
 
-@testset "Tests: Minimum depth Rev CU3(0,0,π/4) gate" begin
+@testset "QC_model Tests: Minimum depth Rev CU3(0,0,π/4) gate" begin
 
     params = Dict{String, Any}(
     "num_qubits" => 3, 
@@ -122,7 +121,7 @@ end
     "decomposition_type" => "exact_optimal"                  
     )
 
-    result_qc = QCO.run_QCModel(params, CBC)
+    result_qc = QCO.run_QCModel(params, MIP_SOLVER)
 
     data = QCO.get_data(params)
 
@@ -135,7 +134,7 @@ end
     @test data["gates_dict"]["4"]["qubit_loc"] == "qubit_3_1"
 end
 
-@testset "Tests: Minimum depth RX, RY, RZ gate decomposition" begin
+@testset "QC_model Tests: Minimum depth RX, RY, RZ gate decomposition" begin
 
     params = Dict{String, Any}(
     "num_qubits" => 2, 
@@ -152,7 +151,7 @@ end
     model_options = Dict{Symbol, Any}(:model_type => "balas_formulation",
                                       :commute_gate_constraints => true)
 
-    result_qc = QCO.run_QCModel(params, CBC; options =  model_options)
+    result_qc = QCO.run_QCModel(params, MIP_SOLVER; options =  model_options)
 
     @test result_qc["termination_status"] == MOI.OPTIMAL
     @test result_qc["primal_status"]      == MOI.FEASIBLE_POINT
@@ -160,7 +159,7 @@ end
 
 end
 
-@testset "Tests: Minimum depth CRX, CRY, CRZ gate decomposition" begin
+@testset "QC_model Tests: Minimum depth CRX, CRY, CRZ gate decomposition" begin
 
     params = Dict{String, Any}(
     "num_qubits" => 3, 
@@ -177,7 +176,7 @@ end
     model_options = Dict{Symbol, Any}(:model_type => "balas_formulation",
                                       :commute_gate_constraints => true)
 
-    result_qc = QCO.run_QCModel(params, CBC; options =  model_options)
+    result_qc = QCO.run_QCModel(params, MIP_SOLVER; options =  model_options)
 
     @test result_qc["termination_status"] == MOI.OPTIMAL
     @test result_qc["primal_status"]      == MOI.FEASIBLE_POINT
@@ -185,7 +184,7 @@ end
 
 end
 
-@testset "Tests: Minimum depth Rev CRX, CRY, CRZ gate decomposition" begin
+@testset "QC_model Tests: Minimum depth Rev CRX, CRY, CRZ gate decomposition" begin
 
     params = Dict{String, Any}(
     "num_qubits" => 3, 
@@ -202,7 +201,7 @@ end
     model_options = Dict{Symbol, Any}(:model_type => "balas_formulation",
                                       :commute_gate_constraints => true)
 
-    result_qc = QCO.run_QCModel(params, CBC; options =  model_options)
+    result_qc = QCO.run_QCModel(params, MIP_SOLVER; options =  model_options)
 
     @test result_qc["termination_status"] == MOI.OPTIMAL
     @test result_qc["primal_status"]      == MOI.FEASIBLE_POINT
@@ -211,7 +210,7 @@ end
 end
 
 
-@testset "Tests: 3-qubit RX gate decomposition" begin
+@testset "QC_model Tests: 3-qubit RX gate decomposition" begin
 
     params = Dict{String, Any}(
     "num_qubits" => 3, 
@@ -227,7 +226,7 @@ end
 
     model_options = Dict{Symbol, Any}(:model_type => "balas_formulation")
 
-    result_qc = QCO.run_QCModel(params, CBC; options =  model_options)
+    result_qc = QCO.run_QCModel(params, MIP_SOLVER; options =  model_options)
     
     data = QCO.get_data(params)
 
@@ -243,7 +242,7 @@ end
 
 end
 
-@testset "Tests: 3-qubit CRX gate decomposition" begin
+@testset "QC_model Tests: 3-qubit CRX gate decomposition" begin
 
     params = Dict{String, Any}(
         "num_qubits" => 3, 
@@ -259,7 +258,7 @@ end
 
     model_options = Dict{Symbol, Any}(:model_type => "balas_formulation")
 
-    result_qc = QCO.run_QCModel(params, CBC; options =  model_options)
+    result_qc = QCO.run_QCModel(params, MIP_SOLVER; options =  model_options)
     
     data = QCO.get_data(params)
 
@@ -275,7 +274,7 @@ end
 
 end
 
-@testset "Tests: feasibility problem" begin
+@testset "QC_model Tests: feasibility problem" begin
     params = Dict{String, Any}(
         "num_qubits" => 2, 
         "maximum_depth" => 3,    
@@ -287,7 +286,7 @@ end
     
     model_options = Dict{Symbol, Any}(:all_valid_constraints => -1)
 
-    result_qc = QCO.run_QCModel(params, CBC; options =  model_options)
+    result_qc = QCO.run_QCModel(params, MIP_SOLVER; options =  model_options)
 
     @test result_qc["termination_status"] == MOI.OPTIMAL
     @test result_qc["primal_status"]      == MOI.FEASIBLE_POINT
@@ -301,7 +300,7 @@ end
   
 end
 
-@testset "Tests: JuMP set_start_value for z_bin_var variables" begin
+@testset "QC_model Tests: JuMP set_start_value for z_bin_var variables" begin
     function input_circuit()
         # [(depth, gate)]
         return [(1, "CNot_2_1"), 
@@ -321,7 +320,7 @@ end
     "decomposition_type" => "exact_optimal"
     )
 
-    result_qc = QCO.run_QCModel(params, CBC)
+    result_qc = QCO.run_QCModel(params, MIP_SOLVER)
     @test result_qc["termination_status"] == MOI.OPTIMAL
     @test result_qc["primal_status"]      == MOI.FEASIBLE_POINT
     @test isapprox(sum(result_qc["solution"]["z_bin_var"][6,:]), 1, atol=tol_0)
@@ -329,7 +328,7 @@ end
     
 end
 
-@testset "Tests: Involutory gate constraints" begin
+@testset "QC_model Tests: Involutory gate constraints" begin
     
     params = Dict{String, Any}(
     "num_qubits" => 2,
@@ -353,14 +352,14 @@ end
         end
     end
 
-    result_qc = QCO.run_QCModel(params, CBC)
+    result_qc = QCO.run_QCModel(params, MIP_SOLVER)
     @test result_qc["termination_status"] == MOI.OPTIMAL
     @test result_qc["primal_status"]      == MOI.FEASIBLE_POINT
     @test isapprox(result_qc["objective"], 4, atol = tol_0)
     
 end
 
-@testset "Tests: TIME_LIMIT for building results dict and log" begin
+@testset "QC_model Tests: TIME_LIMIT for building results dict and log" begin
     
     params = Dict{String, Any}(
     "num_qubits" => 2,
@@ -373,13 +372,13 @@ end
 
     model_options = Dict{Symbol, Any}(:time_limit => 1)
 
-    result_qc = QCO.run_QCModel(params, CBC; options = model_options)
+    result_qc = QCO.run_QCModel(params, MIP_SOLVER; options = model_options)
     @test result_qc["termination_status"] == MOI.TIME_LIMIT
     @test result_qc["primal_status"]      == MOI.NO_SOLUTION
     
 end
 
-@testset "Tests: constraint_redundant_gate_product_pairs" begin
+@testset "QC_model Tests: constraint_redundant_gate_product_pairs" begin
     
     function target_gate()
         T1 = QCO.get_full_sized_gate("U3_2", 2, angle = [0,π/2,π])
@@ -405,14 +404,14 @@ end
     @test isapprox(data["gates_dict"]["1"]["matrix"] * data["gates_dict"]["4"]["matrix"], data["gates_dict"]["2"]["matrix"], atol = tol_0)
     @test isapprox(data["gates_dict"]["5"]["matrix"] * data["gates_dict"]["7"]["matrix"], data["gates_dict"]["6"]["matrix"], atol = tol_0)
 
-    result_qc = QCO.run_QCModel(params, CBC)
+    result_qc = QCO.run_QCModel(params, MIP_SOLVER)
     @test result_qc["termination_status"] == MOI.OPTIMAL
     @test result_qc["primal_status"]      == MOI.FEASIBLE_POINT
     @test isapprox(result_qc["objective"], 2.0, atol = tol_0)
 
 end
 
-@testset "Tests: constraint_idempotent_gates" begin
+@testset "QC_model Tests: constraint_idempotent_gates" begin
     params = Dict{String, Any}(
     "num_qubits" => 2, 
     "maximum_depth" => 3,    
@@ -428,7 +427,7 @@ end
 
     model_options = Dict{Symbol, Any}(:idempotent_gate_constraints => true)
     
-    result_qc = QCO.run_QCModel(params, CBC; options = model_options)
+    result_qc = QCO.run_QCModel(params, MIP_SOLVER; options = model_options)
 
     @test result_qc["termination_status"] == MOI.OPTIMAL
     @test result_qc["primal_status"]      == MOI.FEASIBLE_POINT
@@ -442,7 +441,7 @@ end
 
 end
 
-@testset "Tests: constraint_convex_hull_complex_gates" begin
+@testset "QC_model Tests: constraint_convex_hull_complex_gates" begin
     
     function target_gate()
         num_qubits = 4
@@ -458,23 +457,25 @@ end
 
     params = Dict{String, Any}(
     "num_qubits" => 4,
-    "maximum_depth" => 2,
+    "maximum_depth" => 6,
     "elementary_gates" => ["CV_1_4", "CV_2_4", "CV_3_4", "CVdagger_3_4", "CNot_1_2", "CNot_2_3", "Identity"],
     "target_gate" => target_gate()
     )
     
     model_options = Dict{Symbol, Any}(:relax_integrality => true,
                                       :convex_hull_gate_constraints => true,
-                                      :optimizer_log => false)
+                                      :fix_unitary_variables => true,
+                                      :optimizer_log => false,
+                                      )
 
-    result_qc = QCO.run_QCModel(params, CBC; options = model_options)
+    result_qc = QCO.run_QCModel(params, MIP_SOLVER; options = model_options)
 
     @test result_qc["termination_status"] == MOI.OPTIMAL
     @test result_qc["primal_status"]      == MOI.FEASIBLE_POINT
-    @test isapprox(result_qc["objective"], 0.8275862068965518, atol = tol_0)
+    @test isapprox(result_qc["objective"], 0.6, atol = tol_0)
 end
 
-@testset "Tests: RGate decomposition" begin
+@testset "QC_model Tests: RGate decomposition" begin
     
     num_qubits = 2
     GR1      = QCO.GRGate(num_qubits, π/6, π/3)
@@ -492,8 +493,8 @@ end
     "decomposition_type" => "exact_optimal",
     )
     
-    model_options = Dict{Symbol, Any}(:optimizer_log => false)
-    result_qc = QCO.run_QCModel(params, CBC; options = model_options)
+    model_options = Dict{Symbol, Any}(:optimizer_log => false, :fix_unitary_variables => false)
+    result_qc = QCO.run_QCModel(params, MIP_SOLVER; options = model_options)
 
     @test result_qc["termination_status"] == MOI.OPTIMAL
     @test result_qc["primal_status"]      == MOI.FEASIBLE_POINT
@@ -503,7 +504,7 @@ end
     @test isapprox(sum(z_sol[9, :]),   1, atol = tol_0)
 end
 
-@testset "Tests: GRGate decomposition for Pauli-X gate, and unitary constraints" begin
+@testset "QC_model Tests: GRGate decomposition for Pauli-X gate, and unitary constraints" begin
 
     params = Dict{String, Any}(
     "num_qubits" => 2,
@@ -518,7 +519,7 @@ end
     )
     
     model_options = Dict{Symbol, Any}(:optimizer_log => false, :unitary_constraints => true)
-    result_qc = QCO.run_QCModel(params, CBC; options = model_options)
+    result_qc = QCO.run_QCModel(params, MIP_SOLVER; options = model_options)
 
     @test result_qc["termination_status"] == MOI.OPTIMAL
     @test result_qc["primal_status"]      == MOI.FEASIBLE_POINT
@@ -527,7 +528,7 @@ end
     @test isapprox(sum(z_sol[4, :]), 2, atol = tol_0) # test for number of CZ gates 
 end
 
-@testset "Tests: Decomposition using exact_feasible option" begin
+@testset "QC_model Tests: Decomposition using exact_feasible option" begin
     params = Dict{String, Any}(
         "num_qubits" => 2,
         "maximum_depth" => 5,
@@ -538,10 +539,82 @@ end
         )
 
     model_options = Dict{Symbol, Any}(:optimizer_log => false)
-    result_qc = QCO.run_QCModel(params, CBC; options = model_options)
+    result_qc = QCO.run_QCModel(params, MIP_SOLVER; options = model_options)
     @test result_qc["termination_status"] == MOI.OPTIMAL
     @test result_qc["primal_status"]      == MOI.FEASIBLE_POINT
     @test isapprox(result_qc["objective"], 0.0, atol = tol_0)
     z_sol = result_qc["solution"]["z_bin_var"]
     @test isapprox(sum(z_sol[2:4, :]), 4, atol = tol_0)
+end
+
+@testset "QC_model Tests: Approximate decomposition using outer approximation" begin
+    params = Dict{String, Any}(
+        "num_qubits" => 2,
+        "maximum_depth" => 5,
+        "elementary_gates" => ["H_1", "H_2", "CNot_1_2", "Identity"],
+        "target_gate" => QCO.CNotRevGate(),
+        "objective" => "minimize_depth",
+        "decomposition_type" => "approximate",
+        )
+    model_options = Dict{Symbol, Any}(:optimizer_log => false)
+    result_qc = QCO.run_QCModel(params, MIP_SOLVER; options = model_options)
+    @test result_qc["termination_status"] == MOI.OPTIMAL
+    @test result_qc["primal_status"]      == MOI.FEASIBLE_POINT
+    @test isapprox(result_qc["objective"], 5.0, atol = tol_0)
+end
+
+@testset "QC_model Tests: Approximate decomposition using outer approximation" begin
+    params = Dict{String, Any}(
+        "num_qubits" => 2,
+        "maximum_depth" => 5,
+        "elementary_gates" => ["H_1", "H_2", "CNot_1_2", "Identity"],
+        "target_gate" => QCO.CNotRevGate(),
+        "objective" => "minimize_depth",
+        "decomposition_type" => "approximate",
+        )
+    model_options = Dict{Symbol, Any}(:optimizer_log => false)
+    result_qc = QCO.run_QCModel(params, MIP_SOLVER; options = model_options)
+    @test result_qc["termination_status"] == MOI.OPTIMAL
+    @test result_qc["primal_status"]      == MOI.FEASIBLE_POINT
+    @test isapprox(result_qc["objective"], 5.0, atol = tol_0)
+
+    # Testing approximate decomposition for balas_formulation
+    model_options = Dict{Symbol, Any}(:optimizer_log => false, :model_type => "balas_formulation")
+    result_qc = QCO.run_QCModel(params, MIP_SOLVER; options = model_options)
+    @test result_qc["termination_status"] == MOI.OPTIMAL
+    @test result_qc["primal_status"]      == MOI.FEASIBLE_POINT
+    @test isapprox(result_qc["objective"], 5.0, atol = tol_0)
+
+    # Testing approximate decomposition for feasibility case
+    params["elementary_gates"] = ["H_1", "H_2", "CNot_1_2"]
+    model_options = Dict{Symbol, Any}(:optimizer_log => false)
+    result_qc = QCO.run_QCModel(params, MIP_SOLVER; options = model_options)
+    @test result_qc["termination_status"] == MOI.OPTIMAL
+    @test result_qc["primal_status"]      == MOI.FEASIBLE_POINT
+    @test isapprox(result_qc["objective"], 0.0, atol = tol_0)
+
+    # Testing approximate decomposition for minimizing CNOT gates
+    params["elementary_gates"] = ["H_1", "H_2", "CNot_1_2", "Identity"]
+    params["objective"] = "minimize_cnot"
+    model_options = Dict{Symbol, Any}(:optimizer_log => false)
+    result_qc = QCO.run_QCModel(params, MIP_SOLVER; options = model_options)
+    @test result_qc["termination_status"] == MOI.OPTIMAL
+    @test result_qc["primal_status"]      == MOI.FEASIBLE_POINT
+    @test isapprox(result_qc["objective"], 1.0, atol = tol_0)
+
+    # Testing approximate decomposition for case when slack_var-s are fixed based on U_var-s
+    params = Dict{String, Any}(
+    "num_qubits" => 3,
+    "maximum_depth" => 7,
+    "elementary_gates" => ["CV_1_2", "CV_2_3", "CV_1_3", "CVdagger_1_2", "CVdagger_2_3", "CVdagger_1_3", "CNot_1_2", "CNot_3_2", "CNot_2_3", "CNot_1_3", "Identity"],
+    "target_gate" => QCO.CSwapGate(), #also Fredkin
+    "objective" => "minimize_depth",
+    "decomposition_type" => "approximate"
+    )
+    model_options = Dict{Symbol, Any}(:optimizer_log => false, :relax_integrality => true, :fix_unitary_variables => true)
+    result_qc = QCO.run_QCModel(params, MIP_SOLVER; options = model_options)
+    @test result_qc["termination_status"] == MOI.OPTIMAL
+    @test result_qc["primal_status"]      == MOI.FEASIBLE_POINT
+    @test isapprox(result_qc["objective"], 0.38281250, atol = tol_0)
+
 end
\ No newline at end of file
diff --git a/test/relaxations_tests.jl b/test/relaxations_tests.jl
index 5ed3c0b..aab5e1b 100644
--- a/test/relaxations_tests.jl
+++ b/test/relaxations_tests.jl
@@ -1,7 +1,7 @@
-@testset "Tests: relaxation_bilinear" begin
+@testset "Relaxation tests: relaxation_bilinear" begin
     QCO.silence()
     
-    m = JuMP.Model(CBC)
+    m = JuMP.Model(MIP_SOLVER)
 
     LB = [-1, -2.5, 0, -3, 0]
     UB = [2.5, 0, 1, 3.2, 1]
diff --git a/test/runtests.jl b/test/runtests.jl
index fd09a73..30aa82b 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,17 +1,22 @@
-using QuantumCircuitOpt
-using Memento
 using JuMP
-using LinearAlgebra
-using Test
-using Cbc
+using Test 
 
-const QCO = QuantumCircuitOpt
-const LA = LinearAlgebra
+import QuantumCircuitOpt as QCO
+import Memento
+import LinearAlgebra as LA
+import HiGHS
 
 # Suppress warnings during testing
 QCO.logger_config!("error")
 
-const CBC = JuMP.optimizer_with_attributes(Cbc.Optimizer,  MOI.Silent() => true)
+const HIGHS = MOI.OptimizerWithAttributes(
+    HiGHS.Optimizer,
+    "presolve" => "on",
+    "log_to_console" => false,
+)
+
+const MIP_SOLVER = HIGHS
+
 tol_0 = 1E-6
 
 @testset "QuantumCircuitOpt" begin
diff --git a/test/utility_tests.jl b/test/utility_tests.jl
index 0609f1d..04b574d 100644
--- a/test/utility_tests.jl
+++ b/test/utility_tests.jl
@@ -172,4 +172,57 @@ end
     @test c == 2.5
 
     QCO._get_constraint_slope_intercept(v1, v1)
+end
+
+@testset "Tests: U_var fixed indices" begin
+    num_qubits = 2
+    maximum_depth = 3
+
+    # Input gates
+    H1 = QCO.complex_to_real_gate(QCO.get_full_sized_gate("H_1", num_qubits))
+    CNot_1_2 = QCO.complex_to_real_gate(QCO.get_full_sized_gate("CNot_1_2", num_qubits))
+    Id = QCO.complex_to_real_gate(QCO.IGate(num_qubits))
+
+    # Assuming only H1 in elementary gates
+    n_r = size(H1)[1]
+    gates = zeros(n_r,n_r,1)
+    gates[:,:,1] = H1
+
+    U_fixed_idx = QCO._get_unitary_variables_fixed_indices(gates, maximum_depth)
+    
+    # Depth = 1
+    for idx in keys(U_fixed_idx[1])
+        @test isapprox(U_fixed_idx[1][idx]["value"], gates[idx[1],idx[2],1], atol = tol_0)
+    end
+
+    # Depth = 2 (product of H1 is an identity matrix)
+    @test length(keys(U_fixed_idx[2])) == n_r * n_r
+    for idx in keys(U_fixed_idx[2])
+        @test isapprox(U_fixed_idx[2][idx]["value"], Id[idx[1], idx[2]], atol = tol_0)
+    end
+
+    # Assuming H1 and CNot_1_2 in elementary gates
+    gates = zeros(n_r,n_r,2)
+    gates[:,:,1] = H1
+    gates[:,:,2] = CNot_1_2
+    U_fixed_idx = QCO._get_unitary_variables_fixed_indices(gates, maximum_depth)
+    
+    # Depth = 1        
+    sum_mat = abs.(gates[:,:,1]) + abs.(gates[:,:,2])
+    for idx in keys(U_fixed_idx[1])
+        if isapprox(U_fixed_idx[1][idx]["value"], 0, atol = tol_0)
+            @test isapprox(U_fixed_idx[1][idx]["value"], sum_mat[idx[1], idx[2]], atol = tol_0)
+        end
+    end
+    
+    # Depth = 2
+    prod_mat_1 = gates[:,:,1] * gates[:,:,2]
+    prod_mat_2 = gates[:,:,2] * gates[:,:,1]
+    for idx in keys(U_fixed_idx[2])
+        if isapprox(U_fixed_idx[2][idx]["value"], 0, atol = tol_0)
+            @test isapprox(U_fixed_idx[2][idx]["value"], prod_mat_1[idx[1], idx[2]], atol = tol_0)
+            @test isapprox(U_fixed_idx[2][idx]["value"], prod_mat_2[idx[1], idx[2]], atol = tol_0)
+        end
+    end
+
 end
\ No newline at end of file