Python: multivariate nichtlineare Solver mit Einschränkungen

Question

Gegeben eine Funktion f(x), die einen Eingabevektor x nimmt und einen Vektor der gleichen Länge zurückgibt, wie können Sie die Nullstellen der Funktion Einstellung Einschränkungen auf x finden? (Z. B. ein Bereich für jede Komponente von x.)

Zu meiner Überraschung konnte ich nicht viele nützliche Informationen darüber finden. In der Scipy-Liste für Optimization and Root finding algorithms scheint es einige Optionen für Skalarfunktionen wie brentq zu geben. Ich kann jedoch keine Algorithmen finden, die eine solche Option für den multivariaten Fall unterstützen.

Natürlich könnte man einen Workaround wie quadrieren jede Komponente des zurückgegebenen Vektors und dann einen der Minimizer wie differential_evolution (das ist das einzige, was ich eigentlich denke). Ich kann mir nicht vorstellen, dass dies eine gute Strategie ist, da sie die quadratische Konvergenz von Newtons Algorithmus tötet. Auch finde ich es wirklich überraschend, dass es dafür keine Option zu geben scheint, da es ein wirklich häufiges Problem sein muss. Habe ich etwas verpasst?

Answer 1

Ein (nicht besonders schön, aber hoffentlich arbeiten) Option, um dieses Problem eine Funktion dem Solver geben würde arbeiten, die nur Wurzeln im eingeschränkten Bereich hat und das in einer Art und Weise fortgesetzt wird, um sicherzustellen, dass der Solver ist zurück in die richtige Region geschoben (ein bisschen wie here aber in mehreren Dimensionen).

Was könnte man dies (für rechteckige Einschränkungen zumindest) zu erreichen, tun ein constrainedFunction zu implementieren ist, die linear von der Grenze Wert Ihrer Funktion beginnend fortgesetzt wird:

import numpy as np 

def constrainedFunction(x, f, lower, upper, minIncr=0.001): 
    x = np.asarray(x) 
    lower = np.asarray(lower) 
    upper = np.asarray(upper) 
    xBorder = np.where(x<lower, lower, x) 
    xBorder = np.where(x>upper, upper, xBorder) 
    fBorder = f(xBorder) 
    distFromBorder = (np.sum(np.where(x<lower, lower-x, 0.)) 
         +np.sum(np.where(x>upper, x-upper, 0.))) 
    return (fBorder + (fBorder 
         +np.where(fBorder>0, minIncr, -minIncr))*distFromBorder) 

können Sie diese Funktion übergeben ein x Wert, die Funktion f, die Sie fortsetzen möchten, sowie zwei Arrays lower und upper der gleichen Form wie x gibt die unteren und oberen Grenzen in allen Dimensionen. Jetzt können Sie diese Funktion und nicht Ihre ursprüngliche Funktion an den Löser übergeben, um die Wurzeln zu finden.

Die Steilheit der Fortsetzung wird im Moment einfach als Grenzwert genommen, um steile Sprünge für Vorzeichenwechsel an der Grenze zu verhindern. Um Wurzeln außerhalb des beschränkten Bereichs zu verhindern, wird ein kleiner Wert zu positiven/negativen Grenzwerten addiert/subtrahiert. Ich stimme zu, dass dies kein sehr schöner Weg ist, aber es scheint zu funktionieren.

Hier sind zwei Beispiele. Für beide liegt die anfängliche Schätzung außerhalb der eingeschränkten Region, aber eine korrekte Wurzel in der eingeschränkten Region wird gefunden.

Finden der Wurzeln eines mehrdimensionalen Cosinus beschränkt auf [-2, -1] x [1, 2] ergibt:

from scipy import optimize as opt 

opt.root(constrainedFunction, x0=np.zeros(2), 
     args=(np.cos, np.asarray([-2., 1.]), np.asarray([-1, 2.]))) 

gibt:

fjac: array([[ -9.99999975e-01, 2.22992740e-04], 
     [ 2.22992740e-04, 9.99999975e-01]]) 
    fun: array([ 6.12323400e-17, 6.12323400e-17]) 
message: 'The solution converged.' 
    nfev: 11 
    qtf: array([ -2.50050470e-10, -1.98160617e-11]) 
     r: array([-1.00281376, 0.03518108, -0.9971942 ]) 
    status: 1 
success: True 
     x: array([-1.57079633, 1.57079633]) 

Dies gilt auch für Funktionen Arbeiten, sind nicht diagonal:

def f(x): 
    return np.asarray([0., np.cos(x.sum())]) 

opt.root(constrainedFunction, x0=np.zeros(2), 
     args=(f, np.asarray([-2., 2.]), np.asarray([-1, 4.]))) 

gibt:

fjac: array([[ 0.00254922, 0.99999675], 
     [-0.99999675, 0.00254922]]) 
    fun: array([ 0.00000000e+00, 6.12323400e-17]) 
message: 'The solution converged.' 
    nfev: 11 
    qtf: array([ 1.63189544e-11, 4.16007911e-14]) 
     r: array([-0.75738638, -0.99212138, -0.00246647]) 
    status: 1 
success: True 
     x: array([-1.65863336, 3.22942968]) 

Answer 2

Wenn Sie eine Optimierung mit Nebenbedingungen zu handhaben möchten, können Sie die leichte lirbary verwenden, die viel einfacher als scipy.optimize ist

Hier ist der Link zum Paket:

https://pypi.python.org/pypi/facile/1.2

So verwenden Sie die einfache Bibliothek für Ihr Beispiel. Sie müssen verfeinern, was ich hier schreibe, was nur allgemein ist. Wenn Sie Fehler haben, sagen Sie mir was.

import facile 

# Your vector x 

x = [ facile.variable('name', min, max) for i in range(Size) ] 


# I give an example here of your vector being ordered and each component in a range 
# You could as well put in the range where declaring variables 

for i in range(len(x)-1): 
    facile.constraint(x[i] < x[i+1]) 
    facile.constraint(range[i,0] < x[i] < range[i,1]) #Supposed you have a 'range' array where you store the range for each variable 


def function(...) 
# Define here the function you want to find roots of 


# Add as constraint that you want the vector to be a root of function 
facile.constraint(function(x) == 0) 


# Use facile solver 
if facile.solve(x): 
    print [x[i].value() for i in range(len(x))] 
else: 
    print "Impossible to find roots" 

Answer 3

Auf die Gefahr, etwas vorzuschlagen, das Sie vielleicht bereits durchgestrichen haben, glaube ich, dass dies nur mit möglich sein sollte. Der Catch ist, dass die Funktion nur ein Argument haben muss, aber dieses Argument kann ein Vektor/Liste sein.

Also f (x, y) wird nur f (z) wo z = [x, y].

Ein gutes Beispiel, das Sie nützlich finden könnten, wenn Sie nicht gefunden haben, ist here.

Wenn Sie Grenzen auferlegen wollen, wie Sie bereits erwähnt, für ein 2x1-Vektor, könnten Sie verwenden:

# Specify a (lower, upper) tuple for each component of the vector  
bnds = [(0., 1.) for i in len(x)] 

Und das innerhalb minimize als bounds Parameter verwenden.