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)