diff --git a/ql/termstructures/yield/fittedbonddiscountcurve.cpp b/ql/termstructures/yield/fittedbonddiscountcurve.cpp
index d7f466137a5..2d0aa4b279f 100644
--- a/ql/termstructures/yield/fittedbonddiscountcurve.cpp
+++ b/ql/termstructures/yield/fittedbonddiscountcurve.cpp
@@ -130,10 +130,15 @@ namespace QuantLib {
         ext::shared_ptr<OptimizationMethod> optimizationMethod,
         Array l2,
         const Real minCutoffTime,
-        const Real maxCutoffTime)
+        const Real maxCutoffTime,
+        Constraint constraint)
     : constrainAtZero_(constrainAtZero), weights_(weights), l2_(std::move(l2)),
       calculateWeights_(weights.empty()), optimizationMethod_(std::move(optimizationMethod)),
-      minCutoffTime_(minCutoffTime), maxCutoffTime_(maxCutoffTime) {}
+      constraint_(std::move(constraint)),
+      minCutoffTime_(minCutoffTime), maxCutoffTime_(maxCutoffTime) {
+        if (constraint_.empty())
+            constraint_ = NoConstraint();
+    }
 
     void FittedBondDiscountCurve::FittingMethod::init() {
         // yield conventions
@@ -188,7 +193,6 @@ namespace QuantLib {
     void FittedBondDiscountCurve::FittingMethod::calculate() {
 
         FittingCost& costFunction = *costFunction_;
-        Constraint constraint = NoConstraint();
 
         // start with the guess solution, if it exists
         Array x(size(), 0.0);
@@ -221,7 +225,7 @@ namespace QuantLib {
         if(!optimization){
             optimization = ext::make_shared<Simplex>(curve_->simplexLambda_);
         }
-        Problem problem(costFunction, constraint, x);
+        Problem problem(costFunction, constraint_, x);
 
         Real rootEpsilon = curve_->accuracy_;
         Real functionEpsilon =  curve_->accuracy_;
diff --git a/ql/termstructures/yield/fittedbonddiscountcurve.hpp b/ql/termstructures/yield/fittedbonddiscountcurve.hpp
index dd601d0868a..6db9eba8d43 100644
--- a/ql/termstructures/yield/fittedbonddiscountcurve.hpp
+++ b/ql/termstructures/yield/fittedbonddiscountcurve.hpp
@@ -31,6 +31,7 @@
 #include <ql/patterns/lazyobject.hpp>
 #include <ql/math/array.hpp>
 #include <ql/utilities/clone.hpp>
+#include <ql/math/optimization/constraint.hpp>
 
 namespace QuantLib {
 
@@ -209,6 +210,8 @@ namespace QuantLib {
         Array l2() const;
         //! return optimization method being used
         ext::shared_ptr<OptimizationMethod> optimizationMethod() const;
+        //! return optimization contraint
+        const Constraint& constraint() const;
         //! open discountFunction to public
         DiscountFactor discount(const Array& x, Time t) const;
       protected:
@@ -219,7 +222,8 @@ namespace QuantLib {
                           ext::shared_ptr<OptimizationMethod>(),
                       Array l2 = Array(),
                       Real minCutoffTime = 0.0,
-                      Real maxCutoffTime = QL_MAX_REAL);
+                      Real maxCutoffTime = QL_MAX_REAL,
+                      Constraint constraint = NoConstraint());
         //! rerun every time instruments/referenceDate changes
         virtual void init();
         //! discount function called by FittedBondDiscountCurve
@@ -257,6 +261,8 @@ namespace QuantLib {
         EndCriteria::Type errorCode_ = EndCriteria::None;
         // optimization method to be used, if none provided use Simplex
         ext::shared_ptr<OptimizationMethod> optimizationMethod_;
+        // optimization constraint, if none provided use NoConstraint
+        Constraint constraint_;
         // flat extrapolation of instantaneous forward before / after cutoff
         Real minCutoffTime_, maxCutoffTime_;
     };
@@ -329,6 +335,10 @@ namespace QuantLib {
         return optimizationMethod_;
     }
 
+    inline const Constraint& FittedBondDiscountCurve::FittingMethod::constraint() const {
+        return constraint_;
+    }
+
     inline DiscountFactor FittedBondDiscountCurve::FittingMethod::discount(const Array& x, Time t) const {
         if (t < minCutoffTime_) {
             // flat fwd extrapolation before min cutoff time
diff --git a/ql/termstructures/yield/nonlinearfittingmethods.cpp b/ql/termstructures/yield/nonlinearfittingmethods.cpp
index f9f1b217cc9..dcaebe093c6 100644
--- a/ql/termstructures/yield/nonlinearfittingmethods.cpp
+++ b/ql/termstructures/yield/nonlinearfittingmethods.cpp
@@ -33,39 +33,37 @@ namespace QuantLib {
         const Real minCutoffTime,
         const Real maxCutoffTime,
         const Size numCoeffs,
-        const Real fixedKappa)
-    : FittedBondDiscountCurve::FittingMethod(
-          constrainAtZero, weights, optimizationMethod, l2, minCutoffTime, maxCutoffTime),
-          numCoeffs_(numCoeffs), fixedKappa_(fixedKappa) 
-    {
+        const Real fixedKappa,
+        Constraint constraint)
+    : FittedBondDiscountCurve::FittingMethod(constrainAtZero, weights, optimizationMethod, l2,
+                                             minCutoffTime, maxCutoffTime, std::move(constraint)),
+      numCoeffs_(numCoeffs), fixedKappa_(fixedKappa) {
         QL_REQUIRE(ExponentialSplinesFitting::size() > 0, "At least 1 unconstrained coefficient required");
     }
 
-    ExponentialSplinesFitting::ExponentialSplinesFitting(bool constrainAtZero,
+    ExponentialSplinesFitting::ExponentialSplinesFitting(
+        bool constrainAtZero,
         const Array& weights,
         const Array& l2, const Real minCutoffTime, const Real maxCutoffTime,
-        const Size numCoeffs, const Real fixedKappa)
-        : FittedBondDiscountCurve::FittingMethod(constrainAtZero, weights, ext::shared_ptr<OptimizationMethod>(), l2,
-                                                 minCutoffTime, maxCutoffTime),
-          numCoeffs_(numCoeffs),fixedKappa_(fixedKappa) 
-    {
-        QL_REQUIRE(ExponentialSplinesFitting::size() > 0, "At least 1 unconstrained coefficient required");
-    }
+        const Size numCoeffs, const Real fixedKappa,
+        Constraint constraint)
+    : ExponentialSplinesFitting(constrainAtZero, weights, {}, l2,
+                                minCutoffTime, maxCutoffTime,
+                                numCoeffs, fixedKappa, std::move(constraint)) {}
 
-    ExponentialSplinesFitting::ExponentialSplinesFitting(bool constrainAtZero,
-                              const Size numCoeffs,
-                              const Real fixedKappa,
-                              const Array& weights )
-    : FittedBondDiscountCurve::FittingMethod(constrainAtZero, weights, ext::shared_ptr<OptimizationMethod>(), Array(),0.0,QL_MAX_REAL),
-          numCoeffs_(numCoeffs), fixedKappa_(fixedKappa)
-    {
-        QL_REQUIRE(ExponentialSplinesFitting::size() > 0, "At least 1 unconstrained coefficient required");
-    }
+    ExponentialSplinesFitting::ExponentialSplinesFitting(
+        bool constrainAtZero,
+        const Size numCoeffs,
+        const Real fixedKappa,
+        const Array& weights,
+        Constraint constraint)
+    : ExponentialSplinesFitting(constrainAtZero, weights, {}, Array(),
+                                0.0, QL_MAX_REAL,
+                                numCoeffs, fixedKappa, std::move(constraint)) {}
 
     std::unique_ptr<FittedBondDiscountCurve::FittingMethod>
     ExponentialSplinesFitting::clone() const {
-        return std::unique_ptr<FittedBondDiscountCurve::FittingMethod>(
-                                        new ExponentialSplinesFitting(*this)); 
+        return std::make_unique<ExponentialSplinesFitting>(*this); 
     }
 
     Size ExponentialSplinesFitting::size() const {
@@ -107,20 +105,23 @@ namespace QuantLib {
         const ext::shared_ptr<OptimizationMethod>& optimizationMethod,
         const Array& l2,
         const Real minCutoffTime,
-        const Real maxCutoffTime)
-    : FittedBondDiscountCurve::FittingMethod(
-          true, weights, optimizationMethod, l2, minCutoffTime, maxCutoffTime) {}
+        const Real maxCutoffTime,
+        Constraint constraint)
+    : FittedBondDiscountCurve::FittingMethod(true, weights, optimizationMethod, l2,
+                                             minCutoffTime, maxCutoffTime, std::move(constraint)) {}
 
-    NelsonSiegelFitting::NelsonSiegelFitting(const Array& weights,
-                                             const Array& l2,
-                                             const Real minCutoffTime, const Real maxCutoffTime)
-        : FittedBondDiscountCurve::FittingMethod(true, weights, ext::shared_ptr<OptimizationMethod>(), l2,
-                                                 minCutoffTime, maxCutoffTime) {}
+    NelsonSiegelFitting::NelsonSiegelFitting(
+        const Array& weights,
+        const Array& l2,
+        const Real minCutoffTime,
+        const Real maxCutoffTime,
+        Constraint constraint)
+    : NelsonSiegelFitting(weights, {}, l2,
+                          minCutoffTime, maxCutoffTime, std::move(constraint)) {}
 
     std::unique_ptr<FittedBondDiscountCurve::FittingMethod>
     NelsonSiegelFitting::clone() const {
-        return std::unique_ptr<FittedBondDiscountCurve::FittingMethod>(
-                                              new NelsonSiegelFitting(*this));
+        return std::make_unique<NelsonSiegelFitting>(*this);
     }
 
     Size NelsonSiegelFitting::size() const {
@@ -143,19 +144,22 @@ namespace QuantLib {
                                      const ext::shared_ptr<OptimizationMethod>& optimizationMethod,
                                      const Array& l2,
                                      const Real minCutoffTime,
-                                     const Real maxCutoffTime)
-    : FittedBondDiscountCurve::FittingMethod(
-          true, weights, optimizationMethod, l2, minCutoffTime, maxCutoffTime) {}
+                                     const Real maxCutoffTime,
+                                     Constraint constraint)
+    : FittedBondDiscountCurve::FittingMethod(true, weights, optimizationMethod, l2,
+                                             minCutoffTime, maxCutoffTime, std::move(constraint)) {}
 
     SvenssonFitting::SvenssonFitting(const Array& weights,
-        const Array& l2, const Real minCutoffTime, const Real maxCutoffTime)
-        : FittedBondDiscountCurve::FittingMethod(true, weights, ext::shared_ptr<OptimizationMethod>(), l2,
-                                                 minCutoffTime, maxCutoffTime) {}
+                                     const Array& l2,
+                                     const Real minCutoffTime,
+                                     const Real maxCutoffTime,
+                                     Constraint constraint)
+    : SvenssonFitting(weights, {}, l2,
+                      minCutoffTime, maxCutoffTime, std::move(constraint)) {}
 
     std::unique_ptr<FittedBondDiscountCurve::FittingMethod>
     SvenssonFitting::clone() const {
-        return std::unique_ptr<FittedBondDiscountCurve::FittingMethod>(
-                                              new SvenssonFitting(*this));
+        return std::make_unique<SvenssonFitting>(*this);
     }
 
     Size SvenssonFitting::size() const {
@@ -184,9 +188,10 @@ namespace QuantLib {
         const ext::shared_ptr<OptimizationMethod>& optimizationMethod,
         const Array& l2,
         const Real minCutoffTime,
-        const Real maxCutoffTime)
-    : FittedBondDiscountCurve::FittingMethod(
-          constrainAtZero, weights, optimizationMethod, l2, minCutoffTime, maxCutoffTime),
+        const Real maxCutoffTime,
+        Constraint constraint)
+    : FittedBondDiscountCurve::FittingMethod(constrainAtZero, weights, optimizationMethod, l2,
+                                             minCutoffTime, maxCutoffTime, std::move(constraint)),
       splines_(3, knots.size() - 5, knots) {
 
         QL_REQUIRE(knots.size() >= 8,
@@ -208,34 +213,16 @@ namespace QuantLib {
         }
     }
 
-    CubicBSplinesFitting::CubicBSplinesFitting(const std::vector<Time>& knots,
+    CubicBSplinesFitting::CubicBSplinesFitting(
+        const std::vector<Time>& knots,
         bool constrainAtZero,
         const Array& weights,
         const Array& l2,
-        const Real minCutoffTime, const Real maxCutoffTime)
-        : FittedBondDiscountCurve::FittingMethod(constrainAtZero, weights, ext::shared_ptr<OptimizationMethod>(), l2,
-                                                 minCutoffTime, maxCutoffTime),
-        splines_(3, knots.size() - 5, knots) {
-
-        QL_REQUIRE(knots.size() >= 8,
-            "At least 8 knots are required");
-        Size basisFunctions = knots.size() - 4;
-
-        if (constrainAtZero) {
-            size_ = basisFunctions - 1;
-
-            // Note: A small but nonzero N_th basis function at t=0 may
-            // lead to an ill conditioned problem
-            N_ = 1;
-
-            QL_REQUIRE(std::abs(splines_(N_, 0.0)) > QL_EPSILON,
-                "N_th cubic B-spline must be nonzero at t=0");
-        }
-        else {
-            size_ = basisFunctions;
-            N_ = 0;
-        }
-    }
+        const Real minCutoffTime,
+        const Real maxCutoffTime,
+        Constraint constraint)
+    : CubicBSplinesFitting(knots, constrainAtZero, weights, {}, l2,
+                           minCutoffTime, maxCutoffTime, std::move(constraint)) {}
 
     Real CubicBSplinesFitting::basisFunction(Integer i, Time t) const {
         return splines_(i,t);
@@ -243,8 +230,7 @@ namespace QuantLib {
 
     std::unique_ptr<FittedBondDiscountCurve::FittingMethod>
     CubicBSplinesFitting::clone() const {
-        return std::unique_ptr<FittedBondDiscountCurve::FittingMethod>(
-                                             new CubicBSplinesFitting(*this));
+        return std::make_unique<CubicBSplinesFitting>(*this);
     }
 
     Size CubicBSplinesFitting::size() const {
@@ -287,22 +273,26 @@ namespace QuantLib {
         const ext::shared_ptr<OptimizationMethod>& optimizationMethod,
         const Array& l2,
         const Real minCutoffTime,
-        const Real maxCutoffTime)
-    : FittedBondDiscountCurve::FittingMethod(
-          constrainAtZero, weights, optimizationMethod, l2, minCutoffTime, maxCutoffTime),
+        const Real maxCutoffTime,
+        Constraint constraint)
+    : FittedBondDiscountCurve::FittingMethod(constrainAtZero, weights, optimizationMethod, l2,
+                                             minCutoffTime, maxCutoffTime, std::move(constraint)),
       size_(constrainAtZero ? degree : degree + 1) {}
 
-    SimplePolynomialFitting::SimplePolynomialFitting(Natural degree, bool constrainAtZero,
-                                                     const Array& weights, const Array& l2,
-                                                     const Real minCutoffTime, const Real maxCutoffTime)
-        : FittedBondDiscountCurve::FittingMethod(constrainAtZero, weights, 
-                                                 ext::shared_ptr<OptimizationMethod>(), l2, minCutoffTime, maxCutoffTime),
-        size_(constrainAtZero ? degree : degree + 1) {}
+    SimplePolynomialFitting::SimplePolynomialFitting(
+        Natural degree,
+        bool constrainAtZero,
+        const Array& weights,
+        const Array& l2,
+        const Real minCutoffTime,
+        const Real maxCutoffTime,
+        Constraint constraint)
+    : SimplePolynomialFitting(degree, constrainAtZero, weights, {}, l2,
+                              minCutoffTime, maxCutoffTime, std::move(constraint)) {}
 
     std::unique_ptr<FittedBondDiscountCurve::FittingMethod>
     SimplePolynomialFitting::clone() const {
-        return std::unique_ptr<FittedBondDiscountCurve::FittingMethod>(
-                                          new SimplePolynomialFitting(*this));
+        return std::make_unique<SimplePolynomialFitting>(*this);
     }
 
     Size SimplePolynomialFitting::size() const {
@@ -342,8 +332,7 @@ namespace QuantLib {
 
     std::unique_ptr<FittedBondDiscountCurve::FittingMethod>
     SpreadFittingMethod::clone() const {
-        return std::unique_ptr<FittedBondDiscountCurve::FittingMethod>(
-                                          new SpreadFittingMethod(*this));
+        return std::make_unique<SpreadFittingMethod>(*this);
     }
 
     Size SpreadFittingMethod::size() const {
diff --git a/ql/termstructures/yield/nonlinearfittingmethods.hpp b/ql/termstructures/yield/nonlinearfittingmethods.hpp
index e03a765ed93..ce4eca342c4 100644
--- a/ql/termstructures/yield/nonlinearfittingmethods.hpp
+++ b/ql/termstructures/yield/nonlinearfittingmethods.hpp
@@ -51,24 +51,26 @@ namespace QuantLib {
       public:
         ExponentialSplinesFitting(bool constrainAtZero = true,
                                   const Array& weights = Array(),
-                                  const ext::shared_ptr<OptimizationMethod>& optimizationMethod =
-                                      ext::shared_ptr<OptimizationMethod>(),
+                                  const ext::shared_ptr<OptimizationMethod>& optimizationMethod = {},
                                   const Array& l2 = Array(),
                                   Real minCutoffTime = 0.0,
                                   Real maxCutoffTime = QL_MAX_REAL,
                                   Size numCoeffs = 9,
-                                  Real fixedKappa = Null<Real>());
+                                  Real fixedKappa = Null<Real>(),
+                                  Constraint constraint = NoConstraint());
         ExponentialSplinesFitting(bool constrainAtZero,
                                   const Array& weights,
                                   const Array& l2,
                                   Real minCutoffTime = 0.0,
                                   Real maxCutoffTime = QL_MAX_REAL,
                                   Size numCoeffs = 9,
-                                  Real fixedKappa = Null<Real>());
+                                  Real fixedKappa = Null<Real>(),
+                                  Constraint constraint = NoConstraint());
         ExponentialSplinesFitting(bool constrainAtZero, 
                                   Size numCoeffs, 
                                   Real fixedKappa, 
-                                  const Array& weights = Array() );
+                                  const Array& weights = Array(),
+                                  Constraint constraint = NoConstraint());
 
 
         std::unique_ptr<FittedBondDiscountCurve::FittingMethod> clone() const override;
@@ -94,15 +96,16 @@ namespace QuantLib {
         : public FittedBondDiscountCurve::FittingMethod {
       public:
         NelsonSiegelFitting(const Array& weights = Array(),
-                            const ext::shared_ptr<OptimizationMethod>& optimizationMethod =
-                                ext::shared_ptr<OptimizationMethod>(),
+                            const ext::shared_ptr<OptimizationMethod>& optimizationMethod = {},
                             const Array& l2 = Array(),
                             Real minCutoffTime = 0.0,
-                            Real maxCutoffTime = QL_MAX_REAL);
+                            Real maxCutoffTime = QL_MAX_REAL,
+                            Constraint constraint = NoConstraint());
         NelsonSiegelFitting(const Array& weights,
                             const Array& l2,
                             Real minCutoffTime = 0.0,
-                            Real maxCutoffTime = QL_MAX_REAL);
+                            Real maxCutoffTime = QL_MAX_REAL,
+                            Constraint constraint = NoConstraint());
         std::unique_ptr<FittedBondDiscountCurve::FittingMethod> clone() const override;
       private:
         Size size() const override;
@@ -126,15 +129,16 @@ namespace QuantLib {
         : public FittedBondDiscountCurve::FittingMethod {
       public:
         SvenssonFitting(const Array& weights = Array(),
-                        const ext::shared_ptr<OptimizationMethod>& optimizationMethod =
-                            ext::shared_ptr<OptimizationMethod>(),
+                        const ext::shared_ptr<OptimizationMethod>& optimizationMethod = {},
                         const Array& l2 = Array(),
                         Real minCutoffTime = 0.0,
-                        Real maxCutoffTime = QL_MAX_REAL);
+                        Real maxCutoffTime = QL_MAX_REAL,
+                        Constraint constraint = NoConstraint());
         SvenssonFitting(const Array& weights,
                         const Array& l2,
                         Real minCutoffTime = 0.0,
-                        Real maxCutoffTime = QL_MAX_REAL);
+                        Real maxCutoffTime = QL_MAX_REAL,
+                        Constraint constraint = NoConstraint());
         std::unique_ptr<FittedBondDiscountCurve::FittingMethod> clone() const override;
       private:
         Size size() const override;
@@ -167,17 +171,18 @@ namespace QuantLib {
         CubicBSplinesFitting(const std::vector<Time>& knotVector,
                              bool constrainAtZero = true,
                              const Array& weights = Array(),
-                             const ext::shared_ptr<OptimizationMethod>& optimizationMethod =
-                                 ext::shared_ptr<OptimizationMethod>(),
+                             const ext::shared_ptr<OptimizationMethod>& optimizationMethod = {},
                              const Array& l2 = Array(),
                              Real minCutoffTime = 0.0,
-                             Real maxCutoffTime = QL_MAX_REAL);
+                             Real maxCutoffTime = QL_MAX_REAL,
+                             Constraint constraint = NoConstraint());
         CubicBSplinesFitting(const std::vector<Time>& knotVector,
                              bool constrainAtZero,
                              const Array& weights,
                              const Array& l2,
                              Real minCutoffTime = 0.0,
-                             Real maxCutoffTime = QL_MAX_REAL);
+                             Real maxCutoffTime = QL_MAX_REAL,
+                             Constraint constraint = NoConstraint());
         //! cubic B-spline basis functions
         Real basisFunction(Integer i, Time t) const;
         std::unique_ptr<FittedBondDiscountCurve::FittingMethod> clone() const override;
@@ -207,17 +212,18 @@ namespace QuantLib {
         SimplePolynomialFitting(Natural degree,
                                 bool constrainAtZero = true,
                                 const Array& weights = Array(),
-                                const ext::shared_ptr<OptimizationMethod>& optimizationMethod =
-                                    ext::shared_ptr<OptimizationMethod>(),
+                                const ext::shared_ptr<OptimizationMethod>& optimizationMethod = {},
                                 const Array& l2 = Array(),
                                 Real minCutoffTime = 0.0,
-                                Real maxCutoffTime = QL_MAX_REAL);
+                                Real maxCutoffTime = QL_MAX_REAL,
+                                Constraint constraint = NoConstraint());
         SimplePolynomialFitting(Natural degree,
                                 bool constrainAtZero,
                                 const Array& weights,
                                 const Array& l2,
                                 Real minCutoffTime = 0.0,
-                                Real maxCutoffTime = QL_MAX_REAL);
+                                Real maxCutoffTime = QL_MAX_REAL,
+                                Constraint constraint = NoConstraint());
         std::unique_ptr<FittedBondDiscountCurve::FittingMethod> clone() const override;
       private:
         Size size() const override;
diff --git a/test-suite/fittedbonddiscountcurve.cpp b/test-suite/fittedbonddiscountcurve.cpp
index 61a20a98cf8..fd89877bb96 100644
--- a/test-suite/fittedbonddiscountcurve.cpp
+++ b/test-suite/fittedbonddiscountcurve.cpp
@@ -272,6 +272,61 @@ BOOST_AUTO_TEST_CASE(testGuessSize) {
                           ExpectedErrorMessage("wrong size for guess"));
 }
 
+
+
+
+BOOST_AUTO_TEST_CASE(testConstraint) {
+
+    BOOST_TEST_MESSAGE("Testing that fitted bond curves respect passed constraint...");
+
+    class FlatZero : public FittedBondDiscountCurve::FittingMethod {
+      public:
+        FlatZero(Constraint constraint = NoConstraint())
+        : FittedBondDiscountCurve::FittingMethod(true,        // constrainAtZero
+                                                 Array(),     // weights
+                                                 ext::shared_ptr<OptimizationMethod>(),
+                                                 Array(),     // l2
+                                                 0.0,         // minCutoffTime
+                                                 QL_MAX_REAL, //maxCutoffTime
+                                                 std::move(constraint)) {}
+
+        std::unique_ptr<FittedBondDiscountCurve::FittingMethod> clone() const override {
+            return std::make_unique<FlatZero>(*this);
+        }
+
+      private:
+        Size size() const override { return 1; }
+
+        DiscountFactor discountFunction(const Array& x, Time t) const override {
+            Real zeroRate = x[0];
+            DiscountFactor d = std::exp(-zeroRate * t);
+            return d;
+        }
+    };
+
+    Date today = Settings::instance().evaluationDate();
+    auto bond1 = ext::make_shared<ZeroCouponBond>(3, TARGET(), 100.0, today + Period(1, Years));
+    auto bond2 = ext::make_shared<ZeroCouponBond>(3, TARGET(), 100.0, today + Period(2, Years));
+
+    std::vector<ext::shared_ptr<BondHelper>> helpers(2);
+    helpers[0] = ext::make_shared<BondHelper>(makeQuoteHandle(101.0), bond1);
+    helpers[1] = ext::make_shared<BondHelper>(makeQuoteHandle(102.0), bond2);
+
+    Real accuracy = 1e-10;      // default value
+    Size maxIterations = 10000; // default value
+    Array guess = {0.01};       // something positive so that initial value is in feasible region
+    
+    FlatZero unconstrainedMethod;
+    FittedBondDiscountCurve unconstrainedCurve(0, TARGET(), helpers, Actual365Fixed(), unconstrainedMethod, 
+                                               accuracy, maxIterations, guess);
+    BOOST_CHECK_LT(unconstrainedCurve.fitResults().solution()[0], 0.0);
+
+    FlatZero positiveMethod{PositiveConstraint()};
+    FittedBondDiscountCurve positiveCurve(0, TARGET(), helpers, Actual365Fixed(), positiveMethod,
+                                          accuracy, maxIterations, guess);
+    BOOST_CHECK_GT(positiveCurve.fitResults().solution()[0], 0.0);
+}
+
 BOOST_AUTO_TEST_SUITE_END()
 
 BOOST_AUTO_TEST_SUITE_END()