Python解决常见最优化问题的方法和策略

python求解常见的最优化问题

多变量最优化
问题:
彩电商准备推出两种产品,19寸彩电,建议价格339美元(成本195美元),21寸彩电,建议价格399美元(成本225美元),
固定成本400000美元,已知每售出一台彩电,会导致彩电价格下降1美分,而且每售出一台19寸彩电,会导致21寸彩电下降0.4美分;
每售出一台21寸彩电,会导致19寸彩电下降0.3美分,请问每台彩电应该生产多少台才能达到最大收益。

from sympy import *
from sympy.plotting import plot3d

"""
多变量最优化
问题:
彩电商准备推出两种产品,19寸彩电,建议价格339美元(成本195美元),21寸彩电,建议价格399美元(成本225美元),
固定成本400000美元,已知每售出一台彩电,会导致彩电价格下降1美分,而且每售出一台19寸彩电,会导致21寸彩电下降0.4美分;
每售出一台21寸彩电,会导致19寸彩电下降0.3美分,请问每台彩电应该生产多少台才能达到最大收益。
"""
x1 = symbols('x1')
x2 = symbols('x2')
fx = (339 - 0.01 * x1 - 0.003 * x2) * x1 + (399 - 0.01 * x2 - 0.004 * x1) * x2 - (400000 + 195 * x1 + 225 * x2)

# 分别对x1/x2求偏导数
expr1 = diff(fx, x1)
expr2 = diff(fx, x2)
print(expr1)
print(expr2)
# 解方程组x1=4735,x2= 7043
x1x2 = solve([expr1, expr2], [x1, x2])
nx1, nx2 = round(x1x2[x1]), round(x1x2[x2])
# 将x1,x2代入原方程得到最大收益=553641
result = fx.evalf(subs={x1: nx1, x2: nx2})
print(round(result))

# 是一个曲面
plot3d(fx, (x1, 0, 10000), (x2, 0, 10000))

img

1 求解带有约束的最小化问题

img

from scipy.optimize import minimize
e = 1e-10
fun = lambda x:(x[0] - 0.667) / (x[0] + x[1] + x[2] - 2) # 目标函数
cons = ({'type':'eq','fun':lambda x:x[0] * x[1] * x[2] - 1}, #约束条件:xyz=1
        {'type':'ineq','fun':lambda x:x[0] - e},             #x,y,z都分别大于0
        {'type':'ineq','fun':lambda x:x[1] - e},
        {'type':'ineq','fun':lambda x:x[2] - e})

x0 = np.array((1.0,1.0,1.0)) #设置xyz的初始值
res = minimize(fun,x0,method='SLSQP',constraints=cons)
print("最小值:",res.fun)
print("最优解:",res.x)
print("迭代终止是否成功",res.success)
print("迭代终止原因",res.message)

img

2 求解不带约束的最小化问题

使用工具:scipy.optimize.minimize

