2015-08-27 10 views
7

Ich versuche, die folgende Funktion mit Python scipy.optimize zu maximieren. Nach vielen Versuchen scheint es jedoch nicht zu funktionieren. Die Funktion und mein Code sind unten eingefügt. Danke fürs Helfen!Optimierung mit Python (scipy.optimize)

Problem

Maximize [sum (x_i/y_i)**gamma]**(1/gamma) 
subject to the constraint sum x_i = 1; x_i is in the interval (0,1). 

x ist ein Vektor der Wahl Variablen; y ist ein Vektor von Parametern; gamma ist ein Parameter. Die x s müssen zu eins addiert werden. Und jeder x muss im Intervall (0,1) sein.

-Code

def objective_function(x, y): 
    sum_contributions = 0 
    gamma = 0.2 

    for count in xrange(len(x)): 
     sum_contributions += (x[count]/y[count]) ** gamma 
    value = math.pow(sum_contributions, 1/gamma) 
    return -value 

cons = ({'type': 'eq', 'fun': lambda x: np.array([sum(x) - 1])}) 

y = [0.5, 0.3, 0.2] 
initial_x = [0.2, 0.3, 0.5] 

opt = minimize(objective_function, initial_x, args=(y,), method='SLSQP', 
constraints=cons,bounds=[(0, 1)] * len(x)) 
+0

Es könnte hilfreich sein, wenn Sie über spezifischere waren, was viele versuchen, "brachte das Ziel zu beleuchten deiner Frage. – Erik

+0

Ihr Code funktioniert für mich. Was ist das Problem, das du bekommst? Ich bekomme optimale 'x_opt: array ([0.29465573, 0.33480638, 0.37053789])'. Alles, was ich ändern musste, war 'bounds' sollte' len (initial_x) 'oder' len (y) 'enthalten, da' x' in deinem Code nicht definiert ist. – askewchan

+1

@askewchan, Könnte ein numerisches Stabilitätsproblem sein. Auf meinem Mac bekomme ich 'nan''nan''nan',' 'Iterationslimit überschritten'' –

Antwort

3

Manchmal macht numerische Optimierer aus irgendeinem Grund nicht funktionieren. Wir können das Problem etwas anders parametrisieren und es wird einfach funktionieren. (Und vielleicht schneller arbeiten)

Zum Beispiel für Grenzen (0,1), können wir eine Transformationsfunktion, so dass Werte in (-inf, +inf) nach, transformiert wird, werden bis in (0,1)

Ende Wir können einen ähnlichen Trick mit die Gleichheitsbeschränkungen. Zum Beispiel können wir die Dimension von 3 auf 2 reduzieren, da das letzte Element in x1-sum(x) sein muss.

Wenn es immer noch nicht funktioniert, können wir zu einem Optimierer wechseln, der keine Informationen von Derivaten benötigt, wie .

Und auch dort ist Lagrange multiplier.

In [111]: 

def trans_x(x): 
    x1 = x**2/(1+x**2) 
    z = np.hstack((x1, 1-sum(x1))) 
    return z 

def F(x, y, gamma = 0.2): 
    z = trans_x(x) 
    return -(((z/y)**gamma).sum())**(1./gamma) 
In [112]: 

opt = minimize(F, np.array([0., 1.]), args=(np.array(y),), 
       method='Nelder-Mead') 
opt 
Out[112]: 
    status: 0 
    nfev: 96 
success: True 
    fun: -265.27701747828007 
     x: array([ 0.6463264, 0.7094782]) 
message: 'Optimization terminated successfully.' 
    nit: 52 

Das Ergebnis ist:

In [113]: 

trans_x(opt.x) 
Out[113]: 
array([ 0.29465097, 0.33482303, 0.37052601]) 

Und wir können es sichtbar zu machen, mit:

In [114]: 

x1 = np.linspace(0,1) 
y1 = np.linspace(0,1) 
X,Y = np.meshgrid(x1,y1) 
Z = np.array([F(item, y) for item 
       in np.vstack((X.ravel(), Y.ravel())).T]).reshape((len(x1), -1), order='F') 
