diff --git a/cpmpy/expressions/globalfunctions.py b/cpmpy/expressions/globalfunctions.py
index 515947a6c..0d7882197 100644
--- a/cpmpy/expressions/globalfunctions.py
+++ b/cpmpy/expressions/globalfunctions.py
@@ -10,7 +10,7 @@
     ------------------------
 
     If a solver does not support such a global function (see solvers/), then it will be automatically
-    decomposed by calling its `.decompose_comparison()` function.
+    decomposed by calling its `.decompose()` function.
 
     CPMpy GlobalFunctions does not exactly match what is implemented in the solvers.
     Solvers can have specialised implementations for global functions, when used in a comparison, as global constraints.
@@ -42,8 +42,8 @@ class my_global(GlobalFunction):
             def __init__(self, args):
                 super().__init__("my_global", args)
 
-            def decompose_comparison(self):
-                return [self.args[0] + self.args[1]] # your decomposition
+            def decompose(self):
+                return [self.args[0] + self.args[1]], [] # your decomposition
 
     Also, implement `.value()` accordingly.
 
@@ -65,6 +65,7 @@ def decompose_comparison(self):
 
 """
 import warnings  # for deprecation warning
+
 import numpy as np
 import cpmpy as cp
 
@@ -87,9 +88,10 @@ def is_bool(self):
         """
         return False
 
-    def decompose_comparison(self, cmp_op, cmp_rhs):
+    def decompose(self):
         """
             Returns a decomposition into smaller constraints.
+            Returns a numerical expression and a list of defining constraints.
 
             The decomposition might create auxiliary variables
             and use other global constraints as long as
@@ -97,6 +99,20 @@ def decompose_comparison(self, cmp_op, cmp_rhs):
         """
         raise NotImplementedError("Decomposition for", self, "not available")
 
+    def decompose_comparison(self, cmp_op, cmp_rhs):
+        """
+            Returns a decomposition into smaller constraints.
+
+            The decomposition might create auxiliary variables
+            and use other global constraints as long as
+            it does not create a circular dependency.
+        """
+        warnings.warn(f"Deprecated, use {self}.decompose() instead, will be removed in "
+                      "stable version", DeprecationWarning)
+        val, tl = self.decompose()
+        return eval_comparison(cmp_op, val, cmp_rhs), tl
+
+
     def get_bounds(self):
         """
         Returns the bounds of the global function
@@ -126,20 +142,15 @@ def value(self):
         else:
             return min(argvals)
 
-    def decompose_comparison(self, cpm_op, cpm_rhs):
+    def decompose(self):
         """
-        Decomposition if it's part of a comparison
-
-        Returns two lists of constraints:
-
-        1) constraints representing the comparison
-        2) constraints that (totally) define new auxiliary variables needed in the decomposition,
-           they should be enforced toplevel.
+        Decomposition of Minimum constraint.
+        Returns
+        1) a numerical value to replace the constraint, and
+        2) a list of defining constraints, which should be enforced toplevel
         """
-        lb, ub = self.get_bounds()
-        _min = intvar(lb, ub)
-        return [eval_comparison(cpm_op, _min, cpm_rhs)], \
-               [cp.any(x <= _min for x in self.args), cp.all(x >= _min for x in self.args), ]
+        _min = intvar(*self.get_bounds())
+        return _min, [cp.all(x >= _min for x in self.args), cp.any(x <= _min for x in self.args)]
 
     def get_bounds(self):
         """
@@ -164,20 +175,15 @@ def value(self):
         else:
             return max(argvals)
 
-    def decompose_comparison(self, cpm_op, cpm_rhs):
+    def decompose(self):
         """
-        Decomposition if it's part of a comparison
-
-        Returns two lists of constraints:
-
-        1) constraints representing the comparison
-        2) constraints that (totally) define new auxiliary variables needed in the decomposition,
-           they should be enforced toplevel.
+        Decomposition of Maximum constraint.
+        Returns
+        1) a numerical value to replace the constraint, and
+        2) a list of defining constraints, which should be enforced toplevel
         """
-        lb, ub = self.get_bounds()
-        _max = intvar(lb, ub)
-        return [eval_comparison(cpm_op, _max, cpm_rhs)], \
-               [cp.any(x >= _max for x in self.args), cp.all(x <= _max for x in self.args)]
+        _max = intvar(*self.get_bounds())
+        return _max, [cp.all(x <= _max for x in self.args), cp.any(x >= _max for x in self.args)]
 
     def get_bounds(self):
         """
@@ -197,30 +203,25 @@ def __init__(self, expr):
     def value(self):
         return abs(argval(self.args[0]))
 
-    def decompose_comparison(self, cpm_op, cpm_rhs):
+    def decompose(self):
         """
-        Decomposition if it's part of a comparison
-
-        Returns two lists of constraints:
-
-        1) constraints representing the comparison
-        2) constraints that (totally) define new auxiliary variables needed in the decomposition,
-           they should be enforced toplevel.
+        Decomposition of Abs constraint.
+        Returns
+        1) a numerical value to replace the constraint, and
+        2) a list of defining constraints, which should be enforced toplevel
         """
         arg = self.args[0]
         lb, ub = get_bounds(arg)
-        # when argument is exclusively on one side of the sign
-        if lb >= 0:
-            return [eval_comparison(cpm_op, arg, cpm_rhs)], []
-        elif ub <= 0:
-            return [eval_comparison(cpm_op, -arg, cpm_rhs)], []
-        else: # when domain crosses over 0
-            newarg = intvar(*self.get_bounds())
-            is_pos = boolvar()
-            return [eval_comparison(cpm_op, newarg, cpm_rhs)], \
-                    [is_pos == (arg >= 0), is_pos.implies(arg == newarg), (~is_pos).implies(-arg == newarg)]
+        if lb >= 0: # always positive
+            return arg, []
+        if ub <= 0: # always negative
+            return -arg, []
 
+        _abs = intvar(*self.get_bounds())
+        assert _abs.lb == 0
 
+        is_pos = arg >= 0
+        return _abs, [is_pos.implies(arg == _abs), (~is_pos).implies(arg == -_abs)]
 
     def get_bounds(self):
         """
@@ -274,29 +275,22 @@ def value(self):
                                           + "\n Use argval(expr) to get the value of expr with relational semantics.")
         return None # default
 
-    def decompose_comparison(self, cpm_op, cpm_rhs):
+    def decompose(self):
         """
-            `Element(arr,ix)` represents the array lookup itself (a numeric variable)
-            When used in a comparison relation: Element(arr,idx) <CMP_OP> CMP_RHS
-            it is a constraint, and that one can be decomposed.
+        Decomposition of Abs constraint.
+        Returns
+        1) a numerical value to replace the constraint, and
+        2) a list of defining constraints, which should be enforced toplevel
+        """
+        arr, idx = self.args
 
-            Returns two lists of constraints:
+        idx_lb, idx_ub = get_bounds(idx)
+        assert idx_lb >= 0 and idx_ub < len(arr), "Element constraint is unsafe to decompose as it can be partial. Safen first using `cpmpy.transformations.safening.no_partial_functions`"
 
-            1) constraints representing the comparison
-            2) constraints that (totally) define new auxiliary variables needed in the decomposition,
-               they should be enforced toplevel.
+        _elem = intvar(*self.get_bounds())
+
+        return _elem, [implies(idx == i, _elem == arr[i]) for i in range(len(arr))]
 
-        """
-        arr, idx = self.args
-        # Find where the array indices and the bounds of `idx` intersect
-        lb, ub = get_bounds(idx)
-        new_lb, new_ub = max(lb, 0), min(ub, len(arr) - 1)
-        cons=[]
-        # For every `i` in that intersection, post `(idx = i) -> idx=i -> arr[i] <CMP_OP> cpm_rhs`.
-        for i in range(new_lb, new_ub+1):
-            cons.append(implies(idx == i, eval_comparison(cpm_op, arr[i], cpm_rhs)))
-        cons+=[idx >= new_lb, idx <= new_ub]  # also enforce the new bounds 
-        return cons, []  # no auxiliary variables
 
     def __repr__(self):
         return "{}[{}]".format(self.args[0], self.args[1])
@@ -320,12 +314,15 @@ def __init__(self,arr,val):
             raise TypeError("count takes an array and a value as input, not: {} and {}".format(arr,val))
         super().__init__("count", [arr,val])
 
-    def decompose_comparison(self, cmp_op, cmp_rhs):
+    def decompose(self):
         """
