diff --git a/pypfopt/base_optimizer.py b/pypfopt/base_optimizer.py
index 3af4a15..fb97e04 100644
--- a/pypfopt/base_optimizer.py
+++ b/pypfopt/base_optimizer.py
@@ -9,10 +9,14 @@ evaluate return and risk for a given set of portfolio weights.
 import collections
 import json
 import warnings
+from collections.abc import Iterable
+from typing import List
+
 import numpy as np
 import pandas as pd
 import cvxpy as cp
 import scipy.optimize as sco
+
 from . import objective_functions
 from . import exceptions
 
@@ -174,12 +178,11 @@ class BaseConvexOptimizer(BaseOptimizer):
         self._constraints = []
         self._lower_bounds = None
         self._upper_bounds = None
-        self._map_bounds_to_constraints(weight_bounds)
-
         self._opt = None
         self._solver = solver
         self._verbose = verbose
         self._solver_options = solver_options if solver_options else {}
+        self._map_bounds_to_constraints(weight_bounds)
 
     def _map_bounds_to_constraints(self, test_bounds):
         """
@@ -218,8 +221,32 @@ class BaseConvexOptimizer(BaseOptimizer):
                 self._lower_bounds = np.nan_to_num(lower, nan=-1)
                 self._upper_bounds = np.nan_to_num(upper, nan=1)
 
-        self._constraints.append(self._w >= self._lower_bounds)
-        self._constraints.append(self._w <= self._upper_bounds)
+        self.add_constraint(lambda w: w >= self._lower_bounds)
+        self.add_constraint(lambda w: w <= self._upper_bounds)
+
+    def is_parameter_defined(self, parameter_name: str) -> bool:
+        is_defined = False
+        objective_and_constraints = self._constraints + [self._objective] if self._objective is not None else self._constraints
+        for expr in objective_and_constraints:
+            params = [arg for arg in get_all_args(expr) if isinstance(arg, cp.Parameter)]
+            for param in params:
+                if param.name() == parameter_name and not is_defined:
+                    is_defined = True
+                elif param.name() == parameter_name and is_defined:
+                    raise Exception('Parameter name defined multiple times')
+        return is_defined
+
+    def update_parameter_value(self, parameter_name: str, new_value: float) -> None:
+        assert self.is_parameter_defined(parameter_name)
+        was_updated = False
+        objective_and_constraints = self._constraints + [self._objective] if self._objective is not None else self._constraints
+        for expr in objective_and_constraints:
+            params = [arg for arg in get_all_args(expr) if isinstance(arg, cp.Parameter)]
+            for param in params:
+                if param.name() == parameter_name:
+                    param.value = new_value
+                    was_updated = True
+        assert was_updated
 
     def _solve_cvxpy_opt_problem(self):
         """
@@ -229,14 +256,19 @@ class BaseConvexOptimizer(BaseOptimizer):
         :raises exceptions.OptimizationError: if problem is not solvable by cvxpy
         """
         try:
-            self._opt = cp.Problem(cp.Minimize(self._objective), self._constraints)
-
-            if self._solver is not None:
-                self._opt.solve(
-                    solver=self._solver, verbose=self._verbose, **self._solver_options
-                )
+            if self._opt is None:
+                self._opt = cp.Problem(cp.Minimize(self._objective), self._constraints)
+                self._initial_objective = self._objective.id
+                self._initial_constraint_ids = {const.id for const in self._constraints}
             else:
-                self._opt.solve(verbose=self._verbose, **self._solver_options)
+                assert self._objective.id == self._initial_objective, \
+                    "The objective function was changed after the initial optimization. " \
+                    "Please create a new instance instead."
+                assert {const.id for const in self._constraints} == self._initial_constraint_ids, \
+                    "The constraints were changed after the initial optimization. " \
+                    "Please create a new instance instead."
+            self._opt.solve(solver=self._solver, verbose=self._verbose, **self._solver_options)
+
         except (TypeError, cp.DCPError) as e:
             raise exceptions.OptimizationError from e
 
@@ -262,6 +294,9 @@ class BaseConvexOptimizer(BaseOptimizer):
         :param new_objective: the objective to be added
         :type new_objective: cp.Expression (i.e function of cp.Variable)
         """
+        if self._opt is not None:
+            raise Exception('Adding objectives to an already solved problem might have unintended consequences.'
+                            'A new instance should be created for the new set of objectives.')
         self._additional_objectives.append(new_objective(self._w, **kwargs))
 
     def add_constraint(self, new_constraint):
@@ -280,6 +315,9 @@ class BaseConvexOptimizer(BaseOptimizer):
         """
         if not callable(new_constraint):
             raise TypeError("New constraint must be provided as a lambda function")
+        if self._opt is not None:
+            raise Exception('Adding constraints to an already solved problem might have unintended consequences.'
+                            'A new instance should be created for the new set of constraints.')
         self._constraints.append(new_constraint(self._w))
 
     def add_sector_constraints(self, sector_mapper, sector_lower, sector_upper):
@@ -316,10 +354,10 @@ class BaseConvexOptimizer(BaseOptimizer):
             )
         for sector in sector_upper:
             is_sector = [sector_mapper[t] == sector for t in self.tickers]
-            self._constraints.append(cp.sum(self._w[is_sector]) <= sector_upper[sector])
+            self.add_constraint(lambda w: cp.sum(w[is_sector]) <= sector_upper[sector])
         for sector in sector_lower:
             is_sector = [sector_mapper[t] == sector for t in self.tickers]
-            self._constraints.append(cp.sum(self._w[is_sector]) >= sector_lower[sector])
+            self.add_constraint(lambda w: cp.sum(w[is_sector]) >= sector_lower[sector])
 
     def convex_objective(self, custom_objective, weights_sum_to_one=True, **kwargs):
         """
