Question

Many real-world optimization problems have constraints - for example, a set of parameters may have to sum to 1.0 (eqquality constraint), or some parameters may have to be non-negative (inequality constraint). Sometimes, the constraints can be incorporated into the function to be minimized, for example, the non-negativity constraint \(p > 0\) can be removed by substituting \(p = e^q\) and optimizing for \(q\). Using such workarounds, it may be possible to convert a constrained optimization problem into an unconstrained one, and use the methods discussed above to sovle the problem.

Alternatively, we can use optimization methods that allow the speicification of constraints directly in the problem statement as shown in this section. Internally, constraint violation penalties, barriers and Lagrange multpiliers are some of the methods used used to handle these constraints. We use the example provided in the Scipy tutorial to illustrate how to set constraints.

\[f(x) = -(2xy + 2x - x^2 -2y^2)\]

subject to the constraint

\[\begin{split} x^3 - y = 0 \\ y - (x-1)^4 - 2 \ge 0\end{split}\]\[and the bounds\] \[\begin{split}0.5 \le x \le 1.5 \\ 1.5 \le y \le 2.5\end{split}\]

def f(x):
    return -(2*x[0]*x[1] + 2*x[0] - x[0]**2 - 2*x[1]**2)

x = np.linspace(0, 3, 100)
y = np.linspace(0, 3, 100)
X, Y = np.meshgrid(x, y)
Z = f(np.vstack([X.ravel(), Y.ravel()])).reshape((100,100))
plt.contour(X, Y, Z, np.arange(-1.99,10, 1));
plt.plot(x, x**3, 'k:', linewidth=1)
plt.plot(x, (x-1)**4+2, 'k:', linewidth=1)
plt.fill([0.5,0.5,1.5,1.5], [2.5,1.5,1.5,2.5], alpha=0.3)
plt.axis([0,3,0,3])

[0, 3, 0, 3]

To set consttarints, we pass in a dictionary with keys ty;pe , fun and jac . Note that the inequlaity cosntraint assumes a \(C_j x \ge 0\) form. As usual, the jac is optional and will be numerically estimted if not provided.

cons = ({'type': 'eq',
         'fun' : lambda x: np.array([x[0]**3 - x[1]]),
         'jac' : lambda x: np.array([3.0*(x[0]**2.0), -1.0])},
        {'type': 'ineq',
         'fun' : lambda x: np.array([x[1] - (x[0]-1)**4 - 2])})

bnds = ((0.5, 1.5), (1.5, 2.5))

x0 = [0, 2.5]

Unconstrained optimization

ux = opt.minimize(f, x0, constraints=None)
ux

  status: 0
 success: True
    njev: 5
    nfev: 20
hess_inv: array([[ 1. ,  0.5],
      [ 0.5,  0.5]])
     fun: -1.9999999999999987
       x: array([ 2.,  1.])
 message: 'Optimization terminated successfully.'
     jac: array([ 0.,  0.])

Constrained optimization

cx = opt.minimize(f, x0, bounds=bnds, constraints=cons)
cx

 status: 0
success: True
   njev: 5
   nfev: 21
    fun: 2.0499154720925521
      x: array([ 1.2609,  2.0046])
message: 'Optimization terminated successfully.'
    jac: array([-3.4875,  5.4967,  0.    ])
    nit: 5

x = np.linspace(0, 3, 100)
y = np.linspace(0, 3, 100)
X, Y = np.meshgrid(x, y)
Z = f(np.vstack([X.ravel(), Y.ravel()])).reshape((100,100))
plt.contour(X, Y, Z, np.arange(-1.99,10, 1));
plt.plot(x, x**3, 'k:', linewidth=1)
plt.plot(x, (x-1)**4+2, 'k:', linewidth=1)
plt.text(ux['x'][0], ux['x'][1], 'x', va='center', ha='center', size=20, color='blue')
plt.text(cx['x'][0], cx['x'][1], 'x', va='center', ha='center', size=20, color='red')
plt.fill([0.5,0.5,1.5,1.5], [2.5,1.5,1.5,2.5], alpha=0.3)
plt.axis([0,3,0,3]);