-        Count(arr,val) can only be decomposed if it's part of a comparison
+        Decomposition of the Count constraint.
+        Returns
+        1) a numerical value to replace the constraint, and
+        2) a list of defining constraints, which should be enforced toplevel
         """
         arr, val = self.args
-        return [eval_comparison(cmp_op, Operator('sum',[ai==val for ai in arr]), cmp_rhs)], []
+        return cp.sum(a == val for a in arr), []
 
     def value(self):
         arr, val = self.args
@@ -353,13 +350,17 @@ def __init__(self,arr,vals):
             raise TypeError(f"Among takes a set of values as input, not {vals}")
         super().__init__("among", [arr,vals])
 
-    def decompose_comparison(self, cmp_op, cmp_rhs):
+    def decompose(self):
         """
-            Among(arr, vals) can only be decomposed if it's part of a comparison'
+             Decomposition of the Among constraint.
+        Returns
+        1) a numerical value to replace the constraint, and
+        2) a list of defining constraints, which should be enforced toplevel
+
         """
         arr, values = self.args
-        count_for_each_val = [Count(arr, val) for val in values]
-        return [eval_comparison(cmp_op, cp.sum(count_for_each_val), cmp_rhs)], []
+        return cp.sum(Count(arr, val) for val in values), []
+
 
     def value(self):
         return int(sum(np.isin(argvals(self.args[0]), self.args[1])))
@@ -379,31 +380,27 @@ def __init__(self, arr):
             raise ValueError("NValue takes an array as input")
         super().__init__("nvalue", arr)
 
-    def decompose_comparison(self, cmp_op, cpm_rhs):
+    def decompose(self):
         """
-        NValue(arr) can only be decomposed if it's part of a comparison
+        Decomposition of the Count constraint.
+        Returns
+        1) a numerical value to replace the constraint, and
+        2) a list of defining constraints, which should be enforced toplevel
 
         Based on "simple decomposition" from:
-        
+
             Bessiere, Christian, et al. "Decomposition of the NValue constraint."
             International Conference on Principles and Practice of Constraint Programming.
             Berlin, Heidelberg: Springer Berlin Heidelberg, 2010.
         """
-
         lbs, ubs = get_bounds(self.args)
         lb, ub = min(lbs), max(ubs)
 
-        constraints = []
-
-        # introduce boolvar for each possible value
-        bvars = boolvar(shape=(ub+1-lb,)) # shape is tuple to ensure it is a 1D array
+        n_values = 0
+        for v in range(lb, ub+1):
+            n_values += cp.any(a == v for a in self.args)
 
-        args = cpm_array(self.args)
-        # bvar is true if the value is taken by any variable
-        for bv, val in zip(bvars, range(lb, ub+1)):
-            constraints += [cp.any(args == val) == bv]
-
-        return [eval_comparison(cmp_op, cp.sum(bvars), cpm_rhs)], constraints
+        return n_values, []
 
     def value(self):
         return len(set(argval(a) for a in self.args))
@@ -429,9 +426,12 @@ def __init__(self, arr, n):
             raise ValueError(f"NValueExcept takes an integer as second argument, but got {n} of type {type(n)}")
         super().__init__("nvalue_except",[arr, n])
 
-    def decompose_comparison(self, cmp_op, cpm_rhs):
+    def decompose(self):
         """
-        NValue(arr) can only be decomposed if it's part of a comparison
+        Decomposition of the Count constraint.
+        Returns
+        1) a numerical value to replace the constraint, and
+        2) a list of defining constraints, which should be enforced toplevel
 
         Based on "simple decomposition" from:
 
@@ -439,27 +439,19 @@ def decompose_comparison(self, cmp_op, cpm_rhs):
             International Conference on Principles and Practice of Constraint Programming.
             Berlin, Heidelberg: Springer Berlin Heidelberg, 2010.
         """
-
         arr, n = self.args
-        arr = cpm_array(arr)
+        assert is_num(n)
+
         lbs, ubs = get_bounds(arr)
         lb, ub = min(lbs), max(ubs)
 
-        constraints = []
-
-        # introduce boolvar for each possible value
-        bvars = boolvar(shape=(ub+1-lb,)) # shape is tuple to ensure it is a 1D array
-        idx_of_n = n - lb
-        if 0 <= idx_of_n < len(bvars):
-            count_of_vals = cp.sum(bvars[:idx_of_n]) + cp.sum(bvars[idx_of_n+1:])
-        else:
-            count_of_vals = cp.sum(bvars)
-
-        # bvar is true if the value is taken by any variable
-        for bv, val in zip(bvars, range(lb, ub + 1)):
-            constraints += [cp.any(arr == val) == bv]
+        n_values = 0
+        for v in range(lb, ub+1):
+            if v == n:
+                continue
+            n_values += cp.any(a == v for a in arr)
 
-        return [eval_comparison(cmp_op, count_of_vals, cpm_rhs)], constraints
+        return n_values, []
 
     def value(self):
         return len(set(argval(a) for a in self.args[0]) - {self.args[1]})
diff --git a/cpmpy/solvers/TEMPLATE.py b/cpmpy/solvers/TEMPLATE.py
index dd3626779..ad2488e54 100644
--- a/cpmpy/solvers/TEMPLATE.py
+++ b/cpmpy/solvers/TEMPLATE.py
@@ -59,10 +59,12 @@
 from ..expressions.utils import is_num, is_any_list, is_boolexpr
 from ..transformations.get_variables import get_variables
 from ..transformations.normalize import toplevel_list
-from ..transformations.decompose_global import decompose_in_tree
-from ..transformations.flatten_model import flatten_constraint
+from ..transformations.decompose_global import decompose_in_tree, decompose_objective
+from ..transformations.flatten_model import flatten_constraint, flatten_objective
 from ..transformations.comparison import only_numexpr_equality
 from ..transformations.reification import reify_rewrite, only_bv_reifies
+from ..transformations.safening import safen_objective
+
 
 class CPM_template(SolverInterface):
     """
@@ -292,21 +294,28 @@ def objective(self, expr, minimize=True):
 
             are permanently posted to the solver)
         """
-        # make objective function non-nested
-        (flat_obj, flat_cons) = flatten_objective(expr, csemap=self._csemap)
-        self += flat_cons # add potentially created constraints
-        self.user_vars.update(get_variables(flat_obj)) # add objvars to vars
+
+        # save user variables
+        get_variables(expr, self.user_vars)
+
+        # transform objective
+        obj, safe_cons = safen_objective(expr)
+        # [GUIDELINE] all globals which unsupported in a nested context should be decomposed here.
+        obj, decomp_cons = decompose_objective(obj, supported={"min", "max", "element"}, csemap=self._csemap)
+        obj, flat_cons = flatten_objective(obj, csemap=self._csemap)
+
+        self.add(safe_cons + decomp_cons + flat_cons)
 
         # make objective function or variable and post
-        obj = self._make_numexpr(flat_obj)
+        tpl_obj = self._make_numexpr(obj)
         # [GUIDELINE] if the solver interface does not provide a solver native "numeric expression" object,
         #         _make_numexpr may be removed and an objective can be posted as:
         #           self.TPL_solver.MinimizeWeightedSum(obj.args[0], self.solver_vars(obj.args[1]) or similar
 
         if minimize:
-            self.TPL_solver.Minimize(obj)
+            self.TPL_solver.Minimize(tpl_obj)
         else:
-            self.TPL_solver.Maximize(obj)
+            self.TPL_solver.Maximize(tpl_obj)
 
     def has_objective(self):
         return self.TPL_solver.hasObjective()
diff --git a/cpmpy/solvers/choco.py b/cpmpy/solvers/choco.py
index de432e0da..773e6d561 100644
--- a/cpmpy/solvers/choco.py
+++ b/cpmpy/solvers/choco.py
@@ -57,9 +57,9 @@
 from ..expressions.variables import _NumVarImpl, _IntVarImpl, _BoolVarImpl, NegBoolView, intvar
 from ..expressions.globalconstraints import GlobalConstraint
 from ..expressions.utils import is_num, is_int, is_boolexpr, is_any_list, get_bounds, argval, argvals, STAR
-from ..transformations.decompose_global import decompose_in_tree
+from ..transformations.decompose_global import decompose_in_tree, decompose_objective
 from ..transformations.get_variables import get_variables
-from ..transformations.flatten_model import flatten_constraint
+from ..transformations.flatten_model import flatten_constraint, get_or_make_var
 from ..transformations.comparison import only_numexpr_equality
 from ..transformations.linearize import canonical_comparison
 from ..transformations.safening import no_partial_functions
@@ -351,13 +351,23 @@ def objective(self, expr, minimize):
 
             .. note::
                 technical side note: constraints created during conversion of the objective
-                are premanently posted to the solver. Choco accepts variables to maximize or minimize
+                are permanently posted to the solver. Choco accepts variables to maximize or minimize
                 so it is needed to post constraints and create auxiliary variables
         """
 