@@ -349,7 +387,7 @@ class BaseConvexOptimizer(BaseOptimizer):
             self._objective += obj
 
         if weights_sum_to_one:
-            self._constraints.append(cp.sum(self._w) == 1)
+            self.add_constraint(lambda w: cp.sum(w) == 1)
 
         return self._solve_cvxpy_opt_problem()
 
@@ -413,7 +451,7 @@ class BaseConvexOptimizer(BaseOptimizer):
         # Construct constraints
         final_constraints = []
         if weights_sum_to_one:
-            final_constraints.append({"type": "eq", "fun": lambda x: np.sum(x) - 1})
+            final_constraints.append({"type": "eq", "fun": lambda w: np.sum(w) - 1})
         if constraints is not None:
             final_constraints += constraints
 
@@ -493,3 +531,18 @@ def portfolio_performance(
         if verbose:
             print("Annual volatility: {:.1f}%".format(100 * sigma))
         return None, sigma, None
+
+
+def get_all_args(expression: cp.Expression) -> List[cp.Expression]:
+    if expression.args == []:
+        return [expression]
+    else:
+        return list(flatten([get_all_args(arg) for arg in expression.args]))
+
+
+def flatten(l: Iterable) -> Iterable:
+    for el in l:
+        if isinstance(el, Iterable) and not isinstance(el, (str, bytes)):
+            yield from flatten(el)
+        else:
+            yield el
\ No newline at end of file
diff --git a/pypfopt/efficient_frontier/efficient_cdar.py b/pypfopt/efficient_frontier/efficient_cdar.py
index c9d7f1b..b0e9cac 100644
--- a/pypfopt/efficient_frontier/efficient_cdar.py
+++ b/pypfopt/efficient_frontier/efficient_cdar.py
@@ -125,14 +125,7 @@ class EfficientCDaR(EfficientFrontier):
         for obj in self._additional_objectives:
             self._objective += obj
 
-        self._constraints += [
-            self._z >= self._u[1:] - self._alpha,
-            self._u[1:] >= self._u[:-1] - self.returns.values @ self._w,
-            self._u[0] == 0,
-            self._z >= 0,
-            self._u[1:] >= 0,
-        ]
-
+        self._add_cdar_constraints()
         self._make_weight_sum_constraint(market_neutral)
         return self._solve_cvxpy_opt_problem()
 
@@ -150,26 +143,17 @@ class EfficientCDaR(EfficientFrontier):
         :return: asset weights for the optimal portfolio
         :rtype: OrderedDict
         """
-        self._objective = self._alpha + 1.0 / (
-            len(self.returns) * (1 - self._beta)
-        ) * cp.sum(self._z)
 
-        for obj in self._additional_objectives:
-            self._objective += obj
-
-        self._constraints += [
-            self._z >= self._u[1:] - self._alpha,
-            self._u[1:] >= self._u[:-1] - self.returns.values @ self._w,
-            self._u[0] == 0,
-            self._z >= 0,
-            self._u[1:] >= 0,
-        ]
-
-        ret = self.expected_returns.T @ self._w
-        self._constraints.append(ret >= target_return)
-
-        self._make_weight_sum_constraint(market_neutral)
-        return self._solve_cvxpy_opt_problem()
+        update_existing_parameter = self.is_parameter_defined('target_return')
+        if update_existing_parameter:
+            self._validate_market_neutral(market_neutral)
+            self.update_parameter_value('target_return', target_return)
+            return self._solve_cvxpy_opt_problem()
+        else:
+            ret = self.expected_returns.T @ self._w
+            target_return_par = cp.Parameter(value=target_return, name='target_return', nonneg=True)
+            self._constraints.append(ret >= target_return_par)
+            return self.min_cdar(market_neutral)
 
     def efficient_risk(self, target_cdar, market_neutral=False):
         """
@@ -185,27 +169,36 @@ class EfficientCDaR(EfficientFrontier):
         :return: asset weights for the efficient risk portfolio
         :rtype: OrderedDict
         """
-        self._objective = objective_functions.portfolio_return(
-            self._w, self.expected_returns
-        )
-        for obj in self._additional_objectives:
-            self._objective += obj
 
-        cdar = self._alpha + 1.0 / (len(self.returns) * (1 - self._beta)) * cp.sum(
-            self._z
-        )
-        self._constraints += [
-            cdar <= target_cdar,
-            self._z >= self._u[1:] - self._alpha,
-            self._u[1:] >= self._u[:-1] - self.returns.values @ self._w,
-            self._u[0] == 0,
-            self._z >= 0,
-            self._u[1:] >= 0,
-        ]
+        update_existing_parameter = self.is_parameter_defined('target_cdar')
+        if update_existing_parameter:
+            self._validate_market_neutral(market_neutral)
+            self.update_parameter_value('target_cdar', target_cdar)
+        else:
+            self._objective = objective_functions.portfolio_return(
+                self._w, self.expected_returns
+            )
+            for obj in self._additional_objectives:
+                self._objective += obj
 
-        self._make_weight_sum_constraint(market_neutral)
+            cdar = self._alpha + 1.0 / (len(self.returns) * (1 - self._beta)) * cp.sum(
+                self._z
+            )
+            target_cdar_par = cp.Parameter(value=target_cdar, name='target_cdar', nonneg=True)
+            self.add_constraint(lambda _: cdar <= target_cdar_par)
+
+            self._add_cdar_constraints()
+
+            self._make_weight_sum_constraint(market_neutral)
         return self._solve_cvxpy_opt_problem()
 
+    def _add_cdar_constraints(self) -> None:
+        self.add_constraint(lambda _: self._z >= self._u[1:] - self._alpha)
+        self.add_constraint(lambda w: self._u[1:] >= self._u[:-1] - self.returns.values @ w)
+        self.add_constraint(lambda _: self._u[0] == 0)
+        self.add_constraint(lambda _: self._z >= 0)
+        self.add_constraint(lambda _: self._u[1:] >= 0)
+
     def portfolio_performance(self, verbose=False):
         """
         After optimising, calculate (and optionally print) the performance of the optimal
diff --git a/pypfopt/efficient_frontier/efficient_cvar.py b/pypfopt/efficient_frontier/efficient_cvar.py
index fd72750..cce6a40 100644
--- a/pypfopt/efficient_frontier/efficient_cvar.py
+++ b/pypfopt/efficient_frontier/efficient_cvar.py
@@ -57,7 +57,7 @@ class EfficientCVaR(EfficientFrontier):
     ):
         """
         :param expected_returns: expected returns for each asset. Can be None if
-                                optimising for semideviation only.
+                                optimising for conditional value at risk only.
         :type expected_returns: pd.Series, list, np.ndarray
         :param returns: (historic) returns for all your assets (no NaNs).
                                  See ``expected_returns.returns_from_prices``.
@@ -126,10 +126,8 @@ class EfficientCVaR(EfficientFrontier):
         for obj in self._additional_objectives:
             self._objective += obj
 
-        self._constraints += [
-            self._u >= 0.0,
-            self.returns.values @ self._w + self._alpha + self._u >= 0.0,
-        ]
+        self.add_constraint(lambda _: self._u >= 0.0)
+        self.add_constraint(lambda w: self.returns.values @ w + self._alpha + self._u >= 0.0)
 
         self._make_weight_sum_constraint(market_neutral)
         return self._solve_cvxpy_opt_problem()
@@ -148,22 +146,26 @@ class EfficientCVaR(EfficientFrontier):
         :return: asset weights for the optimal portfolio
         :rtype: OrderedDict
         """
