diff --git a/xbanalysis/__init__.py b/xbanalysis/__init__.py
index 57b2abe285bb2a77cbf467c0ae421d5d00bed2db..c2e384f8f77911a5682fb8d84ad8650ac5395175 100644
--- a/xbanalysis/__init__.py
+++ b/xbanalysis/__init__.py
@@ -7,6 +7,7 @@ Package generating XBA problems from SBML.
 
 from .model import *
 from .problems import *
+from .solvers import *
 from .utils import *
 
 from . import model, problems, utils
@@ -15,3 +16,4 @@ __all__ = []
 __all__ += model.__all__
 __all__ += problems.__all__
 __all__ += utils.__all__
+__all__ += solvers.__all__
diff --git a/xbanalysis/_version.py b/xbanalysis/_version.py
index 2a1d8eaf56f1a62f53d2d3165eeb9157eef640b6..120eb33ded3939e49674a86660c986ede1a23ee9 100644
--- a/xbanalysis/_version.py
+++ b/xbanalysis/_version.py
@@ -1,5 +1,5 @@
 """
 Definition of version string.
 """
-__version__ = "0.2.0"
+__version__ = "0.3.0"
 program_name = 'xbanalysis'
diff --git a/xbanalysis/problems/gba_problem.py b/xbanalysis/problems/gba_problem.py
index 3f73fc705276a3400d3b8fe824cb2eabb02e67a0..eaeb62b36f17d1342f30a1032cfd481b06685337 100644
--- a/xbanalysis/problems/gba_problem.py
+++ b/xbanalysis/problems/gba_problem.py
@@ -15,10 +15,10 @@ from xbanalysis.utils.hessians import (hess_from_grad_grad, hess_from_v_hesss, h
 
 # Decorator to cache last result
 def cache_last(func):
-    def wrapper(self, x):
+    def wrapper(self, *args):
         fname = wrapper.__name__
-        if fname not in GbaProblem.last_values or np.array_equal(GbaProblem.last_values[fname][0], x) is False:
-            GbaProblem.last_values[fname] = [x, func(self, x)]
+        if fname not in GbaProblem.last_values or np.array_equal(GbaProblem.last_values[fname][0], args[0]) is False:
+            GbaProblem.last_values[fname] = [args[0], func(self, args[0])]
         return GbaProblem.last_values[fname][1]
     wrapper.__name__ = func.__name__
     return wrapper
@@ -298,7 +298,7 @@ class GbaProblem:
         enz_trans_hesss = types.FunctionType(hesss_code.co_consts[0], self.model_params)
         return enz_trans, enz_trans_jac, enz_trans_hesss
 
-#    @cache_last
+    @cache_last
     def get_growth_rate_0(self, x):
         """Calculate growth rate at x, assuming no protein degradation
 
@@ -323,6 +323,7 @@ class GbaProblem:
         #    self.last_gr0 = gr0
         return gr0
 
+    @cache_last
     def get_growth_rate_0_grad(self, x):
         """Calculate growth rate gradient at x, assuming no protein degradation
 
@@ -345,6 +346,7 @@ class GbaProblem:
         gr0_grad = - gr0 ** 2 * tau_grad
         return np.nan_to_num(gr0_grad, nan=0.0)
 
+    @cache_last
     def get_growth_rate_0_hess(self, x):
         """Calculate growth rate Hessian at x, assuming no protein degradation
 
@@ -389,6 +391,7 @@ class GbaProblem:
         gr0_hess = - 2 * gr0 * hess_a - gr0 ** 2 * (hess_b1 + hess_b2 + hess_b3)
         return np.nan_to_num(gr0_hess, nan=0.0)
 
+    @cache_last
     def get_growth_rate(self, x):
         """Calculate growth rate at x, considering protein degradation
 
@@ -413,6 +416,7 @@ class GbaProblem:
             gr = gr0 * (1.0 + enz_to_mras.dot(enz_to_pstmas))
         return gr
 
+    @cache_last
     def get_growth_rate_grad(self, x):
         """Calculate growth rate gradient at x, considering protein degradation
 
@@ -440,6 +444,7 @@ class GbaProblem:
                        + gr0 * (enz_to_mras.dot(enz_to_pstmas_jac) + enz_to_pstmas.dot(enz_to_mras_jac)))
         return np.nan_to_num(gr_grad, nan=0.0)
 
+    @cache_last
     def get_growth_rate_hess(self, x):
         """Calculate growth rate Hessian at x, considering protein degradation
 
@@ -490,6 +495,7 @@ class GbaProblem:
             gr_hess = hess_a + 2 * hess_b + gr0 * (hess_c1 + hess_c2 + hess_c3 + hess_c4)
         return np.nan_to_num(gr_hess, nan=0.0)
 
+    @cache_last
     def get_alphas(self, x):
         """Calculate ribosomal allocations at x, considering protein degradation
 
@@ -514,6 +520,7 @@ class GbaProblem:
             alphas = alphas0 - enz_to_mras * enz_to_pstmas
         return alphas
 
+    @cache_last
     def heq_mass_balance(self, x):
         """Calculate mass balance equations at x for metabolites
 
@@ -548,6 +555,7 @@ class GbaProblem:
             metab_mb = metabolic_mb + gr * (ps_drain - cm) + synth_enz_mb
         return metab_mb
 
+    @cache_last
     def heq_mass_balance_jac(self, x):
         """Calculate Jacobian of mass balance equations at x for metabolites
 
@@ -601,6 +609,7 @@ class GbaProblem:
             metab_mb_jac = metabolic_jac + dilute_1_jac + dilute_2_jac + synth_enz_mb_jac
         return np.nan_to_num(metab_mb_jac, nan=0.0)
 
+    @cache_last
     def heq_mass_balance_hess(self, x):
         """Calculate Hessians of mass balance equations at x
 
@@ -681,15 +690,19 @@ class GbaProblem:
             metab_mb_hesss = metabolic_hesss + dilute_1_hesss + 2 * dilute_2_hesss + synth_enz_mb_hesss
         return np.nan_to_num(metab_mb_hesss, nan=0.0)
 
+    @cache_last
     def heq_enz_trans(self, x):
         return self.enz_trans(x)
 
+    @cache_last
     def heq_enz_trans_jac(self, x):
         return self.enz_trans_jac(x)
 
+    @cache_last
     def heq_enz_trans_hess(self, x):
         return self.enz_trans_hesss(x)
 
+    @cache_last
     def heq_density(self, x):
         """Calculate density (g/l) constraint at x
 
@@ -712,10 +725,10 @@ class GbaProblem:
             density = self.model_params['mws'][self.var_mask_enz].dot(ce)
         else:
             density = self.model_params['mws'].dot(x)
-        # !!! when using Ipopt, the rho could be specified as a bound!
-        return np.array(self.model_params['rho'] - density)
+        # TODO check if rho should be specified as a bound!
+        return np.array([self.model_params['rho'] - density])
 
-    # noinspection PyUnusedLocal
+    @cache_last
     def heq_density_jac(self, x):
         """Calculate gradient of density constraint at x
 
@@ -732,13 +745,14 @@ class GbaProblem:
         return density_grad
 
     # noinspection PyUnusedLocal
+    @cache_last
     def heq_density_hess(self, x):
         """Calculate hessian of density constraint at x
 
         :param x: optimization variables
         :type x: numpy.ndarray
         :returns: mass balance Hessian at x
-        :rtype: numpy.ndarray (2D), shape [x,x]
+        :rtype: numpy.ndarray (3D), shape [1,x,x]
         """
-        hess = np.zeros((self.n_vars, self.n_vars))
+        hess = np.zeros((1, self.n_vars, self.n_vars))
         return np.tril(hess)
diff --git a/xbanalysis/solvers/__init__.py b/xbanalysis/solvers/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..37b8472d758df8fde32209a4a5daa2c4def29a96
--- /dev/null
+++ b/xbanalysis/solvers/__init__.py
@@ -0,0 +1,5 @@
+"""Subpackage with XBA model classes """
+
+from .gba_ipopt_problem import GbaIpoptProblem
+
+__all__ = ['GbaIpoptProblem']
diff --git a/xbanalysis/solvers/gba_ipopt_problem.py b/xbanalysis/solvers/gba_ipopt_problem.py
new file mode 100644
index 0000000000000000000000000000000000000000..480c30984387b401a5939051a05c57c4f7efe3c1
--- /dev/null
+++ b/xbanalysis/solvers/gba_ipopt_problem.py
@@ -0,0 +1,70 @@
+"""Implementation of GbaIpoptProblem class.
+
+Peter Schubert, HHU Duesseldorf, June 2022
+"""
+
+import numpy as np
+
+
+# based on cyipopt examples
+class GbaIpoptProblem:
+
+    def __init__(self, gba_problem):
+        """Initialize GbaIpoptProblem.
+
+        :param gba_problem: GbaProblem from where to extract relevant parameters
+        :type gba_problem: GbaProblem
+        """
+        self.gba_problem = gba_problem
+        self._n_vars = len(gba_problem.var_initial)
+        self.report_freq = 0
+
+        self.nfev = 0
+        self.njev = 0
+        self.nit = 0
+        tmp_x = np.ones(self._n_vars) * .01
+        self.constraint_dims = np.array([len(gba_problem.heq_mass_balance(tmp_x)),
+                                         len(gba_problem.heq_enz_trans(tmp_x)),
+                                         len(gba_problem.heq_density(tmp_x))])
+
+    def objective(self, x):
+        self.nfev += 1
+        return -self.gba_problem.get_growth_rate(x)
+
+    def gradient(self, x):
+        self.njev += 1
+        return -self.gba_problem.get_growth_rate_grad(x)
+
+    def constraints(self, x):
+        return np.concatenate((self.gba_problem.heq_mass_balance(x),
+                               self.gba_problem.heq_enz_trans(x),
+                               self.gba_problem.heq_density(x)))
+
+    def jacobian(self, x):
+        return (np.vstack((self.gba_problem.heq_mass_balance_jac(x),
+                           self.gba_problem.heq_enz_trans_jac(x),
+                           self.gba_problem.heq_density_jac(x)))
+                ).flatten()
+
+    def hessianstructure(self):
+        return np.tril_indices(self._n_vars)
+
+    def hessian(self, x, lagrange, obj_factor):
+        hess = obj_factor * -self.gba_problem.get_growth_rate_hess(x)
+        hesss = np.vstack((self.gba_problem.heq_mass_balance_hess(x),
+                           self.gba_problem.heq_enz_trans_hess(x),
+                           self.gba_problem.heq_density_hess(x)))
+        hess += (lagrange.reshape((-1, 1, 1)) * hesss).sum(axis=0)
+        return hess[np.tril_indices(self._n_vars)]
+
+    def set_report_freq(self, mod_val=0):
+        self.report_freq = mod_val
+
+    # noinspection PyUnusedLocal
+    def intermediate(self, alg_mod, iter_count, obj_value, inf_pr, inf_du, mu,
+                     d_norm, regularization_size, alpha_du, alpha_pr,
+                     ls_trials):
+        self.nit = iter_count
+        if self.report_freq > 0:
+            if iter_count % self.report_freq == 0:
+                print(f'[{iter_count:5d}] growth rate: {-3600 *obj_value:.5f}')
diff --git a/xbanalysis/utils/opt_results.py b/xbanalysis/utils/opt_results.py
index df73f038418b072c432de7ccba81f9df0d8d682a..ec6b9afce3d641274ffe4a67f39aedaf8352a28e 100644
--- a/xbanalysis/utils/opt_results.py
+++ b/xbanalysis/utils/opt_results.py
@@ -54,7 +54,7 @@ class OptResults:
             self.ress['rho_m'][i] = rho[self.problem.var_mask_metab]
             self.ress['rho_e'][i] = rho[self.problem.var_mask_enz]
 
-            # !!! what volume to consider, possibly the volume of the ribosome ??
+            # TODO extract reaction volume from kinetics to scale fluxes here
             scale_factor = 3600 * 1e3 / self.problem.model_params['rho']
             self.ress['mr_fluxes'][i] = scale_factor * self.problem.mras(cx) / self.problem.var_vols[-1]
             self.ress['psm_fluxes'][i] = scale_factor / (self.problem.pstm_scale * self.problem.pstmas(cx))