+        # save user vars
+        get_variables(expr, self.user_vars)
+
+        # transform objective
+        supported = {"min", "max", "abs", "count", "element", "alldifferent", "alldifferent_except0", "allequal",
+                     "table", 'negative_table', "short_table", "regular", "InDomain", "cumulative", "circuit", "gcc", "inverse", "nvalue", "increasing",
+                     "decreasing","strictly_increasing","strictly_decreasing","lex_lesseq", "lex_less", "among", "precedence"}
+        obj, decomp_cons = decompose_objective(expr, supported=supported, csemap=self._csemap)
+
         # make objective function non-nested
-        obj_var = intvar(*get_bounds(expr))
-        self += obj_var == expr
+        obj_var, obj_cons = get_or_make_var(obj) # do not pass csemap here, we will still transform obj_var == obj...
+
+        self.add(decomp_cons + obj_cons)
 
         self.obj = obj_var
         self.minimize_obj = minimize  # Choco has as default to maximize
diff --git a/cpmpy/solvers/cplex.py b/cpmpy/solvers/cplex.py
index 79d256ca7..c294b9d5c 100644
--- a/cpmpy/solvers/cplex.py
+++ b/cpmpy/solvers/cplex.py
@@ -60,13 +60,15 @@
 from ..expressions.variables import _BoolVarImpl, NegBoolView, _IntVarImpl, _NumVarImpl, intvar
 from ..expressions.globalconstraints import DirectConstraint
 from ..transformations.comparison import only_numexpr_equality
-from ..transformations.decompose_global import decompose_in_tree
+from ..transformations.decompose_global import decompose_in_tree, decompose_objective
 from ..transformations.flatten_model import flatten_constraint, flatten_objective
 from ..transformations.get_variables import get_variables
-from ..transformations.linearize import linearize_constraint, only_positive_bv, only_positive_bv_wsum
+from ..transformations.linearize import linearize_constraint, only_positive_bv, only_positive_bv_wsum, \
+    only_positive_bv_wsum_const
 from ..transformations.normalize import toplevel_list
 from ..transformations.reification import only_implies, reify_rewrite, only_bv_reifies
-from ..transformations.safening import no_partial_functions
+from ..transformations.safening import no_partial_functions, safen_objective
+
 
 class CPM_cplex(SolverInterface):
     """
@@ -148,6 +150,7 @@ def __init__(self, cpm_model=None, subsolver=None):
 
         from docplex.mp.model import Model
         self.cplex_model = Model()
+        self._obj_offset = 0
         super().__init__(name="cplex", cpm_model=cpm_model)
 
     @property
@@ -299,18 +302,23 @@ def objective(self, expr, minimize=True):
                 technical side note: any constraints created during conversion of the objective
                 are premanently posted to the solver
         """
-        # make objective function non-nested
-        (flat_obj, flat_cons) = flatten_objective(expr)
-        flat_obj = only_positive_bv_wsum(flat_obj)  # remove negboolviews
-        get_variables(flat_obj, collect=self.user_vars)  # add potentially created variables
-        self.add(flat_cons)
+        # save user vars
+        get_variables(expr, self.user_vars)
+
+        # transform objective
+        obj, safe_cons = safen_objective(expr)
+        obj, decomp_cons = decompose_objective(obj, supported={"min", "max", "abs"}, csemap=self._csemap)
+        obj, flat_cons = flatten_objective(obj, csemap=self._csemap)
+        obj, self._obj_offset = only_positive_bv_wsum_const(obj) # remove negboolviews
+
+        self.add(safe_cons + decomp_cons + flat_cons)
 
         # make objective function or variable and post
-        obj = self._make_numexpr(flat_obj)
+        cplex_obj = self._make_numexpr(obj)
         if minimize:
-            self.cplex_model.set_objective('min',obj)
+            self.cplex_model.set_objective('min', cplex_obj)
         else:
-            self.cplex_model.set_objective('max', obj)
+            self.cplex_model.set_objective('max', cplex_obj)
 
     def has_objective(self):
         return self.cplex_model.is_optimized()
@@ -361,7 +369,7 @@ def transform(self, cpm_expr):
         # apply transformations, then post internally
         # expressions have to be linearized to fit in MIP model. See /transformations/linearize
         cpm_cons = toplevel_list(cpm_expr)
-        cpm_cons = no_partial_functions(cpm_cons, safen_toplevel={"mod", "div"})  # linearize expects safe exprs
+        cpm_cons = no_partial_functions(cpm_cons, safen_toplevel={"mod", "div", "element"})  # linearize and decompose expect safe exprs
         supported = {"min", "max", "abs", "alldifferent"} # alldiff has a specialized MIP decomp in linearize
         cpm_cons = decompose_in_tree(cpm_cons, supported, csemap=self._csemap)
         cpm_cons = flatten_constraint(cpm_cons, csemap=self._csemap)  # flat normal form
