1. 程式人生 > >利用Python求解帶約束的最優化問題

利用Python求解帶約束的最優化問題

題目:

\mathbf{min_x 10-x_1^2 - x_2^2, subject \; to\; x_2\geq x_1^2, x_1+x_2=0}

1. 利用拉格朗日乘子法

#匯入sympy包,用於求導,方程組求解等等
from sympy import * 

#設定變數
x1 = symbols("x1")
x2 = symbols("x2")
alpha = symbols("alpha")
beta = symbols("beta")

#構造拉格朗日等式
L = 10 - x1*x1 - x2*x2 + alpha * (x1*x1 - x2) + beta * (x1 + x2)

#求導,構造KKT條件
difyL_x1 = diff(L, x1)  #對變數x1求導
difyL_x2 = diff(L, x2)  #對變數x2求導
difyL_beta = diff(L, beta)  #對乘子beta求導
dualCpt = alpha * (x1 * x1 - x2)  #對偶互補條件

#求解KKT等式
aa = solve([difyL_x1, difyL_x2, difyL_beta, dualCpt], [x1, x2, alpha, beta])

#列印結果,還需驗證alpha>=0和不等式約束<=0
for i in aa:
    if i[2] >= 0:
        if (i[0]**2 - i[1]) <= 0:
            print(i)

結果:

(-1, 1, 4, 6)
(0, 0, 0, 0)

2. scipy包裡面的minimize函式求解

from scipy.optimize import minimize
import numpy as np 

from mpl_toolkits.mplot3d import Axes3D
from matplotlib import pyplot as plt 

#目標函式:
def func(args):
    fun = lambda x: 10 - x[0]**2 - x[1]**2
    return fun

#約束條件,包括等式約束和不等式約束
def con(args):
    cons = ({'type': 'ineq', 'fun': lambda x: x[1]-x[0]**2},
            {'type': 'eq', 'fun': lambda x: x[0]+x[1]})
    return cons 

#畫三維模式圖
def draw3D():
    fig = plt.figure()
    ax = Axes3D(fig)
    x_arange = np.arange(-5.0, 5.0)
    y_arange = np.arange(-5.0, 5.0)
    X, Y = np.meshgrid(x_arange, y_arange)
    Z1 = 10 - X**2 - Y**2
    Z2 = Y - X**2
    Z3 = X + Y
    plt.xlabel('x')
    plt.ylabel('y')
    ax.plot_surface(X, Y, Z1, rstride=1, cstride=1, cmap='rainbow')
    ax.plot_surface(X, Y, Z2, rstride=1, cstride=1, cmap='rainbow')
    ax.plot_surface(X, Y, Z3, rstride=1, cstride=1, cmap='rainbow')
    plt.show()

#畫等高線圖
def drawContour():
    x_arange = np.linspace(-3.0, 4.0, 256)
    y_arange = np.linspace(-3.0, 4.0, 256)
    X, Y = np.meshgrid(x_arange, y_arange)
    Z1 = 10 - X**2 - Y**2
    Z2 = Y - X**2
    Z3 = X + Y
    plt.xlabel('x')
    plt.ylabel('y')
    plt.contourf(X, Y, Z1, 8, alpha=0.75, cmap='rainbow')
    plt.contourf(X, Y, Z2, 8, alpha=0.75, cmap='rainbow')
    plt.contourf(X, Y, Z3, 8, alpha=0.75, cmap='rainbow')
    C1 = plt.contour(X, Y, Z1, 8, colors='black')
    C2 = plt.contour(X, Y, Z2, 8, colors='blue')
    C3 = plt.contour(X, Y, Z3, 8, colors='red')
    plt.clabel(C1, inline=1, fontsize=10)
    plt.clabel(C2, inline=1, fontsize=10)
    plt.clabel(C3, inline=1, fontsize=10)
    plt.show()


if __name__ == "__main__":
    args = ()
    args1 = ()
    cons = con(args1)
    x0 = np.array((1.0, 2.0))  #設定初始值,初始值的設定很重要,很容易收斂到另外的極值點中,建議多試幾個值
    
    #求解#
    res = minimize(func(args), x0, method='SLSQP', constraints=cons)
    #####
    print(res.fun)
    print(res.success)
    print(res.x)

    # draw3D()
    drawContour()

結果:

7.99999990708696
True
[-1.00000002  1.00000002]