-        self._objective = self._alpha + 1.0 / (
-            len(self.returns) * (1 - self._beta)
-        ) * cp.sum(self._u)
+        update_existing_parameter = self.is_parameter_defined('target_return')
+        if update_existing_parameter:
+            self._validate_market_neutral(market_neutral)
+            self.update_parameter_value('target_return', target_return)
+        else:
+            self._objective = self._alpha + 1.0 / (
+                len(self.returns) * (1 - self._beta)
+            ) * cp.sum(self._u)
 
-        for obj in self._additional_objectives:
-            self._objective += obj
+            for obj in self._additional_objectives:
+                self._objective += obj
 
-        self._constraints += [
-            self._u >= 0.0,
-            self.returns.values @ self._w + self._alpha + self._u >= 0.0,
-        ]
+            self.add_constraint(lambda _: self._u >= 0.0)
+            self.add_constraint(lambda w: self.returns.values @ w + self._alpha + self._u >= 0.0)
 
-        ret = self.expected_returns.T @ self._w
-        self._constraints.append(ret >= target_return)
+            ret = self.expected_returns.T @ self._w
+            target_return_par = cp.Parameter(name='target_return', value=target_return)
+            self.add_constraint(lambda _: ret >= target_return_par)
 
-        self._make_weight_sum_constraint(market_neutral)
+            self._make_weight_sum_constraint(market_neutral)
         return self._solve_cvxpy_opt_problem()
 
     def efficient_risk(self, target_cvar, market_neutral=False):
@@ -172,7 +174,7 @@ class EfficientCVaR(EfficientFrontier):
         The resulting portfolio will have a CVaR less than the target
         (but not guaranteed to be equal).
 
-        :param target_cvar: the desired maximum semideviation of the resulting portfolio.
+        :param target_cvar: the desired conditional value at risk of the resulting portfolio.
         :type target_cvar: float
         :param market_neutral: whether the portfolio should be market neutral (weights sum to zero),
                                defaults to False. Requires negative lower weight bound.
