Repository: pybamm-team/PyBaMM

Issue:
Adding constraints to algebraic solver
It seems the algebraic solver sometimes fails when a state goes outside an obvious bound (e.g. concentration going negative) which makes the algebraic equations return nan. Opening this issue to brainstorm some solutions, so suggestions welcome:

1. Simply replace `nan` with a very large number. Will try this first and hope it works as this would be the simplest solution. The other two solutions require manually specifying some bounds in the models (e.g. using the existing event functionality)
2. Impose hard constraints, but neither `scipy.optimize.root` nor casadi have a very clean way of providing constraints
3. Impose soft constraints: add a function of `c` with singularities as c-> 0 and c -> 1 to the function being optimized

N.B. the use case here is for fuel cell models, solving directly for the steady state. 

Discussion:

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 4ffe632a8..a0c480ad9 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,6 +4,7 @@
 
 -   Reformatted Landesfeind electrolytes ([#1064](https://github.com/pybamm-team/PyBaMM/pull/1064))
 -   Adapted examples to be run in Google Colab ([#1061](https://github.com/pybamm-team/PyBaMM/pull/1061))
+-   Added some new solvers for algebraic models ([#1059](https://github.com/pybamm-team/PyBaMM/pull/1059))
 -   Added `length_scales` attribute to models ([#1058](https://github.com/pybamm-team/PyBaMM/pull/1058))
 -   Added averaging in secondary dimensions ([#1057](https://github.com/pybamm-team/PyBaMM/pull/1057))
 
diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index 2afdb9102..1e6de70f9 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -160,8 +160,6 @@ class Discretisation(object):
         # Set the y split for variables
         pybamm.logger.info("Set variable slices for {}".format(model.name))
         self.set_variable_slices(variables)
-        # Keep a record of y_slices in the model
-        model.y_slices = self.y_slices_explicit
 
         # now add extrapolated external variables to the boundary conditions
         # if required by the spatial method
@@ -183,6 +181,11 @@ class Discretisation(object):
             # create an empty copy of the original model
             model_disc = model.new_copy()
 
+        # Keep a record of y_slices in the model
+        model_disc.y_slices = self.y_slices_explicit
+        # Keep a record of the bounds in the model
+        model_disc.bounds = self.bounds
+
         model_disc.bcs = self.bcs
 
         pybamm.logger.info("Discretise initial conditions for {}".format(model.name))
@@ -240,11 +243,13 @@ class Discretisation(object):
         variables : iterable of :class:`pybamm.Variables`
             The variables for which to set slices
         """
-        # Set up y_slices
+        # Set up y_slices and bounds
         y_slices = defaultdict(list)
         y_slices_explicit = defaultdict(list)
         start = 0
         end = 0
+        lower_bounds = []
+        upper_bounds = []
         # Iterate through unpacked variables, adding appropriate slices to y_slices
         for variable in variables:
             # Add up the size of all the domains in variable.domain
@@ -261,13 +266,23 @@ class Discretisation(object):
                     for child, mesh in meshes.items():
                         for domain_mesh in mesh:
                             end += domain_mesh.npts_for_broadcast_to_nodes
+                        # Add to slices
                         y_slices[child.id].append(slice(start, end))
                         y_slices_explicit[child].append(slice(start, end))
+                        # Add to bounds
+                        lower_bounds.extend([child.bounds[0]] * (end - start))
+                        upper_bounds.extend([child.bounds[1]] * (end - start))
+                        # Increment start
                         start = end
             else:
                 end += self._get_variable_size(variable)
+                # Add to slices
                 y_slices[variable.id].append(slice(start, end))
                 y_slices_explicit[variable].append(slice(start, end))
+                # Add to bounds
+                lower_bounds.extend([variable.bounds[0]] * (end - start))
+                upper_bounds.extend([variable.bounds[1]] * (end - start))
+                # Increment start
                 start = end
 
         # Convert y_slices back to normal dictionary
@@ -275,6 +290,9 @@ class Discretisation(object):
         # Also keep a record of what the y_slices are, to be stored in the model
         self.y_slices_explicit = dict(y_slices_explicit)
 
+        # Also keep a record of bounds
+        self.bounds = (np.array(lower_bounds), np.array(upper_bounds))
+
         # reset discretised_symbols
         self._discretised_symbols = {}
 
@@ -885,7 +903,7 @@ class Discretisation(object):
                 # Broadcast new_child to the domain specified by symbol.domain
                 # Different discretisations may broadcast differently
                 if symbol.domain == []:
-                    symbol = disc_child * pybamm.Vector(np.array([1]))
+                    symbol = disc_child * pybamm.Vector([1])
                 else:
                     symbol = spatial_method.broadcast(
                         disc_child,
@@ -921,7 +939,7 @@ class Discretisation(object):
             return pybamm.StateVectorDot(
                 *self.y_slices[symbol.get_variable().id],
                 domain=symbol.domain,
-                auxiliary_domains=symbol.auxiliary_domains
+                auxiliary_domains=symbol.auxiliary_domains,
             )
 
         elif isinstance(symbol, pybamm.Variable):
@@ -972,7 +990,7 @@ class Discretisation(object):
                 return pybamm.StateVector(
                     *y_slices,
                     domain=symbol.domain,
-                    auxiliary_domains=symbol.auxiliary_domains
+                    auxiliary_domains=symbol.auxiliary_domains,
                 )
 
         elif isinstance(symbol, pybamm.SpatialVariable):
diff --git a/pybamm/expression_tree/array.py b/pybamm/expression_tree/array.py
index 41bed0708..95622571b 100644
--- a/pybamm/expression_tree/array.py
+++ b/pybamm/expression_tree/array.py
@@ -13,8 +13,9 @@ class Array(pybamm.Symbol):
     Parameters
     ----------
 
-    entries : numpy.array
-        the array associated with the node
+    entries : numpy.array or list
+        the array associated with the node. If a list is provided, it is converted to a
+        numpy array
     name : str, optional
         the name of the node
     domain : iterable of str, optional
@@ -35,6 +36,8 @@ class Array(pybamm.Symbol):
         auxiliary_domains=None,
         entries_string=None,
     ):
+        if isinstance(entries, list):
+            entries = np.array(entries)
         if entries.ndim == 1:
             entries = entries[:, np.newaxis]
         if name is None:
diff --git a/pybamm/expression_tree/concatenations.py b/pybamm/expression_tree/concatenations.py
index 611b45773..456029314 100644
--- a/pybamm/expression_tree/concatenations.py
+++ b/pybamm/expression_tree/concatenations.py
@@ -127,7 +127,7 @@ class NumpyConcatenation(Concatenation):
         # so that we can concatenate them
         for i, child in enumerate(children):
             if child.evaluates_to_number():
-                children[i] = child * pybamm.Vector(np.array([1]))
+                children[i] = child * pybamm.Vector([1])
         super().__init__(
             *children,
             name="numpy concatenation",
diff --git a/pybamm/expression_tree/matrix.py b/pybamm/expression_tree/matrix.py
index cd613266d..8aed40bc1 100644
--- a/pybamm/expression_tree/matrix.py
+++ b/pybamm/expression_tree/matrix.py
@@ -2,7 +2,7 @@
 # Matrix class
 #
 import pybamm
-
+import numpy as np
 from scipy.sparse import issparse
 
 
@@ -21,6 +21,8 @@ class Matrix(pybamm.Array):
         auxiliary_domains=None,
         entries_string=None,
     ):
+        if isinstance(entries, list):
+            entries = np.array(entries)
         if name is None:
             name = "Matrix {!s}".format(entries.shape)
             if issparse(entries):
diff --git a/pybamm/expression_tree/state_vector.py b/pybamm/expression_tree/state_vector.py
index 47a6c28e1..7bca90e6c 100644
--- a/pybamm/expression_tree/state_vector.py
+++ b/pybamm/expression_tree/state_vector.py
@@ -26,7 +26,7 @@ class StateVectorBase(pybamm.Symbol):
         List of boolean arrays representing slices. Default is None, in which case the
         evaluation_array is computed from y_slices.
 
-    *Extends:* :class:`Array`
+    *Extends:* :class:`pybamm.Symbol`
     """
 
     def __init__(
@@ -212,7 +212,7 @@ class StateVector(StateVectorBase):
         List of boolean arrays representing slices. Default is None, in which case the
         evaluation_array is computed from y_slices.
 
-    *Extends:* :class:`Array`
+    *Extends:* :class:`pybamm.StateVectorBase`
     """
 
     def __init__(
diff --git a/pybamm/expression_tree/variable.py b/pybamm/expression_tree/variable.py
index 3be2f0b61..44055c651 100644
--- a/pybamm/expression_tree/variable.py
+++ b/pybamm/expression_tree/variable.py
@@ -26,20 +26,33 @@ class VariableBase(pybamm.Symbol):
         collector'. For the DFN, the particle concentration would be a Variable with
         domain 'negative particle', secondary domain 'negative electrode' and tertiary
         domain 'current collector'
+    bounds : tuple, optional
+        Physical bounds on the variable
 
     *Extends:* :class:`Symbol`
     """
 
-    def __init__(self, name, domain=None, auxiliary_domains=None):
+    def __init__(self, name, domain=None, auxiliary_domains=None, bounds=None):
         if domain is None:
             domain = []
         if auxiliary_domains is None:
             auxiliary_domains = {}
         super().__init__(name, domain=domain, auxiliary_domains=auxiliary_domains)
+        if bounds is None:
+            bounds = (-np.inf, np.inf)
+        else:
+            if bounds[0] >= bounds[1]:
+                raise ValueError(
+                    "Invalid bounds {}. ".format(bounds)
+                    + "Lower bound should be strictly less than upper bound."
+                )
+        self.bounds = bounds
 
     def new_copy(self):
         """ See :meth:`pybamm.Symbol.new_copy()`. """
-        return self.__class__(self.name, self.domain, self.auxiliary_domains)
+        return self.__class__(
+            self.name, self.domain, self.auxiliary_domains, self.bounds
+        )
 
     def _evaluate_for_shape(self):
         """ See :meth:`pybamm.Symbol.evaluate_for_shape_using_domain()` """
@@ -59,21 +72,24 @@ class Variable(VariableBase):
 
     name : str
         name of the node
-    domain : iterable of str
+        domain : iterable of str, optional
         list of domains that this variable is valid over
-    auxiliary_domains : dict
+    auxiliary_domains : dict, optional
         dictionary of auxiliary domains ({'secondary': ..., 'tertiary': ...}). For
         example, for the single particle model, the particle concentration would be a
         Variable with domain 'negative particle' and secondary auxiliary domain 'current
         collector'. For the DFN, the particle concentration would be a Variable with
         domain 'negative particle', secondary domain 'negative electrode' and tertiary
         domain 'current collector'
-
+    bounds : tuple, optional
+        Physical bounds on the variable
     *Extends:* :class:`Symbol`
     """
 
-    def __init__(self, name, domain=None, auxiliary_domains=None):
-        super().__init__(name, domain=domain, auxiliary_domains=auxiliary_domains)
+    def __init__(self, name, domain=None, auxiliary_domains=None, bounds=None):
+        super().__init__(
+            name, domain=domain, auxiliary_domains=auxiliary_domains, bounds=bounds
+        )
 
     def diff(self, variable):
         if variable.id == self.id:
@@ -110,11 +126,13 @@ class VariableDot(VariableBase):
         collector'. For the DFN, the particle concentration would be a Variable with
         domain 'negative particle', secondary domain 'negative electrode' and tertiary
         domain 'current collector'
-
+    bounds : tuple, optional
+        Physical bounds on the variable. Included for compatibility with `VariableBase`,
+        but ignored.
     *Extends:* :class:`Symbol`
     """
 
-    def __init__(self, name, domain=None, auxiliary_domains=None):
+    def __init__(self, name, domain=None, auxiliary_domains=None, bounds=None):
         super().__init__(name, domain=domain, auxiliary_domains=auxiliary_domains)
 
     def get_variable(self):
diff --git a/pybamm/expression_tree/vector.py b/pybamm/expression_tree/vector.py
index 7714ab176..a03f177f4 100644
--- a/pybamm/expression_tree/vector.py
+++ b/pybamm/expression_tree/vector.py
@@ -21,6 +21,8 @@ class Vector(pybamm.Array):
         auxiliary_domains=None,
         entries_string=None,
     ):
+        if isinstance(entries, list):
+            entries = np.array(entries)
         # make sure that entries are a vector (can be a column vector)
         if entries.ndim == 1:
             entries = entries[:, np.newaxis]
diff --git a/pybamm/models/standard_variables.py b/pybamm/models/standard_variables.py
index bf0b99558..87e0c7b03 100644
--- a/pybamm/models/standard_variables.py
+++ b/pybamm/models/standard_variables.py
@@ -2,27 +2,33 @@
 # Standard variables for the models
 #
 import pybamm
+import numpy as np
 
 # Electrolyte concentration
 c_e_n = pybamm.Variable(
     "Negative electrolyte concentration",
     domain="negative electrode",
     auxiliary_domains={"secondary": "current collector"},
+    bounds=(0, np.inf),
 )
 c_e_s = pybamm.Variable(
     "Separator electrolyte concentration",
     domain="separator",
     auxiliary_domains={"secondary": "current collector"},
+    bounds=(0, np.inf),
 )
 c_e_p = pybamm.Variable(
     "Positive electrolyte concentration",
     domain="positive electrode",
     auxiliary_domains={"secondary": "current collector"},
+    bounds=(0, np.inf),
 )
 c_e = pybamm.Concatenation(c_e_n, c_e_s, c_e_p)
 
 c_e_av = pybamm.Variable(
-    "X-averaged electrolyte concentration", domain="current collector"
+    "X-averaged electrolyte concentration",
+    domain="current collector",
+    bounds=(0, np.inf),
 )
 
 # Electrolyte potential
@@ -105,6 +111,7 @@ c_s_n = pybamm.Variable(
         "secondary": "negative electrode",
         "tertiary": "current collector",
     },
+    bounds=(0, 1),
 )
 c_s_p = pybamm.Variable(
     "Positive particle concentration",
@@ -113,32 +120,41 @@ c_s_p = pybamm.Variable(
         "secondary": "positive electrode",
         "tertiary": "current collector",
     },
+    bounds=(0, 1),
 )
 c_s_n_xav = pybamm.Variable(
     "X-averaged negative particle concentration",
     domain="negative particle",
     auxiliary_domains={"secondary": "current collector"},
+    bounds=(0, 1),
 )
 c_s_p_xav = pybamm.Variable(
     "X-averaged positive particle concentration",
     domain="positive particle",
     auxiliary_domains={"secondary": "current collector"},
+    bounds=(0, 1),
 )
 c_s_n_surf = pybamm.Variable(
     "Negative particle surface concentration",
     domain="negative electrode",
     auxiliary_domains={"secondary": "current collector"},
+    bounds=(0, 1),
 )
 c_s_p_surf = pybamm.Variable(
     "Positive particle surface concentration",
     domain="positive electrode",
     auxiliary_domains={"secondary": "current collector"},
+    bounds=(0, 1),
 )
 c_s_n_surf_xav = pybamm.Variable(
-    "X-averaged negative particle surface concentration", domain="current collector"
+    "X-averaged negative particle surface concentration",
+    domain="current collector",
+    bounds=(0, 1),
 )
 c_s_p_surf_xav = pybamm.Variable(
-    "X-averaged positive particle surface concentration", domain="current collector"
+    "X-averaged positive particle surface concentration",
+    domain="current collector",
+    bounds=(0, 1),
 )
 
 
@@ -147,26 +163,31 @@ eps_n = pybamm.Variable(
     "Negative electrode porosity",
     domain="negative electrode",
     auxiliary_domains={"secondary": "current collector"},
+    bounds=(0, 1),
 )
 eps_s = pybamm.Variable(
     "Separator porosity",
     domain="separator",
     auxiliary_domains={"secondary": "current collector"},
+    bounds=(0, 1),
 )
 eps_p = pybamm.Variable(
     "Positive electrode porosity",
     domain="positive electrode",
     auxiliary_domains={"secondary": "current collector"},
+    bounds=(0, 1),
 )
 eps = pybamm.Concatenation(eps_n, eps_s, eps_p)
 
 # Piecewise constant (for asymptotic models)
 eps_n_pc = pybamm.Variable(
-    "X-averaged negative electrode porosity", domain="current collector"
+    "X-averaged negative electrode porosity", domain="current collector", bounds=(0, 1),
+)
+eps_s_pc = pybamm.Variable(
+    "X-averaged separator porosity", domain="current collector", bounds=(0, 1)
 )
-eps_s_pc = pybamm.Variable("X-averaged separator porosity", domain="current collector")
 eps_p_pc = pybamm.Variable(
-    "X-averaged positive electrode porosity", domain="current collector"
+    "X-averaged positive electrode porosity", domain="current collector", bounds=(0, 1),
 )
 
 eps_piecewise_constant = pybamm.Concatenation(
@@ -219,4 +240,3 @@ L_outer = pybamm.Variable(
     domain=["negative electrode"],
     auxiliary_domains={"secondary": "current collector"},
 )
-
diff --git a/pybamm/solvers/algebraic_solver.py b/pybamm/solvers/algebraic_solver.py
index da93e50d9..dd4e0d0b2 100644
--- a/pybamm/solvers/algebraic_solver.py
+++ b/pybamm/solvers/algebraic_solver.py
@@ -19,7 +19,9 @@ class AlgebraicSolver(pybamm.BaseSolver):
     Parameters
     ----------
     method : str, optional
-        The method to use to solve the system (default is "lm")
+        The method to use to solve the system (default is "lm"). If it starts with
+        "lsq", least-squares minimization is used. The method for least-squares can be
+        specified in the form "lsq_methodname"
     tol : float, optional
         The tolerance for the solver (default is 1e-6).
     extra_options : dict, optional
@@ -123,14 +125,68 @@ class AlgebraicSolver(pybamm.BaseSolver):
                 y_alg[:, idx] = y0_alg
             # Otherwise calculate new y0
             else:
-                sol = optimize.root(
-                    root_fun,
-                    y0_alg,
-                    method=self.method,
-                    tol=self.tol,
-                    jac=jac_fn,
-                    options=self.extra_options,
-                )
+                # Methods which use least-squares are specified as either "lsq", which
+                # uses the default method, or with "lsq__methodname"
+                if self.method.startswith("lsq"):
+
+                    if self.method == "lsq":
+                        method = "trf"
+                    else:
+                        method = self.method[5:]
+                    if jac_fn is None:
+                        jac_fn = "2-point"
+                    sol = optimize.least_squares(
+                        root_fun,
+                        y0_alg,
+                        method=method,
+                        ftol=self.tol,
+                        jac=jac_fn,
+                        bounds=model.bounds,
+                        **self.extra_options,
+                    )
+                # Methods which use minimize are specified as either "minimize", which
+                # uses the default method, or with "minimize__methodname"
+                elif self.method.startswith("minimize"):
+                    # Adapt the root function for minimize
+                    def root_norm(y):
+                        return np.sum(root_fun(y) ** 2)
+
+                    if jac_fn is None:
+                        jac_norm = None
+                    else:
+
+                        def jac_norm(y):
+                            return np.sum(2 * root_fun(y) * jac_fn(y), 0)
+
+                    if self.method == "minimize":
+                        method = None
+                    else:
+                        method = self.method[10:]
+                    extra_options = self.extra_options
+                    if np.any(model.bounds[0] != -np.inf) or np.any(
+                        model.bounds[1] != np.inf
+                    ):
+                        bounds = [
+                            (lb, ub) for lb, ub in zip(model.bounds[0], model.bounds[1])
+                        ]
+                        extra_options["bounds"] = bounds
+                    sol = optimize.minimize(
+                        root_norm,
+                        y0_alg,
+                        method=method,
+                        tol=self.tol,
+                        jac=jac_norm,
+                        **extra_options,
+                    )
+                else:
+                    sol = optimize.root(
+                        root_fun,
+                        y0_alg,
+                        method=self.method,
+                        tol=self.tol,
+                        jac=jac_fn,
+                        options=self.extra_options,
+                    )
 
                 if sol.success and np.all(abs(sol.fun) < self.tol):
                     # update initial guess for the next iteration
@@ -143,12 +199,10 @@ class AlgebraicSolver(pybamm.BaseSolver):
                     )
                 else:
                     raise pybamm.SolverError(
-                        """
-                        Could not find acceptable solution: solver terminated
-                        successfully, but maximum solution error ({})
-                        above tolerance ({})
-                        """.format(
-                            np.max(sol.fun), self.tol
+                        "Could not find acceptable solution: solver terminated "
+                        "successfully, but maximum solution error "
+                        "({}) above tolerance ({})".format(
+                            np.max(abs(sol.fun)), self.tol
                         )
                     )
 
diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index 159118bf7..321709781 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -51,7 +51,7 @@ class BaseSolver(object):
                 "max_steps has been deprecated, and should be set using the "
                 "solver-specific extra-options dictionaries instead"
             )
-        self.models_set_up = set()
+        self.models_set_up = {}
 
         # Defaults, can be overwritten by specific solver
         self.name = "Base solver"
@@ -114,7 +114,7 @@ class BaseSolver(object):
         "Returns a copy of the solver"
         new_solver = copy.copy(self)
         # clear models_set_up
-        new_solver.models_set_up = set()
+        new_solver.models_set_up = {}
         return new_solver
 
     def set_up(self, model, inputs=None):
@@ -504,9 +504,23 @@ class BaseSolver(object):
         if model not in self.models_set_up:
             self.set_up(model, ext_and_inputs)
             set_up_time = timer.time()
-            self.models_set_up.add(model)
+            self.models_set_up.update(
+                {model: {"initial conditions": model.concatenated_initial_conditions}}
+            )
         else:
-            set_up_time = 0
+            ics_set_up = self.models_set_up[model]["initial conditions"]
+            # Check that initial conditions have not been updated
+            if ics_set_up.id == model.concatenated_initial_conditions.id:
+                set_up_time = 0
+            else:
+                # If the new initial conditions are different, set up again
+                # Doing the whole setup again might be slow, but no need to prematurely
+                # optimize this
+                self.set_up(model, ext_and_inputs)
+                self.models_set_up[model][
+                    "initial conditions"
+                ] = model.concatenated_initial_conditions
+                set_up_time = timer.time()
 
         # (Re-)calculate consistent initial conditions
         self._set_initial_conditions(model, ext_and_inputs, update_rhs=True)
diff --git a/pybamm/solvers/casadi_algebraic_solver.py b/pybamm/solvers/casadi_algebraic_solver.py
index 229cc7f6d..e286936f2 100644
--- a/pybamm/solvers/casadi_algebraic_solver.py
+++ b/pybamm/solvers/casadi_algebraic_solver.py
@@ -68,6 +68,7 @@ class CasadiAlgebraicSolver(pybamm.BaseSolver):
         # equations will be equal to the initial condition provided. This allows this
         # solver to be used for initialising the DAE solvers
         if model.rhs == {}:
+            len_rhs = 0
             y0_diff = casadi.DM()
             y0_alg = y0
         else:
@@ -86,12 +87,25 @@ class CasadiAlgebraicSolver(pybamm.BaseSolver):
         t_p_sym = casadi.vertcat(t_sym, p_sym)
         alg = model.casadi_algebraic(t_sym, y_sym, p_sym)
 
+        # Set constraints vector in the casadi format
+        # Constrain the unknowns. 0 (default): no constraint on ui, 1: ui >= 0.0,
+        # -1: ui <= 0.0, 2: ui > 0.0, -2: ui < 0.0.
+        constraints = np.zeros_like(model.bounds[0], dtype=int)
+        # If the lower bound is positive then the variable must always be positive
+        constraints[model.bounds[0] >= 0] = 1
+        # If the upper bound is negative then the variable must always be negative
+        constraints[model.bounds[1] <= 0] = -1
+
         # Set up rootfinder
         roots = casadi.rootfinder(
             "roots",
             "newton",
             dict(x=y_alg_sym, p=t_p_sym, g=alg),
-            {**self.extra_options, "abstol": self.tol},
+            {
+                **self.extra_options,
+                "abstol": self.tol,
+                "constraints": list(constraints[len_rhs:]),
+            },
         )
         for idx, t in enumerate(t_eval):
             # Evaluate algebraic with new t and previous y0, if it's already close
@@ -146,7 +160,7 @@ class CasadiAlgebraicSolver(pybamm.BaseSolver):
                         successfully, but maximum solution error ({})
                         above tolerance ({})
                         """.format(
-                            casadi.mmax(fun), self.tol
+                            casadi.mmax(casadi.fabs(fun)), self.tol
                         )
                     )
 
diff --git a/pybamm/spatial_methods/finite_volume.py b/pybamm/spatial_methods/finite_volume.py
index c2e9683cb..8b00a6bfc 100644
--- a/pybamm/spatial_methods/finite_volume.py
+++ b/pybamm/spatial_methods/finite_volume.py
@@ -641,7 +641,7 @@ class FiniteVolume(pybamm.SpatialMethod):
         n = submesh.npts
         second_dim_repeats = self._get_auxiliary_domain_repeats(symbol.domains)
 
-        bcs_vector = pybamm.Vector(np.array([]))  # starts empty
+        bcs_vector = pybamm.Vector([])  # starts empty
 
         lbc_value, lbc_type = bcs["left"]
         rbc_value, rbc_type = bcs["right"]


Test patch:
diff --git a/tests/unit/test_discretisations/test_discretisation.py b/tests/unit/test_discretisations/test_discretisation.py
index f827a5a13..7ceff8746 100644
--- a/tests/unit/test_discretisations/test_discretisation.py
+++ b/tests/unit/test_discretisations/test_discretisation.py
@@ -131,6 +131,7 @@ class TestDiscretise(unittest.TestCase):
         self.assertEqual(disc.y_slices[a.id][0], slice(0, 10, None))
 
         self.assertEqual(model.y_slices[a][0], slice(0, 10, None))
+        self.assertEqual(model.bounds, disc.bounds)
 
         b_test = np.ones((10, 1))
         np.testing.assert_array_equal(
@@ -265,7 +266,7 @@ class TestDiscretise(unittest.TestCase):
         np.testing.assert_array_equal(y[disc.y_slices[c.id][0]], c_true)
 
         # Several variables
-        d = pybamm.Variable("d", domain=whole_cell)
+        d = pybamm.Variable("d", domain=whole_cell, bounds=(0, 1))
         jn = pybamm.Variable("jn", domain=["negative electrode"])
         variables = [c, d, jn]
         disc.set_variable_slices(variables)
@@ -274,6 +275,12 @@ class TestDiscretise(unittest.TestCase):
             disc.y_slices,
             {c.id: [slice(0, 100)], d.id: [slice(100, 200)], jn.id: [slice(200, 240)]},
         )
+        np.testing.assert_array_equal(
+            disc.bounds[0], [-np.inf] * 100 + [0] * 100 + [-np.inf] * 40
+        )
+        np.testing.assert_array_equal(
+            disc.bounds[1], [np.inf] * 100 + [1] * 100 + [np.inf] * 40
+        )
         d_true = 4 * combined_submesh.nodes
         jn_true = mesh["negative electrode"].nodes ** 3
         y = np.concatenate([c_true, d_true, jn_true])
@@ -297,6 +304,12 @@ class TestDiscretise(unittest.TestCase):
                 jp.id: [slice(265, 300)],
             },
         )
+        np.testing.assert_array_equal(
+            disc.bounds[0], [-np.inf] * 100 + [0] * 100 + [-np.inf] * 100
+        )
+        np.testing.assert_array_equal(
+            disc.bounds[1], [np.inf] * 100 + [1] * 100 + [np.inf] * 100
+        )
         d_true = 4 * combined_submesh.nodes
         jn_true = mesh["negative electrode"].nodes ** 3
         y = np.concatenate([c_true, d_true, jn_true])
@@ -338,12 +351,12 @@ class TestDiscretise(unittest.TestCase):
         self.assertIsInstance(scal_disc, pybamm.Scalar)
         self.assertEqual(scal_disc.value, scal.value)
         # vector
-        vec = pybamm.Vector(np.array([1, 2, 3, 4]))
+        vec = pybamm.Vector([1, 2, 3, 4])
         vec_disc = disc.process_symbol(vec)
         self.assertIsInstance(vec_disc, pybamm.Vector)
         np.testing.assert_array_equal(vec_disc.entries, vec.entries)
         # matrix
-        mat = pybamm.Matrix(np.array([[1, 2, 3, 4], [5, 6, 7, 8]]))
+        mat = pybamm.Matrix([[1, 2, 3, 4], [5, 6, 7, 8]])
         mat_disc = disc.process_symbol(mat)
         self.assertIsInstance(mat_disc, pybamm.Matrix)
         np.testing.assert_array_equal(mat_disc.entries, mat.entries)
@@ -839,39 +852,39 @@ class TestDiscretise(unittest.TestCase):
         with self.assertRaises(pybamm.ModelError):
             disc.process_model(model)
 
-    # def test_process_model_concatenation(self):
-    #     # concatenation of variables as the key
-    #     cn = pybamm.Variable("c", domain=["negative electrode"])
-    #     cs = pybamm.Variable("c", domain=["separator"])
-    #     cp = pybamm.Variable("c", domain=["positive electrode"])
-    #     c = pybamm.Concatenation(cn, cs, cp)
-    #     N = pybamm.grad(c)
-    #     model = pybamm.BaseModel()
-    #     model.rhs = {c: pybamm.div(N)}
-    #     model.initial_conditions = {c: pybamm.Scalar(3)}
-
-    #     model.boundary_conditions = {
-    #         c: {"left": (0, "Neumann"), "right": (0, "Neumann")}
-    #     }
-    #     model.check_well_posedness()
-
-    #     # create discretisation
-    #     disc = get_discretisation_for_testing()
-    #     mesh = disc.mesh
-
-    #     combined_submesh = mesh.combine_submeshes(
-    #         "negative electrode", "separator", "positive electrode"
-    #     )
-
-    #     disc.process_model(model)
-    #     y0 = model.concatenated_initial_conditions.evaluate()
-    #     np.testing.assert_array_equal(
-    #         y0, 3 * np.ones_like(combined_submesh.nodes[:, np.newaxis])
-    #     )
-
-    #     # grad and div are identity operators here
-    #     np.testing.assert_array_equal(y0, model.concatenated_rhs.evaluate(None, y0))
-    #     model.check_well_posedness()
+    def test_process_model_concatenation(self):
+        # concatenation of variables as the key
+        cn = pybamm.Variable("c", domain=["negative electrode"])
+        cs = pybamm.Variable("c", domain=["separator"])
+        cp = pybamm.Variable("c", domain=["positive electrode"])
+        c = pybamm.Concatenation(cn, cs, cp)
+        N = pybamm.grad(c)
+        model = pybamm.BaseModel()
+        model.rhs = {c: pybamm.div(N)}
+        model.initial_conditions = {c: pybamm.Scalar(3)}
+
+        model.boundary_conditions = {
+            c: {"left": (0, "Neumann"), "right": (0, "Neumann")}
+        }
+        model.check_well_posedness()
+
+        # create discretisation
+        disc = get_discretisation_for_testing()
+        mesh = disc.mesh
+
+        combined_submesh = mesh.combine_submeshes(
+            "negative electrode", "separator", "positive electrode"
+        )
+
+        disc.process_model(model)
+        y0 = model.concatenated_initial_conditions.evaluate()
+        np.testing.assert_array_equal(
+            y0, 3 * np.ones_like(combined_submesh.nodes[:, np.newaxis])
+        )
+
+        # grad and div are identity operators here
+        np.testing.assert_array_equal(y0, model.concatenated_rhs.evaluate(None, y0))
+        model.check_well_posedness()
 
     def test_process_model_not_inplace(self):
         # concatenation of variables as the key
@@ -1060,16 +1073,19 @@ class TestDiscretise(unittest.TestCase):
             "a",
             domain=["negative electrode"],
             auxiliary_domains={"secondary": "current collector"},
+            bounds=(-5, -2),
         )
         b = pybamm.Variable(
             "b",
             domain=["separator"],
             auxiliary_domains={"secondary": "current collector"},
+            bounds=(6, 10),
         )
         c = pybamm.Variable(
             "c",
             domain=["positive electrode"],
             auxiliary_domains={"secondary": "current collector"},
+            bounds=(0, 1),
         )
         conc = pybamm.Concatenation(a, b, c)
         disc.set_variable_slices([conc])
@@ -1082,6 +1098,12 @@ class TestDiscretise(unittest.TestCase):
         self.assertEqual(
             disc.y_slices[c.id], [slice(65, 100), slice(165, 200), slice(265, 300)]
         )
+        np.testing.assert_array_equal(
+            disc.bounds[0], ([-5] * 40 + [6] * 25 + [0] * 35) * 3
+        )
+        np.testing.assert_array_equal(
+            disc.bounds[1], ([-2] * 40 + [10] * 25 + [1] * 35) * 3
+        )
         expr = disc.process_symbol(conc)
         self.assertIsInstance(expr, pybamm.DomainConcatenation)
 
diff --git a/tests/unit/test_expression_tree/test_array.py b/tests/unit/test_expression_tree/test_array.py
index 550abd84c..e2adab03b 100644
--- a/tests/unit/test_expression_tree/test_array.py
+++ b/tests/unit/test_expression_tree/test_array.py
@@ -12,6 +12,12 @@ class TestArray(unittest.TestCase):
         arr = pybamm.Array(np.array([1, 2, 3]))
         self.assertEqual(arr.name, "Array of shape (3, 1)")
 
+    def test_list_entries(self):
+        vect = pybamm.Array([1, 2, 3])
+        np.testing.assert_array_equal(vect.entries, np.array([[1], [2], [3]]))
+        vect = pybamm.Array([[1], [2], [3]])
+        np.testing.assert_array_equal(vect.entries, np.array([[1], [2], [3]]))
+
     def test_linspace(self):
         x = np.linspace(0, 1, 100)[:, np.newaxis]
         y = pybamm.linspace(0, 1, 100)
diff --git a/tests/unit/test_expression_tree/test_binary_operators.py b/tests/unit/test_expression_tree/test_binary_operators.py
index 02a4b5cea..71d1689c3 100644
--- a/tests/unit/test_expression_tree/test_binary_operators.py
+++ b/tests/unit/test_expression_tree/test_binary_operators.py
@@ -336,7 +336,7 @@ class TestIsZero(unittest.TestCase):
 
         a = pybamm.Matrix(np.zeros((10, 10)))
         b = pybamm.Matrix(np.ones((10, 10)))
-        c = pybamm.Matrix(np.array([1, 0, 0]))
+        c = pybamm.Matrix([1, 0, 0])
         self.assertTrue(pybamm.is_matrix_zero(a))
         self.assertFalse(pybamm.is_matrix_zero(b))
         self.assertFalse(pybamm.is_matrix_zero(c))
diff --git a/tests/unit/test_expression_tree/test_concatenations.py b/tests/unit/test_expression_tree/test_concatenations.py
index 3e2e42b2b..e87f2aaa2 100644
--- a/tests/unit/test_expression_tree/test_concatenations.py
+++ b/tests/unit/test_expression_tree/test_concatenations.py
@@ -18,9 +18,9 @@ class TestConcatenations(unittest.TestCase):
         self.assertEqual(conc.children[0].name, "a")
         self.assertEqual(conc.children[1].name, "b")
         self.assertEqual(conc.children[2].name, "c")
-        d = pybamm.Vector(np.array([2]))
-        e = pybamm.Vector(np.array([1]))
-        f = pybamm.Vector(np.array([3]))
+        d = pybamm.Vector([2])
+        e = pybamm.Vector([1])
+        f = pybamm.Vector([3])
         conc2 = pybamm.Concatenation(d, e, f)
         with self.assertRaises(TypeError):
             conc2.evaluate()
diff --git a/tests/unit/test_expression_tree/test_matrix.py b/tests/unit/test_expression_tree/test_matrix.py
index 0d2300ed5..f4eb8988d 100644
--- a/tests/unit/test_expression_tree/test_matrix.py
+++ b/tests/unit/test_expression_tree/test_matrix.py
@@ -19,6 +19,12 @@ class TestMatrix(unittest.TestCase):
         self.assertEqual(self.mat.shape, (3, 3))
         self.assertEqual(self.mat.size, 9)
 
+    def test_list_entry(self):
+        mat = pybamm.Matrix([[1, 2, 0], [0, 1, 0], [0, 0, 1]])
+        np.testing.assert_array_equal(
+            mat.entries, np.array([[1, 2, 0], [0, 1, 0], [0, 0, 1]])
+        )
+
     def test_matrix_evaluate(self):
         np.testing.assert_array_equal(self.mat.evaluate(), self.A)
 
diff --git a/tests/unit/test_expression_tree/test_operations/test_copy.py b/tests/unit/test_expression_tree/test_operations/test_copy.py
index ef6d3a063..564eee5e6 100644
--- a/tests/unit/test_expression_tree/test_operations/test_copy.py
+++ b/tests/unit/test_expression_tree/test_operations/test_copy.py
@@ -14,7 +14,7 @@ class TestCopy(unittest.TestCase):
         v_n = pybamm.Variable("v", "negative electrode")
         x_n = pybamm.standard_spatial_vars.x_n
         v_s = pybamm.Variable("v", "separator")
-        vec = pybamm.Vector(np.array([1, 2, 3, 4, 5]))
+        vec = pybamm.Vector([1, 2, 3, 4, 5])
         mesh = get_mesh_for_testing()
 
         for symbol in [
diff --git a/tests/unit/test_expression_tree/test_operations/test_evaluate.py b/tests/unit/test_expression_tree/test_operations/test_evaluate.py
index 8ff7eaaf4..3f08f1137 100644
--- a/tests/unit/test_expression_tree/test_operations/test_evaluate.py
+++ b/tests/unit/test_expression_tree/test_operations/test_evaluate.py
@@ -106,7 +106,7 @@ class TestEvaluate(unittest.TestCase):
         # test matrix
         constant_symbols = OrderedDict()
         variable_symbols = OrderedDict()
-        A = pybamm.Matrix(np.array([[1, 2], [3, 4]]))
+        A = pybamm.Matrix([[1, 2], [3, 4]])
         pybamm.find_symbols(A, constant_symbols, variable_symbols)
         self.assertEqual(len(variable_symbols), 0)
         self.assertEqual(list(constant_symbols.keys())[0], A.id)
@@ -335,7 +335,7 @@ class TestEvaluate(unittest.TestCase):
             self.assertEqual(result, expr.evaluate(t=t, y=y))
 
         # test something with a matrix multiplication
-        A = pybamm.Matrix(np.array([[1, 2], [3, 4]]))
+        A = pybamm.Matrix([[1, 2], [3, 4]])
         expr = A @ pybamm.StateVector(slice(0, 2))
         evaluator = pybamm.EvaluatorPython(expr)
         for t, y in zip(t_tests, y_tests):
@@ -343,7 +343,7 @@ class TestEvaluate(unittest.TestCase):
             np.testing.assert_allclose(result, expr.evaluate(t=t, y=y))
 
         # test something with a heaviside
-        a = pybamm.Vector(np.array([1, 2]))
+        a = pybamm.Vector([1, 2])
         expr = a <= pybamm.StateVector(slice(0, 2))
         evaluator = pybamm.EvaluatorPython(expr)
         for t, y in zip(t_tests, y_tests):
@@ -357,7 +357,7 @@ class TestEvaluate(unittest.TestCase):
             np.testing.assert_allclose(result, expr.evaluate(t=t, y=y))
 
         # test something with a minimum or maximum
-        a = pybamm.Vector(np.array([1, 2]))
+        a = pybamm.Vector([1, 2])
         expr = pybamm.minimum(a, pybamm.StateVector(slice(0, 2)))
         evaluator = pybamm.EvaluatorPython(expr)
         for t, y in zip(t_tests, y_tests):
@@ -378,7 +378,7 @@ class TestEvaluate(unittest.TestCase):
             self.assertEqual(result, expr.evaluate(t=t, y=y))
 
         # test something with a sparse matrix multiplication
-        A = pybamm.Matrix(np.array([[1, 2], [3, 4]]))
+        A = pybamm.Matrix([[1, 2], [3, 4]])
         B = pybamm.Matrix(scipy.sparse.csr_matrix(np.array([[1, 0], [0, 4]])))
         C = pybamm.Matrix(scipy.sparse.coo_matrix(np.array([[1, 0], [0, 4]])))
         expr = A @ B @ C @ pybamm.StateVector(slice(0, 2))
@@ -388,8 +388,8 @@ class TestEvaluate(unittest.TestCase):
             np.testing.assert_allclose(result, expr.evaluate(t=t, y=y))
 
         # test numpy concatenation
-        a = pybamm.Vector(np.array([[1], [2]]))
-        b = pybamm.Vector(np.array([[3]]))
+        a = pybamm.Vector([[1], [2]])
+        b = pybamm.Vector([[3]])
         expr = pybamm.NumpyConcatenation(a, b)
         evaluator = pybamm.EvaluatorPython(expr)
         for t, y in zip(t_tests, y_tests):
diff --git a/tests/unit/test_expression_tree/test_operations/test_jac_2D.py b/tests/unit/test_expression_tree/test_operations/test_jac_2D.py
index 3bd80db91..5fea85ee7 100644
--- a/tests/unit/test_expression_tree/test_operations/test_jac_2D.py
+++ b/tests/unit/test_expression_tree/test_operations/test_jac_2D.py
@@ -222,6 +222,10 @@ class TestJacobian(unittest.TestCase):
         a = 2 * pybamm.PrimaryBroadcast(curr_coll_vector, a_dom)
         b = pybamm.PrimaryBroadcast(curr_coll_vector, b_dom)
         c = 3 * pybamm.PrimaryBroadcast(curr_coll_vector, c_dom)
+        # Add bounds for compatibility with the discretisation
+        a.bounds = (-np.inf, np.inf)
+        b.bounds = (-np.inf, np.inf)
+        c.bounds = (-np.inf, np.inf)
 
         conc = pybamm.Concatenation(a, b, c)
         disc.set_variable_slices([conc])
diff --git a/tests/unit/test_expression_tree/test_operations/test_simplify.py b/tests/unit/test_expression_tree/test_operations/test_simplify.py
index b74259d89..3aac02b6f 100644
--- a/tests/unit/test_expression_tree/test_operations/test_simplify.py
+++ b/tests/unit/test_expression_tree/test_operations/test_simplify.py
@@ -336,7 +336,7 @@ class TestSimplify(unittest.TestCase):
         b = pybamm.Matrix(np.ones((2, 2)))
 
         # matrix multiplication
-        A = pybamm.Matrix(np.array([[1, 0], [0, 1]]))
+        A = pybamm.Matrix([[1, 0], [0, 1]])
         self.assertIsInstance((a @ A).simplify(), pybamm.Matrix)
         np.testing.assert_array_equal(
             (a @ A).simplify().evaluate().toarray(), np.zeros((2, 2))
@@ -347,8 +347,8 @@ class TestSimplify(unittest.TestCase):
         )
 
         # matrix * matrix
-        m1 = pybamm.Matrix(np.array([[2, 0], [0, 2]]))
-        m2 = pybamm.Matrix(np.array([[3, 0], [0, 3]]))
+        m1 = pybamm.Matrix([[2, 0], [0, 2]])
+        m2 = pybamm.Matrix([[3, 0], [0, 3]])
         v = pybamm.StateVector(slice(0, 2))
         v2 = pybamm.StateVector(slice(2, 4))
 
@@ -412,8 +412,8 @@ class TestSimplify(unittest.TestCase):
             )
 
         # matrix * vector
-        m1 = pybamm.Matrix(np.array([[2, 0], [0, 2]]))
-        v1 = pybamm.Vector(np.array([1, 1]))
+        m1 = pybamm.Matrix([[2, 0], [0, 2]])
+        v1 = pybamm.Vector([1, 1])
 
         for expr in [(m1 @ v1).simplify()]:
             self.assertIsInstance(expr, pybamm.Vector)
diff --git a/tests/unit/test_expression_tree/test_variable.py b/tests/unit/test_expression_tree/test_variable.py
index 8e4a8be11..b023f1e43 100644
--- a/tests/unit/test_expression_tree/test_variable.py
+++ b/tests/unit/test_expression_tree/test_variable.py
@@ -33,6 +33,18 @@ class TestVariable(unittest.TestCase):
         self.assertNotEqual(a1.id, a3.id)
         self.assertNotEqual(a1.id, a4.id)
 
+    def test_variable_bounds(self):
+        var = pybamm.Variable("var")
+        self.assertEqual(var.bounds, (-np.inf, np.inf))
+
+        var = pybamm.Variable("var", bounds=(0, 1))
+        self.assertEqual(var.bounds, (0, 1))
+
+        with self.assertRaisesRegex(ValueError, "Invalid bounds"):
+            pybamm.Variable("var", bounds=(1, 0))
+        with self.assertRaisesRegex(ValueError, "Invalid bounds"):
+            pybamm.Variable("var", bounds=(1, 1))
+
 
 class TestVariableDot(unittest.TestCase):
     def test_variable_init(self):
diff --git a/tests/unit/test_expression_tree/test_vector.py b/tests/unit/test_expression_tree/test_vector.py
index 924876dd8..52c414b78 100644
--- a/tests/unit/test_expression_tree/test_vector.py
+++ b/tests/unit/test_expression_tree/test_vector.py
@@ -21,6 +21,12 @@ class TestVector(unittest.TestCase):
         vect1d = pybamm.Vector(np.array([1, 2, 3]))
         np.testing.assert_array_equal(self.vect.entries, vect1d.entries)
 
+    def test_list_entries(self):
+        vect = pybamm.Vector([1, 2, 3])
+        np.testing.assert_array_equal(vect.entries, np.array([[1], [2], [3]]))
+        vect = pybamm.Vector([[1], [2], [3]])
+        np.testing.assert_array_equal(vect.entries, np.array([[1], [2], [3]]))
+
     def test_vector_evaluate(self):
         np.testing.assert_array_equal(self.vect.evaluate(), self.x)
 
diff --git a/tests/unit/test_solvers/test_algebraic_solver.py b/tests/unit/test_solvers/test_algebraic_solver.py
index 45c858327..4e7d660b5 100644
--- a/tests/unit/test_solvers/test_algebraic_solver.py
+++ b/tests/unit/test_solvers/test_algebraic_solver.py
@@ -135,6 +135,7 @@ class TestAlgebraicSolver(unittest.TestCase):
 
         # Test without jacobian
         model.use_jacobian = False
+        solver.models_set_up = set()
         solution_no_jac = solver.solve(model)
         np.testing.assert_array_equal(
             model.variables["var1"].evaluate(t=None, y=solution_no_jac.y), sol[:100]
@@ -143,6 +144,112 @@ class TestAlgebraicSolver(unittest.TestCase):
             model.variables["var2"].evaluate(t=None, y=solution_no_jac.y), sol[100:]
         )
 
+    def test_model_solver_least_squares(self):
+        # Create model
+        model = pybamm.BaseModel()
+        whole_cell = ["negative electrode", "separator", "positive electrode"]
+        var1 = pybamm.Variable("var1", domain=whole_cell)
+        var2 = pybamm.Variable("var2", domain=whole_cell)
+        model.algebraic = {var1: var1 - 3, var2: 2 * var1 - var2}
+        model.initial_conditions = {var1: pybamm.Scalar(1), var2: pybamm.Scalar(4)}
+        model.variables = {"var1": var1, "var2": var2}
+        disc = get_discretisation_for_testing()
+        disc.process_model(model)
+
+        sol = np.concatenate((np.ones(100) * 3, np.ones(100) * 6))[:, np.newaxis]
+
+        # Solve
+        solver = pybamm.AlgebraicSolver("lsq")
+        solution = solver.solve(model)
+        np.testing.assert_array_almost_equal(
+            model.variables["var1"].evaluate(t=None, y=solution.y), sol[:100]
+        )
+        np.testing.assert_array_almost_equal(
+            model.variables["var2"].evaluate(t=None, y=solution.y), sol[100:]
+        )
+
+        # Test without jacobian and with a different method
+        model.use_jacobian = False
+        solver = pybamm.AlgebraicSolver("lsq__trf")
+        solution_no_jac = solver.solve(model)
+        np.testing.assert_array_almost_equal(
+            model.variables["var1"].evaluate(t=None, y=solution_no_jac.y), sol[:100]
+        )
+        np.testing.assert_array_almost_equal(
+            model.variables["var2"].evaluate(t=None, y=solution_no_jac.y), sol[100:]
+        )
+
+    def test_model_solver_minimize(self):
+        # Create model
+        model = pybamm.BaseModel()
+        whole_cell = ["negative electrode", "separator", "positive electrode"]
+        var1 = pybamm.Variable("var1", domain=whole_cell)
+        var2 = pybamm.Variable("var2", domain=whole_cell)
+        model.algebraic = {var1: var1 - 3, var2: 2 * var1 - var2}
+        model.initial_conditions = {var1: pybamm.Scalar(1), var2: pybamm.Scalar(4)}
+        model.variables = {"var1": var1, "var2": var2}
+        disc = get_discretisation_for_testing()
+        disc.process_model(model)
+
+        sol = np.concatenate((np.ones(100) * 3, np.ones(100) * 6))[:, np.newaxis]
+
+        # Solve
+        solver = pybamm.AlgebraicSolver("minimize", tol=1e-8)
+        solution = solver.solve(model)
+        np.testing.assert_array_almost_equal(
+            model.variables["var1"].evaluate(t=None, y=solution.y), sol[:100]
+        )
+        np.testing.assert_array_almost_equal(
+            model.variables["var2"].evaluate(t=None, y=solution.y), sol[100:]
+        )
+
+        # Test without jacobian and with a different method
+        model.use_jacobian = False
+        solver = pybamm.AlgebraicSolver("minimize__BFGS")
+        solution_no_jac = solver.solve(model)
+        np.testing.assert_array_almost_equal(
+            model.variables["var1"].evaluate(t=None, y=solution_no_jac.y), sol[:100]
+        )
+        np.testing.assert_array_almost_equal(
+            model.variables["var2"].evaluate(t=None, y=solution_no_jac.y), sol[100:]
+        )
+
+    def test_model_solver_least_squares_with_bounds(self):
+        # Note: we need a better test case to test this functionality properly
+        # Create model
+        model = pybamm.BaseModel()
+        var1 = pybamm.Variable("var1", bounds=(0, 10))
+        model.algebraic = {var1: pybamm.sin(var1) + 1}
+        model.initial_conditions = {var1: pybamm.Scalar(3)}
+        model.variables = {"var1": var1}
+
+        # Solve
+        solver = pybamm.AlgebraicSolver("lsq", tol=1e-5)
+        solution = solver.solve(model)
+        np.testing.assert_array_almost_equal(
+            model.variables["var1"].evaluate(t=None, y=solution.y),
+            3 * np.pi / 2,
+            decimal=2,
+        )
+
+    def test_model_solver_minimize_with_bounds(self):
+        # Note: we need a better test case to test this functionality properly
+        # Create model
+        model = pybamm.BaseModel()
+        var1 = pybamm.Variable("var1", bounds=(0, 10))
+        model.algebraic = {var1: pybamm.sin(var1) + 1}
+        model.initial_conditions = {var1: pybamm.Scalar(3)}
+        model.variables = {"var1": var1}
+
+        # Solve
+        solver = pybamm.AlgebraicSolver("minimize", tol=1e-16)
+        solution = solver.solve(model)
+        np.testing.assert_array_almost_equal(
+            model.variables["var1"].evaluate(t=None, y=solution.y),
+            3 * np.pi / 2,
+            decimal=4,
+        )
+
     def test_model_solver_with_time(self):
         # Create model
         model = pybamm.BaseModel()
diff --git a/tests/unit/test_solvers/test_base_solver.py b/tests/unit/test_solvers/test_base_solver.py
index 3633cf6bf..35e6cf25d 100644
--- a/tests/unit/test_solvers/test_base_solver.py
+++ b/tests/unit/test_solvers/test_base_solver.py
@@ -117,6 +117,7 @@ class TestBaseSolver(unittest.TestCase):
                     "alg", [t, y, p], [self.algebraic_eval(t, y, p)]
                 )
                 self.convert_to_format = "casadi"
+                self.bounds = (np.array([-np.inf]), np.array([np.inf]))
 
             def rhs_eval(self, t, y, inputs):
                 return np.array([])
@@ -151,6 +152,7 @@ class TestBaseSolver(unittest.TestCase):
                     "alg", [t, y, p], [self.algebraic_eval(t, y, p)]
                 )
                 self.convert_to_format = "casadi"
+                self.bounds = (-np.inf * np.ones(4), np.inf * np.ones(4))
 
             def rhs_eval(self, t, y, inputs):
                 return y[0:1]
@@ -197,6 +199,7 @@ class TestBaseSolver(unittest.TestCase):
                     "alg", [t, y, p], [self.algebraic_eval(t, y, p)]
                 )
                 self.convert_to_format = "casadi"
+                self.bounds = (np.array([-np.inf]), np.array([np.inf]))
 
             def rhs_eval(self, t, y, inputs):
                 return np.array([])
diff --git a/tests/unit/test_solvers/test_casadi_algebraic_solver.py b/tests/unit/test_solvers/test_casadi_algebraic_solver.py
index 208860d04..1d8422318 100644
--- a/tests/unit/test_solvers/test_casadi_algebraic_solver.py
+++ b/tests/unit/test_solvers/test_casadi_algebraic_solver.py
@@ -58,6 +58,7 @@ class TestCasadiAlgebraicSolver(unittest.TestCase):
             p = casadi.MX.sym("p")
             rhs = {}
             casadi_algebraic = casadi.Function("alg", [t, y, p], [y ** 2 + 1])
+            bounds = (np.array([-np.inf]), np.array([np.inf]))
 
             def algebraic_eval(self, t, y, inputs):
                 # algebraic equation has no real root
@@ -104,6 +105,22 @@ class TestCasadiAlgebraicSolver(unittest.TestCase):
             sol[1, :],
         )
 
+    def test_model_solver_with_bounds(self):
+        # Note: we need a better test case to test this functionality properly
+        # Create model
+        model = pybamm.BaseModel()
+        var1 = pybamm.Variable("var1", bounds=(0, 10))
+        model.algebraic = {var1: pybamm.sin(var1) + 1}
+        model.initial_conditions = {var1: pybamm.Scalar(3)}
+        model.variables = {"var1": var1}
+
+        # Solve
+        solver = pybamm.CasadiAlgebraicSolver(tol=1e-12)
+        solution = solver.solve(model)
+        np.testing.assert_array_almost_equal(
+            model.variables["var1"].evaluate(t=None, y=solution.y), 3 * np.pi / 2
+        )
+
     def test_solve_with_input(self):
         # Simple system: a single algebraic equation
         var = pybamm.Variable("var")
diff --git a/tests/unit/test_solvers/test_scikits_solvers.py b/tests/unit/test_solvers/test_scikits_solvers.py
index c984d7f40..6ab8fe2dc 100644
--- a/tests/unit/test_solvers/test_scikits_solvers.py
+++ b/tests/unit/test_solvers/test_scikits_solvers.py
@@ -42,7 +42,7 @@ class TestScikitsSolvers(unittest.TestCase):
 
         # Create custom model so that custom ics
         class Model:
-            mass_matrix = pybamm.Matrix(np.array([[1.0, 0.0], [0.0, 0.0]]))
+            mass_matrix = pybamm.Matrix([[1.0, 0.0], [0.0, 0.0]])
             y0 = np.array([0.0, 1.0])
             terminate_events_eval = []
             timescale_eval = 1
@@ -92,7 +92,7 @@ class TestScikitsSolvers(unittest.TestCase):
 
         # Create custom model so that custom mass matrix can be used
         class Model:
-            mass_matrix = pybamm.Matrix(np.array([[4.0, 0.0], [0.0, 0.0]]))
+            mass_matrix = pybamm.Matrix([[4.0, 0.0], [0.0, 0.0]])
             y0 = np.array([0.0, 0.0])
             terminate_events_eval = []
             timescale_eval = 1
diff --git a/tests/unit/test_solvers/test_scipy_solver.py b/tests/unit/test_solvers/test_scipy_solver.py
index 1d6cfd43a..1d0aae34b 100644
--- a/tests/unit/test_solvers/test_scipy_solver.py
+++ b/tests/unit/test_solvers/test_scipy_solver.py
@@ -314,6 +314,30 @@ class TestScipySolver(unittest.TestCase):
             solution.y[0], 1 * np.exp(-0.1 * solution.t), decimal=5
         )
 
+    def test_model_solver_manually_update_initial_conditions(self):
+        # Create model
+        model = pybamm.BaseModel()
+        var1 = pybamm.Variable("var1")
+        model.rhs = {var1: -var1}
+        model.initial_conditions = {var1: 1}
+
+        # Solve
+        solver = pybamm.ScipySolver(rtol=1e-8, atol=1e-8)
+        t_eval = np.linspace(0, 5, 100)
+        solution = solver.solve(model, t_eval)
+        np.testing.assert_array_almost_equal(
+            solution.y[0], 1 * np.exp(-solution.t), decimal=5
+        )
+
+        # Change initial conditions and solve again
+        model.concatenated_initial_conditions = pybamm.NumpyConcatenation(
+            pybamm.Vector([[2]])
+        )
+        solution = solver.solve(model, t_eval)
+        np.testing.assert_array_almost_equal(
+            solution.y[0], 2 * np.exp(-solution.t), decimal=5
+        )
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")