Penalty function for nonlinear inequality constraints in Mystic are evaluated outside of bounds

I would like to use mystic solver to solve the following nonlinear optimisation problem with nonlinear constraints. Here the code:

import numpy as np
import matplotlib.pyplot as plt
from math import sqrt
from mystic.solvers import diffev2, fmin, fmin_powell
from mystic.monitors import VerboseMonitor
from mystic.penalty import quadratic_inequality, quadratic_equality

def pos_scale(c, q):
    return 1.0 / (1 + c*sqrt(q))

def omega_scaled(w, c, q):
    return min(w, pos_scale(c, q))

def constraints(q1, q2, w1, w2, c1, c2, fx1, fx2):
    #print('{} {}'.format(q1, q2))
    v1 = omega_scaled(w1, c1, q1)*q1*fx1
    v2 = omega_scaled(w2, c2, q2)*q2*fx2
    return v1 + v2

constraints_f = lambda q1, q2: constraints(q1, q2, 0.95, 0.92, 0.06, 0.05, 10000, 1000)
constraints_v = np.vectorize(constraints_f)

def cost(q1, q2, w1, w2, c1, c2, fx1, fx2):
    v1 = (1-omega_scaled(w1, c1, q1))*q1*fx1
    v2 = (1-omega_scaled(w2, c2, q2))*q2*fx2
    return v1 + v2

cost_f = lambda q1, q2: cost(q1, q2, 0.95, 0.92, 0.06, 0.04, 10000, 1000)
cost_v = np.vectorize(cost_f)

objective = lambda x: cost_v(x[0], x[1]).item()

def penalty_value(x, target):
    return target - constraints_f(x[0], x[1])

@quadratic_inequality(penalty_value, kwds={'target': 200000.0})
def penalty(x):
  return 0.0

I model the constraints with a quadratic penalty. The constraints require the domain to be the positive orthant.

mon = VerboseMonitor(10)
bounds=[(0, 50), (0, 300)]
result = fmin(objective, x0=[15, 150], bounds=bounds, penalty=penalty, 
              npop=10, gtol=200, disp=False, full_output=True, itermon=mon, maxiter=500)

result


Generation 0 has ChiSquare: 77606.160271
Generation 10 has ChiSquare: 62080.449073
Generation 20 has ChiSquare: 55726.285526
Generation 30 has ChiSquare: 55505.829370
Generation 40 has ChiSquare: 55478.612377
Generation 50 has ChiSquare: 55475.462051
Generation 60 has ChiSquare: 55474.597220
Generation 70 has ChiSquare: 55474.532390
Generation 80 has ChiSquare: 55474.530891
Generation 90 has ChiSquare: 55474.530773
STOP("CandidateRelativeTolerance with {'xtol': 0.0001, 'ftol': 0.0001}")
(array([21.50326424, 42.0783277 ]), 55474.53077292251, 98, 177, 0)

The solver finds the optimal solution if I use a reasonable starting value. However, another starting value (or another solver like the Powell solver) triggers in the step search a call to the constraints function outside of the valid domain.

How can I best force the constraints function in the penalty to be restricted to the bounds that I provide to the solver? Should that not be checked by the solver itself? Or do I need to take care of that myself also in the constraints function?

To visualize the solution:

fig = plt.figure(figsize=(12,6))
left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
ax = fig.add_axes([left, bottom, width, height])

q1 = np.linspace(0.1, 50, 100)
q2 = np.linspace(1, 300, 100)
X, Y = np.meshgrid(q1, q2)

Z = constraints_v(X, Y)
cp1 = plt.contour(X, Y, Z, 20, colors='black', linestyles='dashed')
cp2 = plt.contour(X, Y, Z, [200000], colors='white', linestyles='solid')
plt.clabel(cp2, inline=True, fontsize=12)

Z = cost_v(X, Y)
cp3 = plt.contourf(X, Y, Z, 25)
plt.colorbar(cp3)

sol = list(result[0])
plt.plot(sol[0], sol[1], 'go--', linewidth=2, markersize=14)

I'm the mystic author. Penalty functions are essentially soft constraints, so they can be violated. When they are, they will add a penalty to the cost. If you want to restrict the input values explicitly, then you want a hard constraint... given with the constraints keyword. So, add the following to ensure that the candidate solutions are always chosen from the positive orthant (i.e. from within your chosen bounds).

...[snip]...

import mystic.symbolic as ms
from mystic.constraints import and_
bounds = '''
x0 >= 0
x0 <= 50
x1 >= 0
x1 <= 300
'''
eqn = ms.simplify(bounds, all=True)
cons = ms.generate_constraint(ms.generate_solvers(eqn), join=and_)

mon = VerboseMonitor(10)
bounds=[(0, 50), (0, 300)]
result = fmin_powell(objective, x0=[15, 150], bounds=bounds, penalty=penalty,
                     npop=10, gtol=200, disp=False, full_output=True,
                     constraints=cons, itermon=mon, maxiter=500)

I have a feature request open to do this automatically, or at least with a toggle/flag -- but currently, you need to manually add the bounds to the hard constraints.