Z = np.fliplr(Z) 
Z = np.flipud(Z) 
plt.contourf(X, Y, Z, 50) 
plt.colorbar() 

enter image description here

+0

Hi @CT Zhu, geniale Lösung. Was macht das "x1 = x ** 2/(1 + x ** 2)"? –

+1

für 'x' in' (-inf, + inf) ', 'transformiere es in' x1' in '(0, 1)'. Ähnlich wie das Setzen von Grenzen, da das OP verlangt, dass x im Intervall (0,1) sein muss. –

+0

Es wandelt tatsächlich -inf in 1 um! Es funktioniert nur für (0, + inf) –

1

Selbst schwere diese Fragen ist ein bisschen veraltet Ich wollte hinzufügen eine alternative Lösung, die für andere hilfreich sein könnte, wenn sie in Zukunft auf diese Frage stößt .

Es macht unser Problem ist analytisch lösbar. Sie können durch Aufschreiben der Lagrange-Funktion des (gleichungsrestringierte) Optimierungsproblems starten:

L = \sum_i (x_i/y_i)^\gamma - \lambda (\sum x_i - 1) 

Die optimale Lösung wird durch Setzen der erste Ableitung dieser Lagrange-Funktion auf Null gefunden:

0 = \partial L/\partial x_i = \gamma x_i^{\gamma-1}/\y_i - \lambda 
=> x_i \propto y_i^{\gamma/(\gamma - 1)} 

Mit dieser Einsicht

In [4]: 
def analytical(y, gamma=0.2): 
    x = y**(gamma/(gamma-1.0)) 
    x /= np.sum(x) 
    return x 
xanalytical = analytical(y) 
xanalytical, objective_function(xanalytical, y) 
Out [4]: 
(array([ 0.29466774, 0.33480719, 0.37052507]), -265.27701765929692) 

CT Zhus Lösung ist elegant, aber es könnte die Positivität Einschränkung für die dritte coo verletzt: kann das Optimierungsproblem einfach und effizient durch gelöst werden raten.Für gamma = 0.2 scheint dies kein Problem in der Praxis zu sein, aber für verschiedene Gammas laufen Sie leicht in Schwierigkeiten:

In [5]: 
y = [0.2, 0.1, 0.8] 
opt = minimize(F, np.array([0., 1.]), args=(np.array(y), 2.0), 
       method='Nelder-Mead') 
trans_x(opt.x), opt.fun 
Out [5]: 
(array([ 1., 1., -1.]), -11.249999999999998) 

für andere Optimierungsprobleme mit den gleichen Wahrscheinlichkeit simplex Einschränkungen wie Ihr Problem, aber für die es keine analytische Lösung, könnte es sich lohnen, in projizierte Gradientenmethoden oder ähnliches zu schauen. Diese Methoden nutzen die Tatsache, dass es einen schnellen Algorithmus für die Projektion eines beliebigen Punktes auf diesen Satz gibt, siehe https://en.wikipedia.org/wiki/Simplex#Projection_onto_the_standard_simplex.

(Um den vollständigen Code und eine bessere Darstellung der Gleichungen nehmen einen Blick auf die Jupyter Notebook http://nbviewer.jupyter.org/github/andim/pysnippets/blob/master/optimization-simplex-constraints.ipynb zu sehen)

+0

Andi, danke. Aber das Optimierungsproblem, das ich lösen wollte, unterscheidet sich etwas von dem, das Sie mit analytischen Mitteln lösen. Hier ist das vollständige Problem: http://math.stackexchange.com/questions/2168338/explicit-solution-to-optimization-problem Ich würde mich freuen, über analytische Lösungen und schnelle numerische Näherungswerte zu hören. Danke: - – user58925

+0

Wenn der dritte Parameter positiv bleiben muss, würde so etwas funktionieren? '# a = dritte Koordinate (Positivitätsbedingung)' 'a = abs (a)' – mikey

Verwandte Themen