diff --git a/python/wrap/solvers.cpp b/python/wrap/solvers.cpp
index 9e56444..e317802 100644
--- a/python/wrap/solvers.cpp
+++ b/python/wrap/solvers.cpp
@@ -1,142 +1,143 @@
 /**
  *  @file
  *  @section LICENSE
  *
  *  Copyright (©) 2016-19 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program is distributed in the hope that it will be useful,
  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  *  GNU Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #include "beck_teboulle.hh"
 #include "condat.hh"
 #include "contact_solver.hh"
 #include "ep_solver.hh"
 #include "epic.hh"
 #include "kato.hh"
 #include "kato_saturated.hh"
 #include "polonsky_keer_rey.hh"
 #include "polonsky_keer_tan.hh"
 #include "wrap.hh"
 /* -------------------------------------------------------------------------- */
 
 namespace tamaas {
 
 namespace wrap {
 
 /* -------------------------------------------------------------------------- */
 using namespace py::literals;
 
 class PyEPSolver : public EPSolver {
 public:
   using EPSolver::EPSolver;
   void solve() override { PYBIND11_OVERLOAD_PURE(void, EPSolver, solve); }
 
   void updateState() override {
     PYBIND11_OVERLOAD(void, EPSolver, updateState);
   }
 };
 
 /* -------------------------------------------------------------------------- */
 void wrapSolvers(py::module& mod) {
   py::class_<ContactSolver>(mod, "ContactSolver")
       // .def(py::init<Model&, const GridBase<Real>&, Real>(), "model"_a,
       //      "surface"_a, "tolerance"_a)
       .def("setMaxIterations", &ContactSolver::setMaxIterations, "max_iter"_a)
       .def("setDumpFrequency", &ContactSolver::setDumpFrequency, "dump_freq"_a)
       .def("addFunctionalTerm", &ContactSolver::addFunctionalTerm)
       .def("solve", py::overload_cast<std::vector<Real>>(&ContactSolver::solve),
            "target_force"_a)
       .def("solve", py::overload_cast<Real>(&ContactSolver::solve),
            "target_normal_pressure"_a);
 
   py::class_<PolonskyKeerRey, ContactSolver> pkr(mod, "PolonskyKeerRey");
   pkr.def(py::init<Model&, const GridBase<Real>&, Real, PolonskyKeerRey::type,
                    PolonskyKeerRey::type>(),
           "model"_a, "surface"_a, "tolerance"_a, "primal_type"_a,
           "constraint_type"_a, py::keep_alive<1, 2>(), py::keep_alive<1, 3>())
       .def("computeError", &PolonskyKeerRey::computeError);
 
   py::enum_<PolonskyKeerRey::type>(pkr, "type")
       .value("gap", PolonskyKeerRey::gap)
       .value("pressure", PolonskyKeerRey::pressure)
       .export_values();
 
   py::class_<KatoSaturated, PolonskyKeerRey>(mod, "KatoSaturated")
-      .def(py::init<Model&, const GridBase<Real>&, Real, Real>(),
-           "model"_a, "surface"_a, "tolerance"_a, "pmax"_a)
+      .def(py::init<Model&, const GridBase<Real>&, Real, Real>(), "model"_a,
+           "surface"_a, "tolerance"_a, "pmax"_a, py::keep_alive<1, 2>(),
+           py::keep_alive<1, 3>())
       .def_property("pmax", &KatoSaturated::getPMax, &KatoSaturated::setPMax);
 
   py::class_<Kato, ContactSolver> kato(mod, "Kato");
   kato.def(py::init<Model&, const GridBase<Real>&, Real, Real>(), "model"_a,
            "surface"_a, "tolerance"_a, "mu"_a, py::keep_alive<1, 2>(),
            py::keep_alive<1, 3>())
       .def("solve", &Kato::solve, "p0"_a, "proj_iter"_a = 50)
       .def("solveRelaxed", &Kato::solveRelaxed, "g0"_a)
       .def("solveRegularized", &Kato::solveRegularized, "p0"_a, "r"_a = 0.01)
       .def("computeCost", &Kato::computeCost, "use_tresca"_a = false);
 
   py::class_<BeckTeboulle, ContactSolver> bt(mod, "BeckTeboulle");
   bt.def(py::init<Model&, const GridBase<Real>&, Real, Real>(), "model"_a,
          "surface"_a, "tolerance"_a, "mu"_a, py::keep_alive<1, 2>(),
          py::keep_alive<1, 3>())
       .def("solve", &BeckTeboulle::solve, "p0"_a)
       .def("computeCost", &BeckTeboulle::computeCost);
 
   py::class_<Condat, ContactSolver> cd(mod, "Condat");
   cd.def(py::init<Model&, const GridBase<Real>&, Real, Real>(), "model"_a,
          "surface"_a, "tolerance"_a, "mu"_a, py::keep_alive<1, 2>(),
          py::keep_alive<1, 3>())
       .def("solve", &Condat::solve, "p0"_a, "grad_step"_a = 0.9)
       .def("computeCost", &Condat::computeCost);
 
   py::class_<PolonskyKeerTan, ContactSolver> pkt(mod, "PolonskyKeerTan");
   pkt.def(py::init<Model&, const GridBase<Real>&, Real, Real>(), "model"_a,
           "surface"_a, "tolerance"_a, "mu"_a, py::keep_alive<1, 2>(),
           py::keep_alive<1, 3>())
       .def("solve", &PolonskyKeerTan::solve, "p0"_a)
       .def("solveTresca", &PolonskyKeerTan::solveTresca, "p0"_a)
       .def("computeCost", &PolonskyKeerTan::computeCost,
            "use_tresca"_a = false);
 
   py::class_<EPSolver, PyEPSolver>(mod, "EPSolver")
-      .def(py::init<Residual&>(), "residual"_a)
+      .def(py::init<Residual&>(), "residual"_a, py::keep_alive<1, 2>())
       .def("solve", &EPSolver::solve)
       .def("getStrainIncrement", &EPSolver::getStrainIncrement,
            py::return_value_policy::reference_internal)
       .def("getResidual", &EPSolver::getResidual,
            py::return_value_policy::reference_internal)
       .def("updateState", &EPSolver::updateState);
 
   py::class_<EPICSolver>(mod, "EPICSolver")
       .def(py::init<ContactSolver&, EPSolver&, Real, Real>(),
            "contact_solver"_a, "elasto_plastic_solver"_a, "tolerance"_a = 1e-10,
-           "relaxation"_a = 0.3)
+           "relaxation"_a = 0.3, py::keep_alive<1, 2>(), py::keep_alive<1, 3>())
       .def("solve", &EPICSolver::solve, "force"_a)
       .def("solve",
            [](EPICSolver& solver, Real pressure) {
              return solver.solve(std::vector<Real>{pressure});
            },
            "normal_pressure"_a)
       .def("acceleratedSolve",
            [](EPICSolver& solver, Real pressure) {
              return solver.acceleratedSolve(std::vector<Real>{pressure});
            },
            "normal_pressure"_a);
 }
 
 }  // namespace wrap
 
 }  // namespace tamaas