How to Solve Optimisation Problems using Penalty Functions in Python
I am working on a implementing a simple quadratic optimisation problem:
$$\min _x \; {\underline{x}}^T Q {\underline{x}}$$ $$s.t. \,\quad {\underline{\mu}}^T{\underline{x}} = R^*$$ $$ \quad \quad \underline{1}^T\underline{x} = 1 $$ I also expect to go on to include the inequality constraint, as an additional complexity, once the above is working.
$$ x_i \geq 0$$ The method I think is simplest, and which I understand best for implementing these constraints, is the penalty function method, where we modify the objective function to 'steer' the optimisation away from forbidden regions. 490
By carefully parameterising the size of the penalties, I have achieved good results using SciPy's built-in Nelder-Mead Simplex algorithm, using the objective function below.
def objective(x):
Q = DF.cov() # Covariance matrix
# Penalty Function method
penalty1 = 0.0005 * abs(np.sum(x)-1) # Large for sum(x) <> 1
penalty2 = 0.05 * abs(R_min - np.matmul(Mus.transpose(), x)) # Large for returns <> R_min
return np.matmul(x.transpose(),np.matmul(Q,x)) + penalty1 + penalty2
Now, I wish to use other optimisation algorithms (in particular BFGS and Newton-CG), which require the gradient and Hessian of the objective function. I have implemented the derivative functions in the unconstrained case, but by adding the penalty terms to the objective (and the derivatives of the penalties to the gradient function) optimisation fails with the following error:
Warning: Desired error not necessarily achieved due to precision loss.
Current function value: 0.000056
Iterations: 0
Function evaluations: 780
Gradient evaluations: 96
(Previously Iterations would be a few hundred). This strictly occurs with penalty1, but not penalty2 by itself, so either my derivative is wrong for penalty1:
penalty1_der = np.sign(x)
Or can I not use the L1 norm in this way? I also tried replacing the constraint with a smoother, quadratic approximation:
penalty1 = np.matmul((x - vector_ones).transpose(), (x - vector_ones))
but unfortunately, although this prevents the error, Minimize() seems to completely ignore my penalty functions (even with vastly increased parameters).
How can I implement my constraints such that I can solve the problem using BFGS/Newton-CG? 509