@@ -180,22 +182,27 @@ class EfficientCVaR(EfficientFrontier):
         :return: asset weights for the efficient risk portfolio
         :rtype: OrderedDict
         """
-        self._objective = objective_functions.portfolio_return(
-            self._w, self.expected_returns
-        )
-        for obj in self._additional_objectives:
-            self._objective += obj
+        update_existing_parameter = self.is_parameter_defined('target_cvar')
+        if update_existing_parameter:
+            self._validate_market_neutral(market_neutral)
+            self.update_parameter_value('target_cvar', target_cvar)
+        else:
+            self._objective = objective_functions.portfolio_return(
+                self._w, self.expected_returns
+            )
+            for obj in self._additional_objectives:
+                self._objective += obj
 
-        cvar = self._alpha + 1.0 / (len(self.returns) * (1 - self._beta)) * cp.sum(
-            self._u
-        )
-        self._constraints += [
-            cvar <= target_cvar,
-            self._u >= 0.0,
-            self.returns.values @ self._w + self._alpha + self._u >= 0.0,
-        ]
+            cvar = self._alpha + 1.0 / (len(self.returns) * (1 - self._beta)) * cp.sum(
+                self._u
+            )
+            target_cvar_par = cp.Parameter(value=target_cvar, name='target_cvar', nonneg=True)
 
-        self._make_weight_sum_constraint(market_neutral)
+            self.add_constraint(lambda _: cvar <= target_cvar_par)
+            self.add_constraint(lambda _: self._u >= 0.0)
+            self.add_constraint(lambda w: self.returns.values @ w + self._alpha + self._u >= 0.0)
+
+            self._make_weight_sum_constraint(market_neutral)
         return self._solve_cvxpy_opt_problem()
 
     def portfolio_performance(self, verbose=False):
diff --git a/pypfopt/efficient_frontier/efficient_frontier.py b/pypfopt/efficient_frontier/efficient_frontier.py
index c1cf460..cd7eceb 100644
--- a/pypfopt/efficient_frontier/efficient_frontier.py
+++ b/pypfopt/efficient_frontier/efficient_frontier.py
@@ -2,8 +2,9 @@
 The ``efficient_frontier`` submodule houses the EfficientFrontier class, which generates
 classical mean-variance optimal portfolios for a variety of objectives and constraints
 """
-
+import copy
 import warnings
+
 import numpy as np
 import pandas as pd
 import cvxpy as cp
@@ -86,6 +87,8 @@ class EfficientFrontier(base_optimizer.BaseConvexOptimizer):
         self.expected_returns = EfficientFrontier._validate_expected_returns(
             expected_returns
         )
+        self._max_return_value = None
+        self._market_neutral = None
 
         if self.expected_returns is None:
             num_assets = len(cov_matrix)
@@ -178,9 +181,11 @@ class EfficientFrontier(base_optimizer.BaseConvexOptimizer):
                 del self._constraints[0]
                 del self._constraints[0]
 
-            self._constraints.append(cp.sum(self._w) == 0)
+
+            self.add_constraint(lambda w: cp.sum(w) == 0)
         else:
-            self._constraints.append(cp.sum(self._w) == 1)
+            self.add_constraint(lambda w: cp.sum(w) == 1)
+        self._market_neutral = is_market_neutral
 
     def min_volatility(self):
         """
@@ -195,7 +200,7 @@ class EfficientFrontier(base_optimizer.BaseConvexOptimizer):
         for obj in self._additional_objectives:
             self._objective += obj
 
-        self._constraints.append(cp.sum(self._w) == 1)
+        self.add_constraint(lambda w: cp.sum(w) == 1)
         return self._solve_cvxpy_opt_problem()
 
     def _max_return(self, return_value=True):
@@ -212,13 +217,10 @@ class EfficientFrontier(base_optimizer.BaseConvexOptimizer):
             self._w, self.expected_returns
         )
 
-        self._constraints.append(cp.sum(self._w) == 1)
+        self.add_constraint(lambda w: cp.sum(w) == 1)
 
         res = self._solve_cvxpy_opt_problem()
 
-        #  Cleanup constraints since this is a helper method
-        del self._constraints[-1]
-
         if return_value:
             return -self._opt.value
         else:
@@ -308,13 +310,18 @@ class EfficientFrontier(base_optimizer.BaseConvexOptimizer):
         if risk_aversion <= 0:
             raise ValueError("risk aversion coefficient must be greater than zero")
 
-        self._objective = objective_functions.quadratic_utility(
-            self._w, self.expected_returns, self.cov_matrix, risk_aversion=risk_aversion
-        )
-        for obj in self._additional_objectives:
-            self._objective += obj
+        update_existing_parameter = self.is_parameter_defined('risk_aversion')
+        if update_existing_parameter:
+            self._validate_market_neutral(market_neutral)
+            self.update_parameter_value('risk_aversion', risk_aversion)
+        else:
+            self._objective = objective_functions.quadratic_utility(
+                self._w, self.expected_returns, self.cov_matrix, risk_aversion=risk_aversion
+            )
+            for obj in self._additional_objectives:
+                self._objective += obj
 
-        self._make_weight_sum_constraint(market_neutral)
+            self._make_weight_sum_constraint(market_neutral)
         return self._solve_cvxpy_opt_problem()
 
     def efficient_risk(self, target_volatility, market_neutral=False):
@@ -345,16 +352,22 @@ class EfficientFrontier(base_optimizer.BaseConvexOptimizer):
                 )
             )
 
-        self._objective = objective_functions.portfolio_return(
-            self._w, self.expected_returns
-        )
-        variance = objective_functions.portfolio_variance(self._w, self.cov_matrix)
+        update_existing_parameter = self.is_parameter_defined('target_variance')
+        if update_existing_parameter:
+            self._validate_market_neutral(market_neutral)
+            self.update_parameter_value('target_variance', target_volatility ** 2)
+        else:
+            self._objective = objective_functions.portfolio_return(
+                self._w, self.expected_returns
+            )
+            variance = objective_functions.portfolio_variance(self._w, self.cov_matrix)
 
-        for obj in self._additional_objectives:
-            self._objective += obj
+            for obj in self._additional_objectives:
+                self._objective += obj
 
-        self._constraints.append(variance <= target_volatility ** 2)
-        self._make_weight_sum_constraint(market_neutral)
+            target_variance = cp.Parameter(name="target_variance", value=target_volatility ** 2, nonneg=True)
+            self.add_constraint(lambda _: variance <= target_variance)
+            self._make_weight_sum_constraint(market_neutral)
         return self._solve_cvxpy_opt_problem()
 
     def efficient_return(self, target_return, market_neutral=False):
@@ -373,23 +386,31 @@ class EfficientFrontier(base_optimizer.BaseConvexOptimizer):
         """
         if not isinstance(target_return, float) or target_return < 0:
             raise ValueError("target_return should be a positive float")
