Solución:
Mientras que la SLSQP algoritmo en scipy.optimize.minimize es bueno, tiene un montón de limitaciones. El primero de los cuales es un QP solucionador, por lo que funciona para ecuaciones que encajan bien en un paradigma de programación cuadrática. Pero, ¿qué sucede si tiene limitaciones funcionales? También, scipy.optimize.minimize no es un optimizador global, por lo que a menudo es necesario comenzar muy cerca de los resultados finales.
Hay un paquete de optimización no lineal restringido (llamado mystic) que ha existido durante casi tanto tiempo como scipy.optimize en sí mismo: lo sugeriría como el método de referencia para manejar cualquier optimización no lineal restringida general.
Por ejemplo, su problema, si entiendo su pseudocódigo, se parece a esto:
import numpy as np
M = 10
N = 3
Q = 10
C = 10
# let's be lazy, and generate s and u randomly...
s = np.random.randint(-Q,Q, size=(M,N,N))
u = np.random.randint(-Q,Q, size=(M,N))
def percentile(p, x):
x = np.sort(x)
p = 0.01 * p * len(x)
if int(p) != p:
return x[int(np.floor(p))]
p = int(p)
return x[p:p+2].mean()
def objective(x, p=5): # inverted objective, to find the max
return -1*percentile(p, [np.dot(np.atleast_2d(u[i]), x)[0] for i in range(0,M-1)])
def constraint(x, p=95, v=C): # 95%(xTsx) - v <= 0
x = np.atleast_2d(x)
return percentile(p, [np.dot(np.dot(x,s[i]),x.T)[0,0] for i in range(0,M-1)]) - v
bounds = [(0,1) for i in range(0,N)]
Entonces, para manejar su problema en mystic, solo necesita especificar los límites y las restricciones.
from mystic.penalty import quadratic_inequality
@quadratic_inequality(constraint, k=1e4)
def penalty(x):
return 0.0
from mystic.solvers import diffev2
from mystic.monitors import VerboseMonitor
mon = VerboseMonitor(10)
result = diffev2(objective, x0=bounds, penalty=penalty, npop=10, gtol=200,
disp=False, full_output=True, itermon=mon, maxiter=M*N*100)
print result[0]
print result[1]
El resultado se parece a esto:
Generation 0 has Chi-Squared: -0.434718
Generation 10 has Chi-Squared: -1.733787
Generation 20 has Chi-Squared: -1.859787
Generation 30 has Chi-Squared: -1.860533
Generation 40 has Chi-Squared: -1.860533
Generation 50 has Chi-Squared: -1.860533
Generation 60 has Chi-Squared: -1.860533
Generation 70 has Chi-Squared: -1.860533
Generation 80 has Chi-Squared: -1.860533
Generation 90 has Chi-Squared: -1.860533
Generation 100 has Chi-Squared: -1.860533
Generation 110 has Chi-Squared: -1.860533
Generation 120 has Chi-Squared: -1.860533
Generation 130 has Chi-Squared: -1.860533
Generation 140 has Chi-Squared: -1.860533
Generation 150 has Chi-Squared: -1.860533
Generation 160 has Chi-Squared: -1.860533
Generation 170 has Chi-Squared: -1.860533
Generation 180 has Chi-Squared: -1.860533
Generation 190 has Chi-Squared: -1.860533
Generation 200 has Chi-Squared: -1.860533
Generation 210 has Chi-Squared: -1.860533
STOP("ChangeOverGeneration with {'tolerance': 0.005, 'generations': 200}")
[-0.17207128 0.73183465 -0.28218955]
-1.86053344078
mystic es muy flexible y puede manejar cualquier tipo de restricciones (por ejemplo, igualdades, desigualdades), incluidas las restricciones simbólicas y funcionales. Especifiqué las restricciones como “penalizaciones” arriba, que es la forma tradicional, en el sentido de que aplican una penalización al objetivo cuando se viola la restricción.
mystic también proporciona transformaciones de kernel no lineales, que restringen el espacio de la solución al reducir el espacio de soluciones válidas (es decir, mediante un mapeo espacial o transformación de kernel).
Como ejemplo, aquí está mystic resolver un problema que rompe muchos solucionadores de QP, ya que las restricciones no tienen la forma de una matriz de restricciones. Está optimizando el diseño de un recipiente a presión.
"Pressure Vessel Design"
def objective(x):
x0,x1,x2,x3 = x
return 0.6224*x0*x2*x3 + 1.7781*x1*x2**2 + 3.1661*x0**2*x3 + 19.84*x0**2*x2
bounds = [(0,1e6)]*4
# with penalty='penalty' applied, solution is:
xs = [0.72759093, 0.35964857, 37.69901188, 240.0]
ys = 5804.3762083
from mystic.symbolic import generate_constraint, generate_solvers, simplify
from mystic.symbolic import generate_penalty, generate_conditions
equations = """
-x0 + 0.0193*x2 <= 0.0
-x1 + 0.00954*x2 <= 0.0
-pi*x2**2*x3 - (4/3.)*pi*x2**3 + 1296000.0 <= 0.0
x3 - 240.0 <= 0.0
"""
cf = generate_constraint(generate_solvers(simplify(equations)))
pf = generate_penalty(generate_conditions(equations), k=1e12)
if __name__ == '__main__':
from mystic.solvers import diffev2
from mystic.math import almostEqual
from mystic.monitors import VerboseMonitor
mon = VerboseMonitor(10)
result = diffev2(objective, x0=bounds, bounds=bounds, constraints=cf, penalty=pf,
npop=40, gtol=50, disp=False, full_output=True, itermon=mon)
assert almostEqual(result[0], xs, rel=1e-2)
assert almostEqual(result[1], ys, rel=1e-2)
Encuentre esto, y aproximadamente 100 ejemplos similares, aquí: https://github.com/uqfoundation/mystic.
Soy el autor, por lo que soy un poco parcial. Sin embargo, el sesgo es muy leve. mystic es a la vez maduro y bien respaldado, y no tiene paralelo en su capacidad para resolver problemas de optimización no lineal con restricciones estrictas.
scipy tiene un paquete espectacular para la optimización no lineal restringida.
Puede comenzar leyendo el optimize doc, pero aquí hay un ejemplo con SLSQP:
minimize(func, [-1.0,1.0], args=(-1.0,), jac=func_deriv, constraints=cons, method='SLSQP', options={'disp': True})
Como también han comentado otros, el paquete de minimización de SciPy es un buen lugar para comenzar. También tenemos una revisión de muchos otros paquetes de optimización en el artículo de Python Gekko (consulte la Sección 4). He incluido un ejemplo a continuación (punto de referencia Hock Schittkowski # 71) que incluye una función objetivo, restricción de igualdad y restricción de desigualdad en Scipy.optimize.minimize.
import numpy as np
from scipy.optimize import minimize
def objective(x):
return x[0]*x[3]*(x[0]+x[1]+x[2])+x[2]
def constraint1(x):
return x[0]*x[1]*x[2]*x[3]-25.0
def constraint2(x):
sum_eq = 40.0
for i in range(4):
sum_eq = sum_eq - x[i]**2
return sum_eq
# initial guesses
n = 4
x0 = np.zeros(n)
x0[0] = 1.0
x0[1] = 5.0
x0[2] = 5.0
x0[3] = 1.0
# show initial objective
print('Initial SSE Objective: ' + str(objective(x0)))
# optimize
b = (1.0,5.0)
bnds = (b, b, b, b)
con1 = {'type': 'ineq', 'fun': constraint1}
con2 = {'type': 'eq', 'fun': constraint2}
cons = ([con1,con2])
solution = minimize(objective,x0,method='SLSQP',
bounds=bnds,constraints=cons)
x = solution.x
# show final objective
print('Final SSE Objective: ' + str(objective(x)))
# print solution
print('Solution')
print('x1 = ' + str(x[0]))
print('x2 = ' + str(x[1]))
print('x3 = ' + str(x[2]))
print('x4 = ' + str(x[3]))
Aquí está el mismo problema con Python Gekko:
from gekko import GEKKO
m = GEKKO()
x1,x2,x3,x4 = m.Array(m.Var,4,lb=1,ub=5)
x1.value = 1; x2.value = 5; x3.value = 5; x4.value = 1
m.Equation(x1*x2*x3*x4>=25)
m.Equation(x1**2+x2**2+x3**2+x4**2==40)
m.Minimize(x1*x4*(x1+x2+x3)+x3)
m.solve(disp=False)
print(x1.value,x2.value,x3.value,x4.value)
También hay un hilo de discusión más completo sobre solucionadores de programación no lineal para Python si SLSQP no puede resolver su problema. El material de mi curso sobre Optimización del diseño de ingeniería está disponible si necesita información adicional sobre los métodos de resolución.