@@ -600,7 +608,7 @@ def solveAll(self, display=None, time_limit=None, solution_limit=None, call_from
 
                 # Translate objective
                 if self.has_objective():
-                    self.objective_value_ = sol_obj_val
+                    self.objective_value_ = sol_obj_val + self._obj_offset
 
                 if display is not None:
                     if isinstance(display, Expression):
diff --git a/cpmpy/solvers/cpo.py b/cpmpy/solvers/cpo.py
index f20913f34..89b9f6c6c 100644
--- a/cpmpy/solvers/cpo.py
+++ b/cpmpy/solvers/cpo.py
@@ -56,7 +56,7 @@
 from ..expressions.utils import is_num, is_any_list, eval_comparison, argval, argvals, get_bounds
 from ..transformations.get_variables import get_variables
 from ..transformations.normalize import toplevel_list
-from ..transformations.decompose_global import decompose_in_tree
+from ..transformations.decompose_global import decompose_in_tree, decompose_objective
 from ..transformations.safening import no_partial_functions
 
 
@@ -379,14 +379,23 @@ def objective(self, expr, minimize=True):
             
                 technical side note: any constraints created during conversion of the objective are permanently posted to the solver
         """
+
+        # save user variables
+        get_variables(expr, self.user_vars)
+
+        supported = {'max', 'min', 'abs', 'element', 'nvalue', 'alldifferent', 'table', 'indomain', 'negative_table'}
+        obj, decomp_cons = decompose_objective(expr, supported=supported)
+        self.add(decomp_cons)
+
         dom = self.get_docp().modeler
         if self.has_objective():
             self.cpo_model.remove(self.cpo_model.get_objective_expression())
-        expr = self._cpo_expr(expr)
+
+        cpo_obj = self._cpo_expr(obj)
         if minimize:
-            self.cpo_model.add(dom.minimize(expr))
+            self.cpo_model.add(dom.minimize(cpo_obj))
         else:
-            self.cpo_model.add(dom.maximize(expr))
+            self.cpo_model.add(dom.maximize(cpo_obj))
 
     def has_objective(self):
         return self.cpo_model.get_objective() is not None
diff --git a/cpmpy/solvers/exact.py b/cpmpy/solvers/exact.py
index 37c00c1b0..6938f714c 100644
--- a/cpmpy/solvers/exact.py
+++ b/cpmpy/solvers/exact.py
@@ -59,11 +59,11 @@
 from ..transformations.comparison import only_numexpr_equality
 from ..transformations.flatten_model import flatten_constraint, flatten_objective
 from ..transformations.get_variables import get_variables
-from ..transformations.decompose_global import decompose_in_tree
+from ..transformations.decompose_global import decompose_in_tree, decompose_objective
 from ..transformations.linearize import linearize_constraint, only_positive_bv, only_positive_bv_wsum
 from ..transformations.reification import only_implies, reify_rewrite, only_bv_reifies
 from ..transformations.normalize import toplevel_list
-from ..transformations.safening import no_partial_functions
+from ..transformations.safening import no_partial_functions, safen_objective
 from ..expressions.globalconstraints import DirectConstraint
 from ..expressions.utils import flatlist, argvals, argval
 from ..exceptions import NotSupportedError
@@ -436,14 +436,19 @@ def objective(self, expr, minimize):
         self.objective_ = expr
         self.objective_is_min_ = minimize
 
-        # make objective function non-nested and with positive BoolVars only
-        (flat_obj, flat_cons) = flatten_objective(expr, csemap=self._csemap)
-        flat_obj = only_positive_bv_wsum(flat_obj)  # remove negboolviews
-        self.user_vars.update(get_variables(flat_obj))  # add objvars to vars
-        self += flat_cons  # add potentially created constraints
+        # save user variables
+        get_variables(expr, self.user_vars)
+
+        # transform objective
+        obj, safe_cons = safen_objective(expr)
+        obj, decomp_cons = decompose_objective(obj, csemap=self._csemap)
+        obj, flat_cons = flatten_objective(obj, csemap=self._csemap)
+        obj = only_positive_bv_wsum(obj)  # remove negboolviews
+
+        self.add(safe_cons + decomp_cons + flat_cons)
 
         # make objective function or variable and post
-        xct_cfvars,xct_rhs = self._make_numexpr(flat_obj,0)
+        xct_cfvars,xct_rhs = self._make_numexpr(obj,0)
 
         # TODO: make this a custom transformation?
         newcfvrs = []
@@ -504,7 +509,7 @@ def transform(self, cpm_expr):
         """
 
         cpm_cons = toplevel_list(cpm_expr)
-        cpm_cons = no_partial_functions(cpm_cons, safen_toplevel={"mod", "div"}) # linearize expects safe exprs
+        cpm_cons = no_partial_functions(cpm_cons, safen_toplevel={"mod", "div", "element"}) # linearize and decompose expects safe exprs
         cpm_cons = decompose_in_tree(cpm_cons, supported=frozenset({'alldifferent', 'abs'}), csemap=self._csemap) # Abs and Alldiff have a specialized MIP decomp
         cpm_cons = flatten_constraint(cpm_cons, csemap=self._csemap)  # flat normal form
         cpm_cons = reify_rewrite(cpm_cons, supported=frozenset(['sum', 'wsum']), csemap=self._csemap)  # constraints that support reification
diff --git a/cpmpy/solvers/gcs.py b/cpmpy/solvers/gcs.py
index 39321f406..4dff7866d 100644
--- a/cpmpy/solvers/gcs.py
+++ b/cpmpy/solvers/gcs.py
@@ -58,12 +58,12 @@
 from ..exceptions import NotSupportedError, GCSVerificationException
 from .solver_interface import SolverInterface, SolverStatus, ExitStatus
 from ..expressions.core import Expression, Comparison, Operator, BoolVal
-from ..expressions.variables import _BoolVarImpl, _IntVarImpl, _NumVarImpl, NegBoolView, boolvar
+from ..expressions.variables import _BoolVarImpl, _IntVarImpl, _NumVarImpl, NegBoolView, boolvar, intvar
 from ..expressions.globalconstraints import GlobalConstraint
 from ..expressions.utils import is_num, argval, argvals
-from ..transformations.decompose_global import decompose_in_tree
+from ..transformations.decompose_global import decompose_in_tree, decompose_objective
 from ..transformations.get_variables import get_variables
-from ..transformations.flatten_model import flatten_constraint, flatten_objective, get_or_make_var
+from ..transformations.flatten_model import flatten_constraint, get_or_make_var
 from ..transformations.safening import no_partial_functions
 
 from ..transformations.normalize import toplevel_list
@@ -385,20 +385,22 @@ def objective(self, expr, minimize=True):
                 technical side note: any constraints created during conversion of the objective
                 are permanently posted to the solver
         """
-        # make objective function non-nested
-        (flat_obj, flat_cons) = flatten_objective(expr, csemap=self._csemap)
-        self += flat_cons # add potentially created constraints
-        self.user_vars.update(get_variables(flat_obj)) # add objvars to vars
 
-        (obj, obj_cons) = get_or_make_var(flat_obj, csemap=self._csemap)
-        self += obj_cons
+        # save variables
+        get_variables(expr, collect=self.user_vars)
 
-        self.objective_var = obj
+        # transform objective
+        obj, decomp_cons = decompose_objective(expr, csemap=self._csemap,
+                                               supported={"min", "max", "abs", "element", "nvalue", "count"})
+        obj_var, obj_cons = get_or_make_var(obj) # do not pass csemap here, we will still transform obj_var == obj...
+        self.add(decomp_cons + obj_cons)
+
+        self.objective_var = obj_var
 
         if minimize:
-            self.gcs.minimise(self.solver_var(obj))  
+            self.gcs.minimise(self.solver_var(obj_var))
         else:
-            self.gcs.maximise(self.solver_var(obj))
+            self.gcs.maximise(self.solver_var(obj_var))
 
     def transform(self, cpm_expr):
         """
diff --git a/cpmpy/solvers/gurobi.py b/cpmpy/solvers/gurobi.py
index 9d1e6f012..df81baa6f 100644
--- a/cpmpy/solvers/gurobi.py
+++ b/cpmpy/solvers/gurobi.py
@@ -52,13 +52,13 @@
 from ..expressions.variables import _BoolVarImpl, NegBoolView, _IntVarImpl, _NumVarImpl, intvar
 from ..expressions.globalconstraints import DirectConstraint
 from ..transformations.comparison import only_numexpr_equality
-from ..transformations.decompose_global import decompose_in_tree
+from ..transformations.decompose_global import decompose_in_tree, decompose_objective
 from ..transformations.flatten_model import flatten_constraint, flatten_objective
 from ..transformations.get_variables import get_variables
 from ..transformations.linearize import linearize_constraint, only_positive_bv, only_positive_bv_wsum
 from ..transformations.normalize import toplevel_list
 from ..transformations.reification import only_implies, reify_rewrite, only_bv_reifies
-from ..transformations.safening import no_partial_functions
+from ..transformations.safening import no_partial_functions, safen_objective
 
 try:
     import gurobipy as gp
@@ -284,18 +284,23 @@ def objective(self, expr, minimize=True):
         """
         from gurobipy import GRB
 
-        # make objective function non-nested
-        (flat_obj, flat_cons) = flatten_objective(expr, csemap=self._csemap)
-        flat_obj = only_positive_bv_wsum(flat_obj)  # remove negboolviews
-        get_variables(flat_obj, collect=self.user_vars)  # add potentially created variables
-        self += flat_cons
+        # save user variables
+        get_variables(expr, self.user_vars)
+
+        # transform objective
+        obj, safe_cons = safen_objective(expr)
+        obj, decomp_cons = decompose_objective(obj, supported={"min", "max", "abs"}, csemap=self._csemap)
+        obj, flat_cons = flatten_objective(obj, csemap=self._csemap)
+        obj = only_positive_bv_wsum(obj)  # remove negboolviews
+
+        self.add(safe_cons + decomp_cons + flat_cons)
 
         # make objective function or variable and post
-        obj = self._make_numexpr(flat_obj)
+        grb_obj = self._make_numexpr(obj)
         if minimize:
-            self.grb_model.setObjective(obj, sense=GRB.MINIMIZE)
+            self.grb_model.setObjective(grb_obj, sense=GRB.MINIMIZE)
         else:
-            self.grb_model.setObjective(obj, sense=GRB.MAXIMIZE)
+            self.grb_model.setObjective(grb_obj, sense=GRB.MAXIMIZE)
         self.grb_model.update()
 
     def has_objective(self):
@@ -347,7 +352,7 @@ def transform(self, cpm_expr):
         # apply transformations, then post internally
         # expressions have to be linearized to fit in MIP model. See /transformations/linearize
         cpm_cons = toplevel_list(cpm_expr)
-        cpm_cons = no_partial_functions(cpm_cons, safen_toplevel={"mod", "div"})  # linearize expects safe exprs
+        cpm_cons = no_partial_functions(cpm_cons, safen_toplevel={"mod", "div", "element"})  # linearize and decompose expect safe exprs
         supported = {"min", "max", "abs", "alldifferent"} # alldiff has a specialized MIP decomp in linearize
         cpm_cons = decompose_in_tree(cpm_cons, supported, csemap=self._csemap)
         cpm_cons = flatten_constraint(cpm_cons, csemap=self._csemap)  # flat normal form
diff --git a/cpmpy/solvers/minizinc.py b/cpmpy/solvers/minizinc.py
index d1227ec91..02cd224b0 100644
--- a/cpmpy/solvers/minizinc.py
+++ b/cpmpy/solvers/minizinc.py
@@ -71,7 +71,7 @@
 from ..expressions.variables import _NumVarImpl, _IntVarImpl, _BoolVarImpl, NegBoolView, cpm_array
 from ..expressions.globalconstraints import DirectConstraint
 from ..expressions.utils import is_num, is_any_list, argvals, argval
-from ..transformations.decompose_global import decompose_in_tree
+from ..transformations.decompose_global import decompose_in_tree, decompose_objective
 from ..exceptions import MinizincPathException, NotSupportedError
 from ..transformations.get_variables import get_variables
 from ..transformations.normalize import toplevel_list
@@ -555,15 +555,27 @@ def objective(self, expr, minimize):
 
             'objective()' can be called multiple times, only the last one is stored
         """
-        # get_variables(expr, collect=self.user_vars)  # add objvars to vars  # all are user vars
+
+        # save user variables
+        get_variables(expr, collect=self.user_vars) # add objvars to vars
+
+        supported = {"min", "max", "abs", "element", "count", "nvalue", "alldifferent", "alldifferent_except0", "allequal",
+                     "inverse", "ite" "xor", "table", "cumulative", "gcc", "increasing", "decreasing",
+                     "no_overlap", "strictly_increasing", "strictly_decreasing", "lex_lesseq", "lex_less", "lex_chain_less",
+                     "lex_chain_lesseq", "among"}
+
+        # transform objective
+        obj, decomp_cons = decompose_objective(expr, supported=supported)
+        self.add(decomp_cons)
 
         # make objective function or variable and post
-        obj = self._convert_expression(expr)
+
+        mzn_obj = self._convert_expression(obj)
         # do not add it to the mzn_model yet, supports only one 'solve' entry
         if minimize:
-            self.mzn_txt_solve = "solve minimize {};\n".format(obj)
+            self.mzn_txt_solve = "solve minimize {};\n".format(mzn_obj)
         else:
-            self.mzn_txt_solve = "solve maximize {};\n".format(obj)
+            self.mzn_txt_solve = "solve maximize {};\n".format(mzn_obj)
 
     def has_objective(self):
         return self.mzn_txt_solve != "solve satisfy;"
diff --git a/cpmpy/solvers/ortools.py b/cpmpy/solvers/ortools.py
index afc68026f..6f28895fd 100644
--- a/cpmpy/solvers/ortools.py
+++ b/cpmpy/solvers/ortools.py
@@ -56,13 +56,13 @@
 from ..expressions.globalconstraints import GlobalConstraint
 from ..expressions.utils import is_num, is_int, eval_comparison, flatlist, argval, argvals, get_bounds, is_true_cst, \
     is_false_cst
-from ..transformations.decompose_global import decompose_in_tree
+from ..transformations.decompose_global import decompose_in_tree, decompose_objective
 from ..transformations.get_variables import get_variables
 from ..transformations.flatten_model import flatten_constraint, flatten_objective, get_or_make_var
 from ..transformations.normalize import toplevel_list
 from ..transformations.reification import only_implies, reify_rewrite, only_bv_reifies
 from ..transformations.comparison import only_numexpr_equality
-from ..transformations.safening import no_partial_functions
+from ..transformations.safening import no_partial_functions, safen_objective
 
 
 class CPM_ortools(SolverInterface):
@@ -339,17 +339,23 @@ def objective(self, expr, minimize):
                 technical side note: any constraints created during conversion of the objective
                 are premanently posted to the solver
         """
-        # make objective function non-nested
-        (flat_obj, flat_cons) = flatten_objective(expr, csemap=self._csemap)
-        self += flat_cons  # add potentially created constraints
-        get_variables(flat_obj, collect=self.user_vars)  # add objvars to vars
+
+        # save user varables
+        get_variables(expr, self.user_vars)
+
+        # transform objective
+        obj, safe_cons = safen_objective(expr)
+        obj, decomp_cons = decompose_objective(obj, supported={"min", "max", "abs", "element"}, csemap=self._csemap)
+        obj, flat_cons = flatten_objective(obj, csemap=self._csemap)
+
+        self.add(safe_cons+decomp_cons+flat_cons)
 
         # make objective function or variable and post
-        obj = self._make_numexpr(flat_obj)
+        ort_obj = self._make_numexpr(obj)
         if minimize:
-            self.ort_model.Minimize(obj)
+            self.ort_model.Minimize(ort_obj)
         else:
-            self.ort_model.Maximize(obj)
+            self.ort_model.Maximize(ort_obj)
 
     def has_objective(self):
         return self.ort_model.HasObjective()
diff --git a/cpmpy/solvers/pindakaas.py b/cpmpy/solvers/pindakaas.py
index 20540903a..64189319e 100755
--- a/cpmpy/solvers/pindakaas.py
+++ b/cpmpy/solvers/pindakaas.py
@@ -231,7 +231,7 @@ def solver_var(self, cpm_var):
 
     def transform(self, cpm_expr):
         cpm_cons = toplevel_list(cpm_expr)
-        cpm_cons = no_partial_functions(cpm_cons)
+        cpm_cons = no_partial_functions(cpm_cons, safen_toplevel={"div", "mod", "element"})
         cpm_cons = decompose_in_tree(cpm_cons, csemap=self._csemap)
         cpm_cons = simplify_boolean(cpm_cons)
         cpm_cons = flatten_constraint(cpm_cons, csemap=self._csemap)  # flat normal form
diff --git a/cpmpy/solvers/pumpkin.py b/cpmpy/solvers/pumpkin.py
index f1f5d0034..31bd1d35f 100644
--- a/cpmpy/solvers/pumpkin.py
+++ b/cpmpy/solvers/pumpkin.py
@@ -52,8 +52,8 @@
 from ..transformations.get_variables import get_variables
 from ..transformations.linearize import canonical_comparison
 from ..transformations.normalize import toplevel_list
-from ..transformations.decompose_global import decompose_in_tree
-from ..transformations.flatten_model import flatten_constraint
+from ..transformations.decompose_global import decompose_in_tree, decompose_objective
+from ..transformations.flatten_model import flatten_constraint, get_or_make_var
 from ..transformations.comparison import only_numexpr_equality
 from ..transformations.reification import reify_rewrite, only_bv_reifies, only_implies
 from ..transformations.safening import no_partial_functions
@@ -290,11 +290,21 @@ def objective(self, expr, minimize=True):
                 technical side note: any constraints created during conversion of the objective
                 are premanently posted to the solver
         """
-        # make objective function non-nested
-        obj_var = intvar(*get_bounds(expr))
-        self += expr == obj_var
 
-        # make objective function or variable and post
+        # save user variables
+        get_variables(expr, self.user_vars)
+
+        # transform objective
+        obj, decomp_cons = decompose_objective(expr, supported={"min", "max", "element", "abs"}, csemap=self._csemap)
+        obj_var, obj_cons = get_or_make_var(obj) # do not pass csemap here, we will still transform obj_var == obj...
+        if expr.is_bool():
+            ivar = intvar(0,1)
+            obj_cons += [ivar == obj_var]
+            obj_var = ivar
+
+        self.add(decomp_cons + obj_cons)
+
+        # save objective function
         self._objective = obj_var
         self.objective_is_min = minimize
 
diff --git a/cpmpy/solvers/pysat.py b/cpmpy/solvers/pysat.py
index def9e6804..835a9281c 100644
--- a/cpmpy/solvers/pysat.py
+++ b/cpmpy/solvers/pysat.py
@@ -70,6 +70,7 @@
 from ..transformations.normalize import toplevel_list, simplify_boolean
 from ..transformations.reification import only_implies, only_bv_reifies
 from ..transformations.int2bool import int2bool, _encode_int_var, _decide_encoding
+from ..transformations.safening import no_partial_functions
 
 
 class CPM_pysat(SolverInterface):
@@ -366,6 +367,7 @@ def transform(self, cpm_expr):
             :return: list of Expression
         """
         cpm_cons = toplevel_list(cpm_expr)
+        cpm_cons = no_partial_functions(cpm_cons, safen_toplevel={"div", "mod", "element"})
         cpm_cons = decompose_in_tree(cpm_cons, supported=frozenset({"alldifferent"}), csemap=self._csemap)
         cpm_cons = simplify_boolean(cpm_cons)
         cpm_cons = flatten_constraint(cpm_cons, csemap=self._csemap)  # flat normal form
diff --git a/cpmpy/solvers/z3.py b/cpmpy/solvers/z3.py
index 3451583d3..d5a003f8b 100644
--- a/cpmpy/solvers/z3.py
+++ b/cpmpy/solvers/z3.py
@@ -56,9 +56,9 @@
 from ..expressions.globalfunctions import GlobalFunction
 from ..expressions.variables import _BoolVarImpl, NegBoolView, _NumVarImpl, _IntVarImpl, intvar
 from ..expressions.utils import is_num, is_any_list, is_bool, is_int, is_boolexpr, eval_comparison
-from ..transformations.decompose_global import decompose_in_tree
+from ..transformations.decompose_global import decompose_in_tree, decompose_objective
 from ..transformations.normalize import toplevel_list
-from ..transformations.safening import no_partial_functions
+from ..transformations.safening import no_partial_functions, safen_objective
 
 
 class CPM_z3(SolverInterface):
@@ -312,17 +312,20 @@ def objective(self, expr, minimize=True):
         if not isinstance(self.z3_solver, z3.Optimize):
             raise NotSupportedError("Use the z3 optimizer for optimization problems")
 
-        if isinstance(expr, GlobalFunction): # not supported by Z3
-            obj_var = intvar(*expr.get_bounds())
-            self += expr == obj_var
-            expr = obj_var
+        # transform objective
+        obj, safe_cons = safen_objective(expr)
+        obj, decomp_cons = decompose_objective(obj, csemap=self._csemap)
 
-        obj = self._z3_expr(expr)
+        self.add(safe_cons + decomp_cons)
+
+        z3_obj = self._z3_expr(obj)
+        if isinstance(z3_obj, z3.BoolRef):
+            z3_obj = z3.If(z3_obj, 1, 0) # must be integer
         if minimize:
-            self.obj_handle = self.z3_solver.minimize(obj)
+            self.obj_handle = self.z3_solver.minimize(z3_obj)
             self._minimize = True # record direction of optimisation
         else:
-            self.obj_handle = self.z3_solver.maximize(obj)
+            self.obj_handle = self.z3_solver.maximize(z3_obj)
             self._minimize = False # record direction of optimisation
 
 
@@ -342,7 +345,7 @@ def transform(self, cpm_expr):
         """
 
         cpm_cons = toplevel_list(cpm_expr)
-        cpm_cons = no_partial_functions(cpm_cons, safen_toplevel={"div", "mod"})
+        cpm_cons = no_partial_functions(cpm_cons, safen_toplevel={"div", "mod", "element"})
         supported = {"alldifferent", "xor", "ite"}  # z3 accepts these reified too
         cpm_cons = decompose_in_tree(cpm_cons, supported, supported, csemap=self._csemap)
         return cpm_cons
diff --git a/cpmpy/transformations/decompose_global.py b/cpmpy/transformations/decompose_global.py
index 0b752080f..328217aa9 100644
--- a/cpmpy/transformations/decompose_global.py
+++ b/cpmpy/transformations/decompose_global.py
@@ -11,7 +11,7 @@
 from ..expressions.core import Expression, Comparison, Operator
 from ..expressions.variables import intvar, cpm_array, NDVarArray
 from ..expressions.utils import is_any_list, eval_comparison
-from ..expressions.python_builtins import all
+from ..expressions.python_builtins import all as cpm_all
 from .flatten_model import flatten_constraint, normalized_numexpr
 
 
@@ -36,9 +36,6 @@ def decompose_in_tree(lst_of_expr, supported=set(), supported_reified=set(), _to
     if _toplevel is None:
         _toplevel = []
 
-    # swap the arguments of a comparison while maintaining its semantics
-    flipmap = {"==": "==", "!=": "!=", "<": ">", "<=": ">=", ">": "<", ">=": "<="}
-
     newlist = []  # decomposed constraints will go here
     for expr in lst_of_expr:
 
@@ -63,9 +60,9 @@ def decompose_in_tree(lst_of_expr, supported=set(), supported_reified=set(), _to
                 newlist.append(expr)
                 continue
 
-            if any(isinstance(a,GlobalFunction) for a in expr.args):
-                expr, base_con = normalized_numexpr(expr, csemap=csemap)
-                _toplevel.extend(base_con)  # should be added toplevel
+            # if any(isinstance(a,GlobalFunction) for a in expr.args):
+            #     expr, base_con = normalized_numexpr(expr, csemap=csemap)
+            #     _toplevel.extend(base_con)  # should be added toplevel
             # recurse into arguments, recreate through constructor (we know it stores no other state)
             args = decompose_in_tree(expr.args, supported, supported_reified, _toplevel, nested=True, csemap=csemap)
             newlist.append(Operator(expr.name, args))
@@ -89,93 +86,45 @@ def decompose_in_tree(lst_of_expr, supported=set(), supported_reified=set(), _to
                     expr.update_args(decompose_in_tree(expr.args, supported, supported_reified, _toplevel, nested=True, csemap=csemap))
                     newlist.append(expr)
 
-            else:
+            else: # unsupported, need to decompose
+                # check if it is in the csemap
+                if csemap is not None and expr in csemap:
+                    newlist.append(csemap[expr])
+                    continue # no need to decompose, just use the variable we already have
+
+                dec = expr.decompose()
+                if not isinstance(dec, tuple):
+                    warnings.warn(f"Decomposing an old-style global ({expr}) that does not return a tuple, which is "
+                                  "deprecated. Support for old-style globals will be removed in stable version",
+                                  DeprecationWarning)
+                    dec = (dec, [])
+                value, define = dec
+                if not is_any_list(value):
+                    value = [value]
+
+                _toplevel.extend(define)  # definitions should be added toplevel
+                # the `decomposed` expression might contain other global constraints, check it
+                decomposed = decompose_in_tree(value, supported, supported_reified, _toplevel, nested=nested, csemap=csemap)
                 if expr.is_bool():
-                    assert isinstance(expr, GlobalConstraint)
-                    # boolean global constraints
-                    dec = expr.decompose()
-                    if not isinstance(dec, tuple):
-                        warnings.warn(f"Decomposing an old-style global ({expr}) that does not return a tuple, which is "
-                                      "deprecated. Support for old-style globals will be removed in stable version",
-                                      DeprecationWarning)
-                        dec = (dec, [])
-                    decomposed, define = dec
-
-                    _toplevel.extend(define)  # definitions should be added toplevel
-                    # the `decomposed` expression might contain other global constraints, check it
-                    decomposed = decompose_in_tree(decomposed, supported, supported_reified, _toplevel, nested=nested, csemap=csemap)
-                    newlist.append(all(decomposed))
+                    value = cpm_all(decomposed)
 
                 else:
-                    # global function, replace by a fresh variable and decompose the equality to this
-                    assert isinstance(expr, GlobalFunction)
-                    lb,ub = expr.get_bounds()
-                    
-                    if csemap is not None and expr in csemap:
-                        aux = csemap[expr]
-                    else:
-                        aux = intvar(lb, ub)
-                        if csemap is not None:
-                            csemap[expr] = aux
-
-                    # NOTE Do we need to decompose here (the expr's args)? As done in the Comparison's section?
-
-                    dec = expr.decompose_comparison("==", aux)
-                    if not isinstance(dec, tuple):
-                        warnings.warn(f"Decomposing an old-style global ({expr}) that does not return a tuple, which is "
-                                      "deprecated. Support for old-style globals will be removed in stable version",
-                                      DeprecationWarning)
-                        dec = (dec, [])
-                    auxdef, otherdef = dec
-
-                    _toplevel.extend(auxdef + otherdef)  # all definitions should be added toplevel
-                    newlist.append(aux)  # replace original expression by aux
+                    assert len(decomposed) == 1, "Global functions should return a single numerical value, not a list"
+                    value = decomposed[0]
+
+                newlist.append(value)
+                if csemap is not None:
+                    csemap[expr] = value
+
 
         elif isinstance(expr, Comparison):
             if not expr.has_subexpr(): # Only recurse if there are nested expressions
                 newlist.append(expr)
                 continue
 
-            # if one of the two children is a (numeric) global constraint, we can decompose the comparison directly
-            # otherwise e.g., min(x,y,z) == a would become `min(x,y,z).decompose_comparison('==',aux) + [aux == a]`
-            lhs, rhs = expr.args
-            exprname = expr.name  # can change when rhs needs decomp
-
-            decomp_lhs = isinstance(lhs, GlobalFunction) and not lhs.name in supported
-            decomp_rhs = isinstance(rhs, GlobalFunction) and not rhs.name in supported
-
-            if not decomp_lhs:
-                if not decomp_rhs:
-                    # nothing special, create a fresh version and recurse into arguments
-                    expr = copy.copy(expr)
-                    expr.update_args(decompose_in_tree(expr.args, supported, supported_reified, _toplevel, nested=True, csemap=csemap))
-                    newlist.append(expr)
-
-                else:
-                    # only rhs needs decomposition, so flip comparison to make lhs needing the decomposition
-                    exprname = flipmap[expr.name]
-                    lhs, rhs = rhs, lhs
-                    decomp_lhs, decomp_rhs = True, False  # continue into next 'if'
-
-            if decomp_lhs:
-                if lhs.has_subexpr():
-                    # recurse into lhs args and decompose nested subexpressions
-                    lhs = copy.copy(lhs)
-                    lhs.update_args(decompose_in_tree(lhs.args, supported, supported_reified, _toplevel, nested=True, csemap=csemap))
-
-                # decompose comparison of lhs and rhs
-                dec = lhs.decompose_comparison(exprname, rhs)
-                if not isinstance(dec, tuple):
-                    warnings.warn(f"Decomposing an old-style global ({lhs}) that does not return a tuple, which is "
-                                  f"deprecated. Support for old-style globals will be removed in stable version",
-                                  DeprecationWarning)
-                    dec = (dec, [])
-                decomposed, define = dec
-
-                _toplevel.extend(define)  # definitions should be added toplevel
-                # the `decomposed` expression (and rhs) might contain other global constraints, check it
-                decomposed = decompose_in_tree(decomposed, supported, supported_reified, _toplevel, nested=True, csemap=csemap)
-                newlist.append(all(decomposed))
+            expr = copy.copy(expr)
+            expr.update_args(decompose_in_tree(expr.args, supported, supported_reified, _toplevel, nested=True, csemap=csemap))
+            newlist.append(expr)
 
         else:  # constants, variables, direct constraints
             newlist.append(expr)
@@ -190,6 +139,16 @@ def decompose_in_tree(lst_of_expr, supported=set(), supported_reified=set(), _to
         return toplevel_list(newlist) + decompose_in_tree(toplevel_list(_toplevel), supported, supported_reified, nested=False, csemap=csemap)
 
 
+def decompose_objective(expr, supported=set(), csemap=None):
+    if is_any_list(expr):
+        raise ValueError(f"Expected a numerical expression as objective but got a list {expr}")
+
+    toplevel = []
+    decomp_expr = decompose_in_tree([expr], supported=supported, _toplevel=toplevel, nested=True, csemap=csemap)
+    assert len(decomp_expr) == 1, f"Expected {expr} to be decomposed into a single expression, but got {decomp_expr}.\nPlease report on github."
+    return decomp_expr[0], toplevel
+
+
 # DEPRECATED!
 # old way of doing decompositions post-flatten
 # will be removed in any future version!
diff --git a/cpmpy/transformations/safening.py b/cpmpy/transformations/safening.py
index f2e9e0f1d..626bcca7c 100644
--- a/cpmpy/transformations/safening.py
+++ b/cpmpy/transformations/safening.py
@@ -6,7 +6,7 @@
 
 from ..expressions.variables import _NumVarImpl, boolvar, intvar, NDVarArray, cpm_array
 from ..expressions.core import Expression, Operator, BoolVal
-from ..expressions.utils import get_bounds, is_num
+from ..expressions.utils import get_bounds, is_num, is_any_list
 from ..expressions.globalfunctions import GlobalFunction, Element
 from ..expressions.globalconstraints import DirectConstraint
 from ..expressions.python_builtins import all as cpm_all
@@ -249,3 +249,12 @@ def _safen_hole(cpm_expr, exclude, idx_to_safen):
     return is_defined, output_var, toplevel
 
 
+def safen_objective(expr):
+    if is_any_list(expr):
+        raise ValueError(f"Expected numerical expression as objective but got a list {expr}")
+
+    toplevel, nbc = [],[]
+    safe_expr = no_partial_functions([expr], _toplevel=toplevel, _nbc=nbc)
+    assert len(safe_expr) == 1
+    return safe_expr[0], toplevel + nbc
+
diff --git a/tests/test_cse.py b/tests/test_cse.py
index 125dfd7c0..35c5ca091 100644
--- a/tests/test_cse.py
+++ b/tests/test_cse.py
@@ -70,23 +70,29 @@ def test_decompose(self):
         csemap = dict()
         decomp = decompose_in_tree([nested_cons], csemap=csemap)
     
-        self.assertEqual(len(decomp), 6)
+        self.assertEqual(len(decomp), 5)
         self.assertEqual(str(decomp[0]), "(b) == ((IV0) + (q) <= 10)")
-        self.assertEqual(str(decomp[1]), "(IV1) == (IV0)") # TODO... this seems stupid, why do we need this (comes from _max in Maximul decomp)?
-        self.assertEqual(str(decomp[2]), "or([(x) >= (IV1), (y) >= (IV1), (z) >= (IV1)])")
-        self.assertEqual(str(decomp[3]), "(x) <= (IV1)")
-        self.assertEqual(str(decomp[4]), "(y) <= (IV1)")
-        self.assertEqual(str(decomp[5]), "(z) <= (IV1)")
+        self.assertEqual(str(decomp[1]), "(x) <= (IV0)")
+        self.assertEqual(str(decomp[2]), "(y) <= (IV0)")
+        self.assertEqual(str(decomp[3]), "(z) <= (IV0)")
+        self.assertEqual(str(decomp[4]), "or([(x) >= (IV0), (y) >= (IV0), (z) >= (IV0)])")
+
 
         # next time we use max([x,y,z]) it should replace the max-constraint with IV0
-        #  ... it seems like we should be able to do more here e.g., cp.max([x,y,z]) != 42 should be replaced with IV0 != 42
-        #  ...      but the current code-flow of decompose_in_tree and .decompose_comparison does not allow this
         nested2 = (q + cp.max([x,y,z]) != 42)
         decomp = decompose_in_tree([nested2], csemap=csemap)
 
         self.assertEqual(len(decomp), 1)
         self.assertEqual(str(decomp[0]), "(q) + (IV0) != 42")
 
+        # also in non-nested cases
+        nested3 = (cp.max([x, y, z]) == 42)
+        decomp = decompose_in_tree([nested3], csemap=csemap)
+
+        self.assertEqual(len(decomp), 1)
+        self.assertEqual(str(decomp[0]), "IV0 == 42")
+
+
     def test_only_numexpr_equality(self):
         x,y,z = cp.intvar(0,10, shape=3, name=tuple("xyz"))
 
diff --git a/tests/test_globalconstraints.py b/tests/test_globalconstraints.py
index 2b8c8e20e..6e9c171e4 100644
--- a/tests/test_globalconstraints.py
+++ b/tests/test_globalconstraints.py
@@ -8,6 +8,8 @@
 from cpmpy.exceptions import TypeError, NotSupportedError
 from cpmpy.expressions.utils import STAR, argvals
 from cpmpy.solvers import CPM_minizinc
+from cpmpy.transformations.decompose_global import decompose_in_tree
+from cpmpy.transformations.safening import no_partial_functions
 
 from utils import skip_on_missing_pblib
 
@@ -1444,13 +1446,23 @@ def test_issue_699(self):
         self.assertTrue(cp.Model(cp.AllDifferentExceptN([x,3,y,0], [3,0]).decompose()).solve())
 
 
-    def test_element_index_dom_mismatched(self):
-        """
-            Check transform of `[0,1,2][x in -1..1] == y in 1..5`
-            Note the index variable has a lower bound *outside* the indexable range, and an upper bound inside AND lower than the indexable range upper bound
-        """
-        constraint=cp.Element([0,1,2], cp.intvar(-1,1, name="x"))
-        self.assertEqual(
-            str(constraint.decompose_comparison("==", cp.intvar(1,5, name="y"))),
-            "([(x == 0) -> (y == 0), (x == 1) -> (y == 1), x >= 0, x <= 1], [])"
-        )
+    # disabled after #793, unsafe element can no longer be decomposed
+    # def test_element_index_dom_mismatched(self):
+    #     """
+    #         Check transform of `[0,1,2][x in -1..1] == y in 1..5`
+    #         Note the index variable has a lower bound *outside* the indexable range, and an upper bound inside AND lower than the indexable range upper bound
+    #     """
+    #     constraint=cp.Element([0,1,2], cp.intvar(-1,1, name="x")) <= cp.intvar(1,5, name="y")
+    #     decomposed = decompose_in_tree(no_partial_functions([constraint], safen_toplevel={"element"}))
+    #     self.assertSetEqual(set(map(str, decomposed)), {
+    #         # safening constraints
+    #         "(BV0) == ((x >= 0) and (x <= 2))",
+    #         "(BV0) -> ((IV0) == (x))",
+    #         "(~BV0) -> (IV0 == 0)",
+    #         "BV0",
+    #         # actual decomposition
+    #         '(IV0 == 0) -> (IV1 == 0)',
+    #         '(IV0 == 1) -> (IV1 == 1)',
+    #         '(IV0 == 2) -> (IV1 == 2)',
+    #         '(IV1) <= (y)'
+    #     })
\ No newline at end of file
diff --git a/tests/test_solvers.py b/tests/test_solvers.py
index b0416e7e9..c687fcce2 100644
--- a/tests/test_solvers.py
+++ b/tests/test_solvers.py
@@ -21,6 +21,7 @@
 from cpmpy import SolverLookup
 from cpmpy.exceptions import MinizincNameException, NotSupportedError
 
+from test_constraints import numexprs
 from utils import skip_on_missing_pblib
 
 pysat_available = CPM_pysat.supported()
@@ -821,10 +822,9 @@ def test_count_mzn(self):
         m = cp.Model([x + y == 2, wsum == 9])
         self.assertTrue(m.solve(solver="minizinc"))
 
-@pytest.mark.parametrize(
-        "solver",
-        [name for name, solver in SolverLookup.base_solvers() if solver.supported()]
-)
+
+solvers = [name for name, solver in SolverLookup.base_solvers() if solver.supported()]
+@pytest.mark.parametrize("solver", solvers)
 class TestSupportedSolvers:
     def test_installed_solvers(self, solver):
         # basic model
@@ -1192,4 +1192,22 @@ def test_optimisation_direction(self, solver):
         # 4) Minimisation - solver
         s = cp.SolverLookup.get(solver, m)
         assert s.solve()
-        assert s.objective_value() == 5
\ No newline at end of file
+        assert s.objective_value() == 5
+
+
+
+@pytest.mark.parametrize(("solver", "expr"), [(s, expr) for s in solvers for expr in numexprs(s)], ids=str)
+def test_objective_numexprs(solver, expr):
+
+    if solver == "z3": pytest.skip("Optimization in z3 is too slow for testing here.")
+
+    model = cp.Model(cp.intvar(0, 10, shape=3) >= 1) # just to have some constraints
+    try:
+        model.minimize(expr)
+        assert model.solve(solver=solver)
+        assert expr.value() < expr.get_bounds()[1] # bounds are not always tight, but should be smaller than ub for sure
+        model.maximize(expr)
+        assert model.solve(solver=solver)
+        assert expr.value() > expr.get_bounds()[0] # bounds are not always tight, but should be larger than lb for sure
+    except NotSupportedError:
+        pytest.skip(reason=f"{solver} does not support optimisation")
diff --git a/tests/test_transf_decompose.py b/tests/test_transf_decompose.py
index afe9c2b83..74673b8ea 100644
--- a/tests/test_transf_decompose.py
+++ b/tests/test_transf_decompose.py
@@ -55,22 +55,22 @@ def test_decompose_num(self):
         bv = boolvar(name="bv")
 
         cons = [min(ivs) <= 1]
-        self.assertEqual(str(decompose_in_tree(cons)),
-                         "[IV0 <= 1, ((x) <= (IV0)) or ((y) <= (IV0)), (x) >= (IV0), (y) >= (IV0)]")
+        self.assertSetEqual(set(map(str,decompose_in_tree(cons))),
+                            {"IV0 <= 1", "((x) <= (IV0)) or ((y) <= (IV0))", "(x) >= (IV0)", "(y) >= (IV0)"})
         # reified
         cons = [bv.implies(min(ivs) <= 1)]
-        self.assertEqual(str(decompose_in_tree(cons)),
-                         "[(bv) -> (IV1 <= 1), ((x) <= (IV1)) or ((y) <= (IV1)), (x) >= (IV1), (y) >= (IV1)]")
+        self.assertSetEqual(set(map(str,decompose_in_tree(cons))),
+                            {"(bv) -> (IV1 <= 1)", "((x) <= (IV1)) or ((y) <= (IV1))", "(x) >= (IV1)", "(y) >= (IV1)"})
         self.assertEqual(str(decompose_in_tree(cons, supported={"min"})),str(cons))
 
         cons = [(min(ivs) <= 1).implies(bv)]
-        self.assertEqual(str(decompose_in_tree(cons)),
-                         "[(IV2 <= 1) -> (bv), ((x) <= (IV2)) or ((y) <= (IV2)), (x) >= (IV2), (y) >= (IV2)]")
+        self.assertSetEqual(set(map(str,decompose_in_tree(cons))),
+                            {"(IV2 <= 1) -> (bv)", "((x) <= (IV2)) or ((y) <= (IV2))", "(x) >= (IV2)", "(y) >= (IV2)"})
         self.assertEqual(str(decompose_in_tree(cons, supported={"min"})), str(cons))
 
         cons = [(min(ivs) <= 1) == (bv)]
-        self.assertEqual(str(decompose_in_tree(cons)),
-                         "[(IV3 <= 1) == (bv), ((x) <= (IV3)) or ((y) <= (IV3)), (x) >= (IV3), (y) >= (IV3)]")
+        self.assertSetEqual(set(map(str,decompose_in_tree(cons))),
+                            {"(IV3 <= 1) == (bv)", "((x) <= (IV3)) or ((y) <= (IV3))", "(x) >= (IV3)", "(y) >= (IV3)"})
         self.assertEqual(str(decompose_in_tree(cons, supported={"min"})), str(cons))
 
 
@@ -79,25 +79,25 @@ def test_decompose_nested(self):
         ivs = [intvar(1,9,name=n) for n in "xyz"]
 
         cons = [AllDifferent(ivs) == 0]
-        self.assertEqual(str(decompose_in_tree(cons)), "[not([and([(x) != (y), (x) != (z), (y) != (z)])])]")
+        self.assertSetEqual(set(map(str,decompose_in_tree(cons))), {"not([and([(x) != (y), (x) != (z), (y) != (z)])])"})
 
         cons = [0 == AllDifferent(ivs)]
-        self.assertEqual(str(decompose_in_tree(cons)), "[not([and([(x) != (y), (x) != (z), (y) != (z)])])]")
+        self.assertSetEqual(set(map(str,decompose_in_tree(cons))), {"not([and([(x) != (y), (x) != (z), (y) != (z)])])"})
 
         cons = [AllDifferent(ivs) == AllEqual(ivs[:-1])]
-        self.assertEqual(str(decompose_in_tree(cons)), "[(and([(x) != (y), (x) != (z), (y) != (z)])) == ((x) == (y))]")
+        self.assertSetEqual(set(map(str,decompose_in_tree(cons))), {"(and([(x) != (y), (x) != (z), (y) != (z)])) == ((x) == (y))"})
 
         cons = [min(ivs) == max(ivs)]
-        self.assertEqual(str(decompose_in_tree(cons, supported={"min"})),
-                         '[(IV0) == (min(x,y,z)), or([(x) >= (IV0), (y) >= (IV0), (z) >= (IV0)]), (x) <= (IV0), (y) <= (IV0), (z) <= (IV0)]')
+        self.assertSetEqual(set(map(str,decompose_in_tree(cons, supported={"min"}))),
+                            {"(min(x,y,z)) == (IV0)", "or([(x) >= (IV0), (y) >= (IV0), (z) >= (IV0)])", "(x) <= (IV0)", "(y) <= (IV0)", "(z) <= (IV0)"})
 
-        self.assertEqual(str(decompose_in_tree(cons, supported={"max"})),
-                         "[(IV1) == (max(x,y,z)), or([(x) <= (IV1), (y) <= (IV1), (z) <= (IV1)]), (x) >= (IV1), (y) >= (IV1), (z) >= (IV1)]")
+        self.assertEqual(set(map(str,decompose_in_tree(cons, supported={"max"}))),
+                         {"(IV1) == (max(x,y,z))", "or([(x) <= (IV1), (y) <= (IV1), (z) <= (IV1)])", "(x) >= (IV1)", "(y) >= (IV1)", "(z) >= (IV1)"})
 
         # numerical in non-comparison context
         cons = [AllEqual([min(ivs[:-1]),ivs[-1]])]
-        self.assertEqual(str(decompose_in_tree(cons, supported={"allequal"})),
-                         "[allequal(IV2,z), (IV3) == (IV2), ((x) <= (IV3)) or ((y) <= (IV3)), (x) >= (IV3), (y) >= (IV3)]")
+        self.assertEqual(set(map(str,decompose_in_tree(cons, supported={"allequal"}))),
+                         {"allequal(IV2,z)", "((x) <= (IV2)) or ((y) <= (IV2))", "(x) >= (IV2)", "(y) >= (IV2)"})
 
         self.assertEqual(str(decompose_in_tree(cons, supported={"min"})),
                          "[(min(x,y)) == (z)]")