-        if target_return > self._max_return():
+        if not self._max_return_value:
+            self._max_return_value = copy.deepcopy(self)._max_return()
+        if target_return > self._max_return_value:
             raise ValueError(
                 "target_return must be lower than the maximum possible return"
             )
 
-        self._objective = objective_functions.portfolio_variance(
-            self._w, self.cov_matrix
-        )
-        ret = objective_functions.portfolio_return(
-            self._w, self.expected_returns, negative=False
-        )
+        update_existing_parameter = self.is_parameter_defined('target_return')
+        if update_existing_parameter:
+            self._validate_market_neutral(market_neutral)
+            self.update_parameter_value('target_return', target_return)
+        else:
+            self._objective = objective_functions.portfolio_variance(
+                self._w, self.cov_matrix
+            )
+            ret = objective_functions.portfolio_return(
+                self._w, self.expected_returns, negative=False
+            )
 
-        for obj in self._additional_objectives:
-            self._objective += obj
+            for obj in self._additional_objectives:
+                self._objective += obj
 
-        self._constraints.append(ret >= target_return)
-        self._make_weight_sum_constraint(market_neutral)
+            target_return_par = cp.Parameter(name='target_return', value=target_return)
+            self.add_constraint(lambda _: ret >= target_return_par)
+            self._make_weight_sum_constraint(market_neutral)
         return self._solve_cvxpy_opt_problem()
 
     def portfolio_performance(self, verbose=False, risk_free_rate=0.02):
@@ -423,3 +444,6 @@ class EfficientFrontier(base_optimizer.BaseConvexOptimizer):
             verbose,
             risk_free_rate,
         )
+
+    def _validate_market_neutral(self, market_neutral: bool) -> None:
+        assert self._market_neutral == market_neutral, 'A new instance must be created when changing market_neutral'
diff --git a/pypfopt/efficient_frontier/efficient_semivariance.py b/pypfopt/efficient_frontier/efficient_semivariance.py
index cb83c73..6ac063e 100644
--- a/pypfopt/efficient_frontier/efficient_semivariance.py
+++ b/pypfopt/efficient_frontier/efficient_semivariance.py
@@ -125,7 +125,7 @@ class EfficientSemivariance(EfficientFrontier):
             self._objective += obj
 
         B = (self.returns.values - self.benchmark) / np.sqrt(self._T)
-        self._constraints.append(B @ self._w - p + n == 0)
+        self.add_constraint(lambda w: B @ w - p + n == 0)
         self._make_weight_sum_constraint(market_neutral)
         return self._solve_cvxpy_opt_problem()
 
@@ -146,17 +146,23 @@ class EfficientSemivariance(EfficientFrontier):
         if risk_aversion <= 0:
             raise ValueError("risk aversion coefficient must be greater than zero")
 
-        p = cp.Variable(self._T, nonneg=True)
-        n = cp.Variable(self._T, nonneg=True)
-        mu = objective_functions.portfolio_return(self._w, self.expected_returns)
-        mu /= self.frequency
-        self._objective = mu + 0.5 * risk_aversion * cp.sum(cp.square(n))
-        for obj in self._additional_objectives:
-            self._objective += obj
+        update_existing_parameter = self.is_parameter_defined('risk_aversion')
+        if update_existing_parameter:
+            self._validate_market_neutral(market_neutral)
+            self.update_parameter_value('risk_aversion', risk_aversion)
+        else:
+            p = cp.Variable(self._T, nonneg=True)
+            n = cp.Variable(self._T, nonneg=True)
+            mu = objective_functions.portfolio_return(self._w, self.expected_returns)
+            mu /= self.frequency
+            risk_aversion_par = cp.Parameter(value=risk_aversion, name='risk_aversion', nonneg=True)
+            self._objective = mu + 0.5 * risk_aversion_par * cp.sum(cp.square(n))
+            for obj in self._additional_objectives:
+                self._objective += obj
 
