I have an optimization problem where I'm trying to find an array that needs to optimize two functions simultaneously.
In the minimal example below I have two known arrays w
and x
and an unknown array y
. I initialize array y
to contains only 1s.
I then specify function np.sqrt(np.sum((x-np.array)**2)
and want to find the array y
where
np.sqrt(np.sum((x-y)**2)
approaches 5
np.sqrt(np.sum((w-y)**2)
approaches 8
The code below can be used to successfully optimize y
with respect to a single array, but I would like to find that the solution that optimizes y
with respect to both x
and y
simultaneously, but am unsure how to specify the two constraints.
y
should only consist of values greater than 0.
Any ideas on how to go about this ?
w = np.array([6, 3, 1, 0, 2])
x = np.array([3, 4, 5, 6, 7])
y = np.array([1, 1, 1, 1, 1])
def func(x, y):
z = np.sqrt(np.sum((x-y)**2)) - 5
return np.zeros(x.shape[0],) + z
r = opt.root(func, x0=y, method='hybr')
print(r.x)
# array([1.97522498 3.47287981 5.1943792 2.10120135 4.09593969])
print(np.sqrt(np.sum((x-r.x)**2)))
# 5.0
One option is to use scipy.optimize.minimize
instead of root
, Here you have multiple solver options and some of them (ie SLSQP
) allow you to specify multiple constraints. Note that I changed the variable names so that x
is the array you want to optimise and y
and z
define the constraints.
from scipy.optimize import minimize
import numpy as np
x0 = np.array([1, 1, 1, 1, 1])
y = np.array([6, 3, 1, 0, 2])
z = np.array([3, 4, 5, 6, 7])
constraint_x = dict(type='ineq',
fun=lambda x: x) # fulfilled if > 0
constraint_y = dict(type='eq',
fun=lambda x: np.linalg.norm(x-y) - 5) # fulfilled if == 0
constraint_z = dict(type='eq',
fun=lambda x: np.linalg.norm(x-z) - 8) # fulfilled if == 0
res = minimize(fun=lambda x: np.linalg.norm(x), constraints=[constraint_y, constraint_z], x0=x0,
method='SLSQP', options=dict(ftol=1e-8)) # default 1e-6
print(res.x) # [1.55517124 1.44981672 1.46921122 1.61335466 2.13174483]
print(np.linalg.norm(res.x-y)) # 5.00000000137866
print(np.linalg.norm(res.x-z)) # 8.000000000930026
This is a minimizer so besides the constraints it also wants a function to minimize, I chose just the norm of y
, but setting the function to a constant (ie lambda x: 1) would have also worked.
Note also that the constraints are not exactly fulfilled, you can increase the accuracy by setting optional argument ftol
to a smaller value ie 1e-10
.
For more information see also the documentation and the corresponding sections for each solver.