在使用scipy.optimize.minimize求解最小化时,优化器可以选择4种,分别是:

  • Nelder-Mead
  • Newton-CG
  • Powell
  • BFGS
  • BFGS算法的拓展阅读:https://machinelearningmastery.com/bfgs-optimization-in-python/

    Nelder-Mean算法的拓展阅读:https://machinelearningmastery.com

    img

    from scipy.optimize import minimize
    from numpy.random import rand
    
    # 两个未知数
    def objective(x):
        return x[0]**2.0 + x[1]**2
    
    r_min,r_max = -5.0,5.0
    
    # 生成随机初始值
    pt = r_min + rand(2) * (r_max - r_min)
    
    result = minimize(objective,pt,method='L-BFGS-B')
    
    print("Status: %s" %result['message'])
    print("Total Evaluation: %d" %result['nfev'])
    
    solution = result['x']
    evaluation = objective(solution)
    print("Solution:f(%s) = %.5f" %(solution,evaluation))
    print(result)
    
    

    img

    3 求解线性优化问题

    使用工具:scipy.optimize.line_search

    对于凸函数寻求线性求解。

    from numpy import arange
    from scipy.optimize import line_search
    from matplotlib import pyplot
    
    # 目标函数
    def objective(x):
        return (-5.0 + x)**2.0
    
    # 目标函数的导函数
    def gradient(x):
        return 2.0 * (-5.0 + x)
    
    # 定义开始点
    point = -5.0
    # 定义移动方位
    direction = 100.0
    # 打印初始条件
    print('start=%.1f, direction=%.1f' % (point, direction))
    
    # 执行线性优化函数,传入目标函数objective,导函数gradident,初始点point,移动方位direction
    result = line_search(objective, gradient, point, direction)
    
    # 得到结果,求解出alpha,前进的最优幅度
    alpha = result[0]
    print('Alpha: %.3f' % alpha)
    print('Function evaluations: %d' % result[1])
    
    # define objective function minima
    end = point + alpha * direction
    # evaluate objective function minima
    print('f(end) = f(%.3f) = %.3f' % (end, objective(end)))
    
    
    # define range
    r_min, r_max = -10.0, 20.0
    # prepare inputs
    inputs = arange(r_min, r_max, 0.1)
    # compute targets
    targets = [objective(x) for x in inputs]
    # plot inputs vs objective
    pyplot.plot(inputs, targets, '--', label='objective')
    # plot start and end of the search
    pyplot.plot([point], [objective(point)], 's', color='g')
    pyplot.plot([end], [objective(end)], 's', color='r')
    pyplot.legend()
    pyplot.show()
    

    img

    4 求解线性规划问题

    使用工具:scipy.optimize.linprog

    img

    import numpy as np
    from scipy.optimize import linprog
    
    C = np.array([0.4,0.5])                                # 目标函数的系数 
    A_ub = np.array([[0.3,0.1],[-0.6,-0.4]])               # 不等式约束的系数
    b_ub = np.array([2.7,-6])                              # 不等式约束的值
    A_eq = np.array([[0.5,0.5]])                           # 等式约束的系数
    b_eq = np.array([6])                                   # 等式约束的值
    r = linprog(C,A_ub,b_ub,A_eq,b_eq,bounds=((0,None),(0,None)))
    print(r)
    

    img

    5 求解全局最优化

    求解全局优化有若干种方法,scipy中提供了三种,分别是:

  • basinhopping() (盆地跳跃算法)
  • differential_evolution() (微分进化算法)
  • dual_annealing() (双重模拟退火算法)
  • img

    5.1 双重模拟退火算法(dual_annealing)

    使用工具:scipy.optimize.dual_annealing

    Similated Annealing算法的拓展阅读:

    https://machinelearningmastery.com/dual-annealing-optimization-with-python/

    https://machinelearningmastery.com

    # 导入优化包
    from scipy.optimize import dual_annealing,basinhopping,differential_evolution
    # simulated annealing global optimization for a multimodal objective function
    
    
    # objective function
    def objective(x):
        return (x[0]**2 + x[1] - 11)**2 + (x[0] + x[1]**2 - 7)**2
    
    # 定义x和y的范围
    r_min, r_max = -5.0, 5.0
    # define the bounds on the search
    bounds = [[r_min, r_max], [r_min, r_max]]
    
    # 执行模拟退火算法
    # perform the simulated annealing search
    result = dual_annealing(objective, bounds)
    # summarize the result
    print('Status : %s' % result['message'])
    print('Total Evaluations: %d' % result['nfev'])
    # evaluate solution
    solution = result['x']
    evaluation = objective(solution)
    print('Solution: f(%s) = %.5f' % (solution, evaluation))
    print(result)
    

    img

    5.2 盆地跳跃算法(basinhopping)

    使用工具:scipy.optimize.basinhopping

    5.2.1 basinhopping求解单目标

    img

    from scipy.optimize import basinhopping
    # 定义目标函数
    func = lambda x:np.cos(14.5 * x - 0.3) + (x + 0.2) * x
    # 定义带求解变量初始值
    x0 = [1.0]
    minimizer_kwargs = {"method":"BFGS"}                 #优化器使用BFGS
    
    ret = basinhopping(func,x0,minimizer_kwargs=minimizer_kwargs,niter=200)
    print("Global minimum x=%.4f ,f(x0) = %.4f" %(ret.x,ret.fun))
    print(ret)
    

    img

    这里的优化器也可以使用:Nelder-Mead、Powell

    5.2.2 basinhopping求解多目标

    img

    # simulated annealing global optimization for a multimodal objective function
    from scipy.optimize import dual_annealing,basinhopping,differential_evolution
    
    def func2d(x):
        f = np.cos(14.5 * x[0] - 0.3) + (x[1] + 0.2) * x[1] + (x[0] + 0.2) * x[0]
        df = np.zeros(2)
        df[0] = -14.5 * np.sin(14.5 * x[0] - 0.3) + 2. * x[0] + 0.2
        df[1] = 2 * x[1] + 0.2
        return f,df
    
    minimizer_kwargs = {"method":"L-BFGS-B","jac":True}
    x0 = [1.0,1.0]
    ret = basinhopping(func2d,x0,minimizer_kwargs=minimizer_kwargs,niter=200)
    print("Global minimum:x = [%.4f,%.4f],f(x0) = %.4f" %(ret.x[0],ret.x[1],ret.fun))
    

    img

    5.3 微分进化算法(differential_evolution)

    使用工具:scipy.optimize.differential_evolution

    5.3.1 不带约束

    img

    from scipy.optimize import differential_evolution
    import numpy as np
    
    # 定义目标函数
    def ackley(x):
        arg1 = -0.2 * np.sqrt(0.5 * (x[0] ** 2 + x[1] ** 2))
        arg2 = 0.5 * (np.cos(2. * np.pi * x[0]) + np.cos(2. * np.pi * x[1]))
        return -20. * np.exp(arg1) - np.exp(arg2) + 20. + np.e
    # 定义x和y的范围
    bounds = [(-5, 5), (-5, 5)]
    
    # 执行微分进化算法
    result = differential_evolution(ackley, bounds)
    result.x, result.fun
    
    

    img

    5.3.2 带约束

    一个约束

    img

    from scipy.optimize import differential_evolution,NonlinearConstraint
    import numpy as np
    
    # 定义目标函数
    def ackley(x):
        arg1 = -0.2 * np.sqrt(0.5 * (x[0] ** 2 + x[1] ** 2))
        arg2 = 0.5 * (np.cos(2. * np.pi * x[0]) + np.cos(2. * np.pi * x[1]))
        return -20. * np.exp(arg1) - np.exp(arg2) + 20. + np.e
    bounds = [(-5, 5), (-5, 5)]
    
    # 定义约束条件
    def constr_f(x):
        return np.array(x[0] + 2*x[1])
    
    # 实例约束条件
    # x1+x2<=4
    nlc = NonlinearConstraint(constr_f,-np.inf,4)
    
    result = differential_evolution(ackley, bounds,constraints=(nlc))
    print(result.x), print(result.fun)
    print(result)
    

    img

    两个约束

    from scipy.optimize import differential_evolution
    import numpy as np
    def ackley(x):
        arg1 = -0.2 * np.sqrt(0.5 * (x[0] ** 2 + x[1] ** 2))
        arg2 = 0.5 * (np.cos(2. * np.pi * x[0]) + np.cos(2. * np.pi * x[1]))
        return -20. * np.exp(arg1) - np.exp(arg2) + 20. + np.e
    bounds = [(-5, 5), (-5, 5)]
    
    # 设置两个约束条件
    def constr_f1(x):
        return np.array(x[0]*x[1])
    def constr_f2(x):
        return np.array(x[0] + 2*x[1])
    
    # x1+x2<=4
    nlc1 = NonlinearConstraint(constr_f1,-np.inf,4)
    nlc2 = NonlinearConstraint(constr_f2,-np.inf,4)
    
    result = differential_evolution(ackley, bounds,constraints=(nlc1,nlc2))
    print(result.x), print(result.fun)
    print(result)
    
    

    img

    6 求解二次规划

    使用工具:cvxopt.matrix、cvxopt.solvers

    6.1 例子1

    img

    from cvxopt import matrix
    from cvxopt import solvers
    
    P = matrix([[1.0,0.0],[0.0,0.0]])                   # P二次型系数矩阵
    q = matrix([3.0,4.0])                               # 一次项的系数矩阵
    G = matrix([[-1.0,0.0,-1.0,2.0,3.0],[0.0,-1.0,-3.0,5.0,4.0]])  # 约束条件系数矩阵
    h = matrix([0.0,0.0,-15.0,100.0,80.0])                         # 约束条件常数矩阵
    
    # 求解二次规划
    sol = solvers.qp(P,q,G,h)
    
    sol['x']
    sol['primal objective']
    
    print(sol['x'])
    

    img

    6.2 例子2

    img

    import pprint
    from cvxopt import matrix, solvers
    P = matrix([[4.0,1.0],[1.0,2.0]])           # 二次型系数矩阵
    q = matrix([1.0,1.0])                       # 一次项系数矩阵
    G = matrix([[-1.0,0.0],[0.0,-1.0]])         # 不等式约束系数矩阵
    h = matrix([0.0,0.0])                       # 不等式约束常数矩阵
    #A = matrix([1.0,1.0],(1,2))#原型为cvxopt.matrix(array,dims),等价于A = matrix([[1.0],[1.0]])
    A = matrix([[1.0],[1.0]])                   # 等式约束常数矩阵
    b = matrix([1.0])                           # 等式约束常数值
    result = solvers.qp(P,q,G,h,A,b)
     
    print('x\n',result['x'])
    

    img

    参考资料

    [1] Python–使用scipy求解带约束的最优化问题_爱吃汉堡薯条的博客-CSDN博客

    [2] 使用scipy.optimize.linprog求解线性规划问题_weixin_43638511的博客-CSDN博客

    [3] A Gentle Introduction to the BFGS Optimization Algorithm – Machine Learning Mastery

    [4] https://machinelearningmastery.com/dual-annealing-optimization-with-python/

    [5] https://machinelearningmastery.com/simulated-annealing-from-scratch-in-python/

    mastery.com/bfgs-optimization-in-python/)

    [4] https://machinelearningmastery.com/dual-annealing-optimization-with-python/

    [5] https://machinelearningmastery.com/simulated-annealing-from-scratch-in-python/

    [6] https://machinelearningmastery.com

    作者:斯凯利.瑞恩

    物联沃分享整理
    物联沃-IOTWORD物联网 » Python解决常见最优化问题的方法和策略

    发表回复