-        B = (self.returns.values - self.benchmark) / np.sqrt(self._T)
-        self._constraints.append(B @ self._w - p + n == 0)
-        self._make_weight_sum_constraint(market_neutral)
+            B = (self.returns.values - self.benchmark) / np.sqrt(self._T)
+            self.add_constraint(lambda w: B @ w - p + n == 0)
+            self._make_weight_sum_constraint(market_neutral)
         return self._solve_cvxpy_opt_problem()
 
     def efficient_risk(self, target_semideviation, market_neutral=False):
@@ -173,21 +179,25 @@ class EfficientSemivariance(EfficientFrontier):
         :return: asset weights for the efficient risk portfolio
         :rtype: OrderedDict
         """
-        self._objective = objective_functions.portfolio_return(
-            self._w, self.expected_returns
-        )
-        for obj in self._additional_objectives:
-            self._objective += obj
+        update_existing_parameter = self.is_parameter_defined('target_semivariance')
+        if update_existing_parameter:
+            self._validate_market_neutral(market_neutral)
+            self.update_parameter_value('target_semivariance', target_semideviation ** 2)
+        else:
+            self._objective = objective_functions.portfolio_return(
+                self._w, self.expected_returns
+            )
+            for obj in self._additional_objectives:
+                self._objective += obj
 
-        p = cp.Variable(self._T, nonneg=True)
-        n = cp.Variable(self._T, nonneg=True)
+            p = cp.Variable(self._T, nonneg=True)
+            n = cp.Variable(self._T, nonneg=True)
 
-        self._constraints.append(
-            self.frequency * cp.sum(cp.square(n)) <= (target_semideviation ** 2)
-        )
-        B = (self.returns.values - self.benchmark) / np.sqrt(self._T)
-        self._constraints.append(B @ self._w - p + n == 0)
-        self._make_weight_sum_constraint(market_neutral)
+            target_semivariance = cp.Parameter(value=target_semideviation**2, name='target_semivariance', nonneg=True)
+            self.add_constraint(lambda _: self.frequency * cp.sum(cp.square(n)) <= target_semivariance)
+            B = (self.returns.values - self.benchmark) / np.sqrt(self._T)
+            self.add_constraint(lambda w: B @ w - p + n == 0)
+            self._make_weight_sum_constraint(market_neutral)
         return self._solve_cvxpy_opt_problem()
 
     def efficient_return(self, target_return, market_neutral=False):
@@ -210,18 +220,24 @@ class EfficientSemivariance(EfficientFrontier):
             raise ValueError(
                 "target_return must be lower than the largest expected return"
             )
-        p = cp.Variable(self._T, nonneg=True)
-        n = cp.Variable(self._T, nonneg=True)
-        self._objective = cp.sum(cp.square(n))
-        for obj in self._additional_objectives:
-            self._objective += obj
 
-        self._constraints.append(
-            cp.sum(self._w @ self.expected_returns) >= target_return
-        )
-        B = (self.returns.values - self.benchmark) / np.sqrt(self._T)
-        self._constraints.append(B @ self._w - p + n == 0)
-        self._make_weight_sum_constraint(market_neutral)
+        update_existing_parameter = self.is_parameter_defined('target_return')
+        if update_existing_parameter:
+            self._validate_market_neutral(market_neutral)
+            self.update_parameter_value('target_return', target_return)
+        else:
+
+            p = cp.Variable(self._T, nonneg=True)
+            n = cp.Variable(self._T, nonneg=True)
+            self._objective = cp.sum(cp.square(n))
+            for obj in self._additional_objectives:
+                self._objective += obj
+
+            target_return_par = cp.Parameter(name='target_return', value=target_return)
+            self.add_constraint(lambda w: cp.sum(w @ self.expected_returns) >= target_return_par)
+            B = (self.returns.values - self.benchmark) / np.sqrt(self._T)
+            self.add_constraint(lambda w: B @ w - p + n == 0)
+            self._make_weight_sum_constraint(market_neutral)
         return self._solve_cvxpy_opt_problem()
 
     def portfolio_performance(self, verbose=False, risk_free_rate=0.02):
diff --git a/pypfopt/objective_functions.py b/pypfopt/objective_functions.py
index 3f6e86d..fb9acb2 100644
--- a/pypfopt/objective_functions.py
+++ b/pypfopt/objective_functions.py
@@ -158,7 +158,8 @@ def quadratic_utility(w, expected_returns, cov_matrix, risk_aversion, negative=T
     mu = w @ expected_returns
     variance = cp.quad_form(w, cov_matrix)
 
-    utility = mu - 0.5 * risk_aversion * variance
+    risk_aversion_par = cp.Parameter(value=risk_aversion, name='risk_aversion', nonneg=True)
+    utility = mu - 0.5 * risk_aversion_par * variance
     return _objective_value(w, sign * utility)
 
 
diff --git a/pypfopt/plotting.py b/pypfopt/plotting.py
index 60f56eb..c9bbc38 100644
--- a/pypfopt/plotting.py
+++ b/pypfopt/plotting.py
@@ -168,15 +168,13 @@ def _plot_ef(ef, ef_param, ef_param_range, ax, show_assets):
 
     # Create a portfolio for each value of ef_param_range
     for param_value in ef_param_range:
-        ef_i = copy.deepcopy(ef)
-
         try:
             if ef_param == "utility":
-                ef_i.max_quadratic_utility(param_value)
+                ef.max_quadratic_utility(param_value)
             elif ef_param == "risk":
-                ef_i.efficient_risk(param_value)
+                ef.efficient_risk(param_value)
             elif ef_param == "return":
-                ef_i.efficient_return(param_value)
+                ef.efficient_return(param_value)
             else:
                 raise NotImplementedError(
                     "ef_param should be one of {'utility', 'risk', 'return'}"
@@ -190,7 +188,7 @@ def _plot_ef(ef, ef_param, ef_param_range, ax, show_assets):
                 )
             )
 
-        ret, sigma, _ = ef_i.portfolio_performance()
+        ret, sigma, _ = ef.portfolio_performance()
         mus.append(ret)
         sigmas.append(sigma)
 
diff --git a/tests/test_base_optimizer.py b/tests/test_base_optimizer.py
index 2179e86..e8c8eaa 100644
--- a/tests/test_base_optimizer.py
+++ b/tests/test_base_optimizer.py
@@ -7,7 +7,7 @@ import pandas as pd
 import pytest
 import cvxpy as cp
 
-from pypfopt import EfficientFrontier
+from pypfopt import EfficientFrontier, objective_functions
 from pypfopt import exceptions
 from pypfopt.base_optimizer import portfolio_performance, BaseOptimizer
 from tests.utilities_for_tests import get_data, setup_efficient_frontier
@@ -56,7 +56,10 @@ def test_weight_bounds_minus_one_to_one():
         *setup_efficient_frontier(data_only=True), weight_bounds=(-1, 1)
     )
     assert ef.max_sharpe()
-    assert ef.min_volatility()
+    ef2 = EfficientFrontier(
+        *setup_efficient_frontier(data_only=True), weight_bounds=(-1, 1)
+    )
+    assert ef2.min_volatility()
 
 
 def test_none_bounds():
@@ -189,12 +192,13 @@ def test_efficient_frontier_init_errors():
 
 
 def test_set_weights():
-    ef = setup_efficient_frontier()
-    w1 = ef.min_volatility()
-    test_weights = ef.weights
-    ef.min_volatility()
-    ef.set_weights(w1)
-    np.testing.assert_array_almost_equal(test_weights, ef.weights)
+    ef1 = setup_efficient_frontier()
+    w1 = ef1.min_volatility()
+    test_weights = ef1.weights
+    ef2 = setup_efficient_frontier()
+    ef2.min_volatility()
+    ef2.set_weights(w1)
+    np.testing.assert_array_almost_equal(test_weights, ef2.weights)
 
 
 def test_save_weights_to_file():
@@ -287,3 +291,38 @@ def test_problem_access():
     ef = setup_efficient_frontier()
     ef.max_sharpe()
     assert isinstance(ef._opt, cp.Problem)
+
+
+def test_exception_immutability():
+    ef = setup_efficient_frontier()
+    ef.efficient_return(0.2)
+
+    with pytest.raises(Exception,
+                       match='Adding constraints to an already solved problem might have unintended consequences'):
+        ef.min_volatility()
+
+    ef = setup_efficient_frontier()
+    ef.efficient_return(0.2)
+    with pytest.raises(Exception,
+                       match='Adding constraints to an already solved problem might have unintended consequences'):
+        ef.add_constraint(lambda w: w >= 0.1)
+
+    ef = setup_efficient_frontier()
+    ef.efficient_return(0.2)
+    prev_w = np.array([1 / ef.n_assets] * ef.n_assets)
+    with pytest.raises(Exception,
+                       match='Adding objectives to an already solved problem might have unintended consequences'):
+        ef.add_objective(objective_functions.transaction_cost, w_prev=prev_w)
+
+    ef = setup_efficient_frontier()
+    ef.efficient_return(0.2)
+    ef._constraints += [ef._w >= 0.1]
+    with pytest.raises(Exception,
+                       match='The constraints were changed after the initial optimization'):
+        ef.efficient_return(0.2)
+
+    ef = setup_efficient_frontier()
+    ef.efficient_return(0.2, market_neutral=True)
+    with pytest.raises(Exception,
+                       match='A new instance must be created when changing market_neutral'):
+        ef.efficient_return(0.2, market_neutral=False)
diff --git a/tests/test_custom_objectives.py b/tests/test_custom_objectives.py
index 97cbad1..26645fe 100644
--- a/tests/test_custom_objectives.py
+++ b/tests/test_custom_objectives.py
@@ -40,7 +40,6 @@ def test_custom_convex_abs_exposure():
     )
 
     ef.add_constraint(lambda x: cp.norm(x, 1) <= 2)
-    ef.min_volatility()
     ef.convex_objective(
         objective_functions.portfolio_variance,
         cov_matrix=ef.cov_matrix,
diff --git a/tests/test_efficient_cdar.py b/tests/test_efficient_cdar.py
index 5a17962..a4ddf76 100644
--- a/tests/test_efficient_cdar.py
+++ b/tests/test_efficient_cdar.py
@@ -350,15 +350,14 @@ def test_efficient_risk_L2_reg():
         atol=1e-4,
     )
 
-    ef2 = setup_efficient_cdar()
-    cd.add_objective(objective_functions.L2_reg, gamma=1)
-    ef2.efficient_risk(0.18)
+    cd2 = setup_efficient_cdar()
+    cd2.efficient_risk(0.18)
 
     # L2_reg should pull close to equal weight
     equal_weight = np.full((cd.n_assets,), 1 / cd.n_assets)
     assert (
             np.abs(equal_weight - cd.weights).sum()
-            < np.abs(equal_weight - ef2.weights).sum()
+            < np.abs(equal_weight - cd2.weights).sum()
     )
 
 
@@ -455,3 +454,13 @@ def test_cdar_errors():
     historical_rets = historical_rets.iloc[:, :-1]
     with pytest.raises(ValueError):
         EfficientCDaR(mu, historical_rets)
+
+
+def test_parametrization():
+    cd = setup_efficient_cdar()
+    cd.efficient_risk(0.08)
+    cd.efficient_risk(0.07)
+
+    cd = setup_efficient_cdar()
+    cd.efficient_return(0.08)
+    cd.efficient_return(0.07)
diff --git a/tests/test_efficient_cvar.py b/tests/test_efficient_cvar.py
index 5645955..30e3829 100644
--- a/tests/test_efficient_cvar.py
+++ b/tests/test_efficient_cvar.py
@@ -348,15 +348,14 @@ def test_efficient_risk_L2_reg():
         atol=1e-4,
     )
 
-    ef2 = setup_efficient_cvar()
-    cv.add_objective(objective_functions.L2_reg, gamma=1)
-    ef2.efficient_risk(0.03)
+    cv2 = setup_efficient_cvar()
+    cv2.efficient_risk(0.19)
 
     # L2_reg should pull close to equal weight
     equal_weight = np.full((cv.n_assets,), 1 / cv.n_assets)
     assert (
         np.abs(equal_weight - cv.weights).sum()
-        < np.abs(equal_weight - ef2.weights).sum()
+        < np.abs(equal_weight - cv2.weights).sum()
     )
 
 
@@ -456,3 +455,13 @@ def test_cvar_errors():
     historical_rets = historical_rets.iloc[:, :-1]
     with pytest.raises(ValueError):
         EfficientCVaR(mu, historical_rets)
+
+
+def test_parametrization():
+    cv = setup_efficient_cvar()
+    cv.efficient_risk(0.19)
+    cv.efficient_risk(0.19)
+
+    cv = setup_efficient_cvar()
+    cv.efficient_return(0.25)
+    cv.efficient_return(0.25)
diff --git a/tests/test_efficient_frontier.py b/tests/test_efficient_frontier.py
index 107af3f..0ded8ab 100644
--- a/tests/test_efficient_frontier.py
+++ b/tests/test_efficient_frontier.py
@@ -203,8 +203,9 @@ def test_min_volatility_L2_reg_many_values():
     ef.min_volatility()
     # Count the number of weights more 1%
     initial_number = sum(ef.weights > 0.01)
-    for _ in range(10):
-        ef.add_objective(objective_functions.L2_reg, gamma=0.05)
+    for gamma_multiplier in range(1, 10):
+        ef = setup_efficient_frontier()
+        ef.add_objective(objective_functions.L2_reg, gamma=0.05*gamma_multiplier)
         ef.min_volatility()
         np.testing.assert_almost_equal(ef.weights.sum(), 1)
         new_number = sum(ef.weights > 0.01)
@@ -388,7 +389,7 @@ def test_max_sharpe_error():
 
     # An unsupported constraint type, which is incidentally meaningless.
     v = cp.Variable((2, 2), PSD=True)
-    ef._constraints.append(v >> np.zeros((2, 2)))
+    ef.add_constraint(lambda _: v >> np.zeros((2, 2)))
     with pytest.raises(TypeError):
         ef.max_sharpe()
 
@@ -824,8 +825,9 @@ def test_max_quadratic_utility():
 
     ret1, var1, _ = ef.portfolio_performance()
     # increasing risk_aversion should lower both vol and return
-    ef.max_quadratic_utility(10)
-    ret2, var2, _ = ef.portfolio_performance()
+    ef2 = setup_efficient_frontier()
+    ef2.max_quadratic_utility(10)
+    ret2, var2, _ = ef2.portfolio_performance()
     assert ret2 < ret1 and var2 < var1
 
 
diff --git a/tests/test_efficient_semivariance.py b/tests/test_efficient_semivariance.py
index 9110dc9..9a96989 100644
--- a/tests/test_efficient_semivariance.py
+++ b/tests/test_efficient_semivariance.py
@@ -465,7 +465,6 @@ def test_efficient_risk_L2_reg():
     )
 
     ef2 = setup_efficient_semivariance()
-    es.add_objective(objective_functions.L2_reg, gamma=1)
     ef2.efficient_risk(0.19)
 
     # L2_reg should pull close to equal weight
@@ -582,3 +581,17 @@ def test_efficient_semivariance_vs_heuristic_weekly():
 
     assert semi_deviation < semi_deviation_ef
     assert mu_es / semi_deviation > mu_ef / semi_deviation_ef
+
+
+def test_parametrization():
+    es = setup_efficient_semivariance()
+    es.efficient_risk(0.19)
+    es.efficient_risk(0.19)
+
+    es = setup_efficient_semivariance()
+    es.efficient_return(0.25)
+    es.efficient_return(0.25)
+
+    es = setup_efficient_semivariance()
+    es.max_quadratic_utility(1)
+    es.max_quadratic_utility(1)